program test_nurbs_curve
use forcad_kinds, only: rk
use forcad_nurbs_curve, only: nurbs_curve, compute_Tgc, compute_dTgc
use forunittest, only: unit_tests
implicit none
integer, parameter :: NTESTS = 113
type(unit_tests) :: ut
real(rk), parameter :: TOL = 1.0e-5_rk
real(rk), parameter :: PI = acos(-1.0_rk)
type(nurbs_curve) :: nurbs, bsp
real(rk), allocatable :: Xc(:,:), Wc(:)
real(rk), allocatable :: Xg(:,:), Xgb(:,:)
real(rk) :: knot(6)
integer, allocatable :: elemConn(:,:)
real(rk), allocatable :: Tgc(:,:), dTgc(:,:), Tgcb(:,:), dTgcb(:,:), d2Tgc(:,:), d2Tgcb(:,:)
real(rk), allocatable :: Tgc1(:), dTgc1(:), Tgc1b(:), dTgc1b(:), d2Tgc1(:), d2Tgc1b(:)
real(rk) :: nearest_Xg(3), nearest_Xt, length, lengthb
integer :: i, id, ti
real(rk), allocatable :: Xt(:)
character(len=*), parameter :: fXc_nurbs = 'vtk/test_nurbs_curve_Xc.vtk'
character(len=*), parameter :: fXg_nurbs = 'vtk/test_nurbs_curve_Xg.vtk'
character(len=*), parameter :: fXth_nurbs = 'vtk/test_nurbs_curve_Xth.vtk'
character(len=*), parameter :: fIgs_nurbs = 'iges/test_nurbs_curve.iges'
character(len=*), parameter :: fXc_bsp = 'vtk/test_bsp_curve_Xc.vtk'
character(len=*), parameter :: fXg_bsp = 'vtk/test_bsp_curve_Xg.vtk'
character(len=*), parameter :: fXth_bsp = 'vtk/test_bsp_curve_Xth.vtk'
character(len=*), parameter :: fIgs_bsp = 'iges/test_bsp_curve.iges'
logical :: okXc_nurbs, okXg_nurbs, okXth_nurbs, okIges_nurbs
logical :: okXc_bsp, okXg_bsp, okXth_bsp, okIges_bsp
call ut%initialize(NTESTS)
ti = 1
! Initialize curve data
allocate(Xc(3, 3))
Xc(1,:) = [0.0_rk, 0.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, 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]
! 1) Set and create NURBS and B-spline curves
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=bsp%get_degree(), nc=bsp%get_nc(), Xc=bsp%get_Xc())
call nurbs%create(res=23)
call bsp%create(res=23)
call ut%test(ti)%check( &
name="set(): degree==2", &
res=nurbs%get_degree(), &
expected=2, &
msg="NURBS degree not 2", &
group="setup"); ti=ti+1
call ut%test(ti)%check( &
name="set(): degree==2 (B-spline)", &
res=bsp%get_degree(), &
expected=2, &
msg="B-spline degree not 2", &
group="setup"); ti=ti+1
call ut%test(ti)%check( &
name="set(): nc==3", &
res=nurbs%get_nc(), &
expected=3, &
msg="NURBS nc not 3", &
group="setup"); ti=ti+1
call ut%test(ti)%check( &
name="set(): nc==3 (B-spline)", &
res=bsp%get_nc(), &
expected=3, &
msg="B-spline nc not 3", &
group="setup"); ti=ti+1
call ut%test(ti)%check( &
name="set(): knot matches", &
res=nurbs%get_knot(), &
expected=knot, &
tol=TOL, &
msg="NURBS knot vector mismatch", &
group="setup"); ti=ti+1
call ut%test(ti)%check( &
name="set(): knot matches (B-spline)", &
res=bsp%get_knot(), &
expected=knot, &
tol=TOL, &
msg="B-spline knot vector mismatch", &
group="setup"); ti=ti+1
! 2) Export tests
call nurbs%export_Xc(fXc_nurbs)
call bsp%export_Xc(fXc_bsp)
call nurbs%export_Xg(fXg_nurbs)
call bsp%export_Xg(fXg_bsp)
call nurbs%export_Xth(fXth_nurbs)
call bsp%export_Xth(fXth_bsp)
call nurbs%export_iges(fIgs_nurbs)
call bsp%export_iges(fIgs_bsp)
inquire(file=fXc_nurbs, exist=okXc_nurbs)
inquire(file=fXg_nurbs, exist=okXg_nurbs)
inquire(file=fXth_nurbs, exist=okXth_nurbs)
inquire(file=fIgs_nurbs, exist=okIges_nurbs)
inquire(file=fXc_bsp, exist=okXc_bsp)
inquire(file=fXg_bsp, exist=okXg_bsp)
inquire(file=fXth_bsp, exist=okXth_bsp)
inquire(file=fIgs_bsp, exist=okIges_bsp)
call ut%test(ti)%check( &
name="export_Xc/Xg/Xth/IGES(): files exist (NURBS)", &
res=merge(1,0,okXc_nurbs .and. okXg_nurbs .and. okXth_nurbs .and. okIges_nurbs), &
expected=1, &
msg="NURBS export did not create all files", &
group="io"); ti=ti+1
call ut%test(ti)%check( &
name="export_Xc/Xg/Xth/IGES(): files exist (B-spline)", &
res=merge(1,0,okXc_bsp .and. okXg_bsp .and. okXth_bsp .and. okIges_bsp), &
expected=1, &
msg="B-spline export did not create all files", &
group="io"); ti=ti+1
! 3) Length tests
call nurbs%cmp_length(length)
call bsp%cmp_length(lengthb)
call ut%test(ti)%check( &
name="cmp_length(): length == 2 (NURBS)", &
res=length, &
expected=2.0_rk, &
tol=TOL, &
msg="NURBS arc length incorrect", &
group="geometry"); ti=ti+1
call ut%test(ti)%check( &
name="cmp_length(): length == 2 (B-spline)", &
res=lengthb, &
expected=2.0_rk, &
tol=TOL, &
msg="B-spline arc length incorrect", &
group="geometry"); ti=ti+1
! 4) Nearest point tests
call nurbs%nearest_point([0.0_rk, 0.0_rk, 0.5_rk], nearest_Xg, nearest_Xt, id)
call ut%test(ti)%check( &
name="nearest_point(): point matches [0,0,0] (NURBS)", &
res=nearest_Xg, &
expected=[0.0_rk, 0.0_rk, 0.0_rk], &
tol=TOL, &
msg="NURBS nearest point mismatch", &
group="nearest"); ti=ti+1
call ut%test(ti)%check( &
name="nearest_point(): param == 0 (NURBS)", &
res=nearest_Xt, &
expected=0.0_rk, &
tol=TOL, &
msg="NURBS nearest param incorrect", &
group="nearest"); ti=ti+1
call ut%test(ti)%check( &
name="nearest_point(): id == 1 (NURBS)", &
res=id, &
expected=1, &
msg="NURBS nearest id incorrect", &
group="nearest"); ti=ti+1
call bsp%nearest_point([0.0_rk, 0.0_rk, 0.5_rk], nearest_Xg, nearest_Xt, id)
call ut%test(ti)%check( &
name="nearest_point(): point matches [0,0,0] (B-spline)", &
res=nearest_Xg, &
expected=[0.0_rk, 0.0_rk, 0.0_rk], &
tol=TOL, &
msg="B-spline nearest point mismatch", &
group="nearest"); ti=ti+1
call ut%test(ti)%check( &
name="nearest_point(): param == 0 (B-spline)", &
res=nearest_Xt, &
expected=0.0_rk, &
tol=TOL, &
msg="B-spline nearest param incorrect", &
group="nearest"); ti=ti+1
call ut%test(ti)%check( &
name="nearest_point(): id == 1 (B-spline)", &
res=id, &
expected=1, &
msg="B-spline nearest id incorrect", &
group="nearest"); ti=ti+1
call nurbs%nearest_point2([0.0_rk, 0.0_rk, 0.5_rk], 1.0e-10_rk, 10, nearest_Xt, nearest_Xg)
call ut%test(ti)%check( &
name="nearest_point2(): point matches [0,0,0] (NURBS)", &
res=nearest_Xg, &
expected=[0.0_rk, 0.0_rk, 0.0_rk], &
tol=TOL, &
msg="NURBS nearest_point2 point mismatch", &
group="nearest"); ti=ti+1
call ut%test(ti)%check( &
name="nearest_point2(): param == 0 (NURBS)", &
res=nearest_Xt, &
expected=0.0_rk, &
tol=TOL, &
msg="NURBS nearest_point2 param incorrect", &
group="nearest"); ti=ti+1
call bsp%nearest_point2([0.0_rk, 0.0_rk, 0.5_rk], 1.0e-10_rk, 10, nearest_Xt, nearest_Xg)
call ut%test(ti)%check( &
name="nearest_point2(): point matches [0,0,0] (B-spline)", &
res=nearest_Xg, &
expected=[0.0_rk, 0.0_rk, 0.0_rk], &
tol=TOL, &
msg="B-spline nearest_point2 point mismatch", &
group="nearest"); ti=ti+1
call ut%test(ti)%check( &
name="nearest_point2(): param == 0 (B-spline)", &
res=nearest_Xt, &
expected=0.0_rk, &
tol=TOL, &
msg="B-spline nearest_point2 param incorrect", &
group="nearest"); ti=ti+1
! 5) Store initial geometry
Xg = nurbs%get_Xg()
Xgb = bsp%get_Xg()
! 6) Set with Xth_dir and continuity
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 nurbs%create(res=23)
call bsp%create(res=23)
call ut%test(ti)%check( &
name="set(Xth_dir): geometry preserved (NURBS)", &
res=nurbs%get_Xg(), &
expected=Xg, &
tol=TOL, &
msg="NURBS geometry changed after set(Xth_dir)", &
group="setup"); ti=ti+1
call ut%test(ti)%check( &
name="set(Xth_dir): geometry preserved (B-spline)", &
res=bsp%get_Xg(), &
expected=Xgb, &
tol=TOL, &
msg="B-spline geometry changed after set(Xth_dir)", &
group="setup"); ti=ti+1
! 7) Set with Xc and Wc
call nurbs%set(Xc, Wc)
call bsp%set(Xc)
call nurbs%create(res=23)
call bsp%create(res=23)
call ut%test(ti)%check( &
name="set(Xc,Wc): geometry preserved (NURBS)", &
res=nurbs%get_Xg(), &
expected=Xg, &
tol=TOL, &
msg="NURBS geometry changed after set(Xc,Wc)", &
group="setup"); ti=ti+1
call ut%test(ti)%check( &
name="set(Xc): geometry preserved (B-spline)", &
res=bsp%get_Xg(), &
expected=Xgb, &
tol=TOL, &
msg="B-spline geometry changed after set(Xc)", &
group="setup"); ti=ti+1
! 8) Create with explicit Xt
call nurbs%create(Xt=nurbs%get_Xt())
call bsp%create(Xt=bsp%get_Xt())
call ut%test(ti)%check( &
name="create(Xt): geometry preserved (NURBS)", &
res=nurbs%get_Xg(), &
expected=Xg, &
tol=TOL, &
msg="NURBS geometry changed after create(Xt)", &
group="sampling"); ti=ti+1
call ut%test(ti)%check( &
name="create(Xt): geometry preserved (B-spline)", &
res=bsp%get_Xg(), &
expected=Xgb, &
tol=TOL, &
msg="B-spline geometry changed after create(Xt)", &
group="sampling"); ti=ti+1
! 9) Export after create(Xt)
call nurbs%export_Xc(fXc_nurbs)
call bsp%export_Xc(fXc_bsp)
call nurbs%export_Xg(fXg_nurbs)
call bsp%export_Xg(fXg_bsp)
call nurbs%export_Xth(fXth_nurbs)
call bsp%export_Xth(fXth_bsp)
call nurbs%export_iges(fIgs_nurbs)
call bsp%export_iges(fIgs_bsp)
inquire(file=fXc_nurbs, exist=okXc_nurbs)
inquire(file=fXg_nurbs, exist=okXg_nurbs)
inquire(file=fXth_nurbs, exist=okXth_nurbs)
inquire(file=fIgs_nurbs, exist=okIges_nurbs)
inquire(file=fXc_bsp, exist=okXc_bsp)
inquire(file=fXg_bsp, exist=okXg_bsp)
inquire(file=fXth_bsp, exist=okXth_bsp)
inquire(file=fIgs_bsp, exist=okIges_bsp)
call ut%test(ti)%check( &
name="export after create(Xt): files exist (NURBS)", &
res=merge(1,0,okXc_nurbs .and. okXg_nurbs .and. okXth_nurbs .and. okIges_nurbs), &
expected=1, &
msg="NURBS export after create(Xt) did not create all files", &
group="io"); ti=ti+1
call ut%test(ti)%check( &
name="export after create(Xt): files exist (B-spline)", &
res=merge(1,0,okXc_bsp .and. okXg_bsp .and. okXth_bsp .and. okIges_bsp), &
expected=1, &
msg="B-spline export after create(Xt) did not create all files", &
group="io"); ti=ti+1
! 10) Getter tests
call ut%test(ti)%check( &
name="get_Xc(): matches input (NURBS)", &
res=nurbs%get_Xc(), &
expected=Xc, &
tol=TOL, &
msg="NURBS Xc mismatch", &
group="getters"); ti=ti+1
call ut%test(ti)%check( &
name="get_Xc(): matches input (B-spline)", &
res=bsp%get_Xc(), &
expected=Xc, &
tol=TOL, &
msg="B-spline Xc mismatch", &
group="getters"); ti=ti+1
call ut%test(ti)%check( &
name="get_Xc(1): matches input (NURBS)", &
res=nurbs%get_Xc(1), &
expected=Xc(1,:), &
tol=TOL, &
msg="NURBS Xc(1) mismatch", &
group="getters"); ti=ti+1
call ut%test(ti)%check( &
name="get_Xc(1): matches input (B-spline)", &
res=bsp%get_Xc(1), &
expected=Xc(1,:), &
tol=TOL, &
msg="B-spline Xc(1) mismatch", &
group="getters"); ti=ti+1
call ut%test(ti)%check( &
name="get_Xc(1,1): matches input (NURBS)", &
res=nurbs%get_Xc(1,1), &
expected=Xc(1,1), &
tol=TOL, &
msg="NURBS Xc(1,1) mismatch", &
group="getters"); ti=ti+1
call ut%test(ti)%check( &
name="get_Xc(1,1): matches input (B-spline)", &
res=bsp%get_Xc(1,1), &
expected=Xc(1,1), &
tol=TOL, &
msg="B-spline Xc(1,1) mismatch", &
group="getters"); ti=ti+1
call ut%test(ti)%check( &
name="get_Xg(): matches initial (NURBS)", &
res=nurbs%get_Xg(), &
expected=Xg, &
tol=TOL, &
msg="NURBS Xg mismatch", &
group="getters"); ti=ti+1
call ut%test(ti)%check( &
name="get_Xg(): matches initial (B-spline)", &
res=bsp%get_Xg(), &
expected=Xgb, &
tol=TOL, &
msg="B-spline Xg mismatch", &
group="getters"); ti=ti+1
call ut%test(ti)%check( &
name="get_Xg(1): matches initial (NURBS)", &
res=nurbs%get_Xg(1), &
expected=Xg(1,:), &
tol=TOL, &
msg="NURBS Xg(1) mismatch", &
group="getters"); ti=ti+1
call ut%test(ti)%check( &
name="get_Xg(1): matches initial (B-spline)", &
res=bsp%get_Xg(1), &
expected=Xgb(1,:), &
tol=TOL, &
msg="B-spline Xg(1) mismatch", &
group="getters"); ti=ti+1
call ut%test(ti)%check( &
name="get_Xg(1,1): matches initial (NURBS)", &
res=nurbs%get_Xg(1,1), &
expected=Xg(1,1), &
tol=TOL, &
msg="NURBS Xg(1,1) mismatch", &
group="getters"); ti=ti+1
call ut%test(ti)%check( &
name="get_Xg(1,1): matches initial (B-spline)", &
res=bsp%get_Xg(1,1), &
expected=Xgb(1,1), &
tol=TOL, &
msg="B-spline Xg(1,1) mismatch", &
group="getters"); ti=ti+1
call ut%test(ti)%check( &
name="get_Wc(): matches input (NURBS)", &
res=nurbs%get_Wc(), &
expected=Wc, &
tol=TOL, &
msg="NURBS Wc mismatch", &
group="getters"); ti=ti+1
call ut%test(ti)%check( &
name="get_Wc(1): matches input (NURBS)", &
res=nurbs%get_Wc(1), &
expected=Wc(1), &
tol=TOL, &
msg="NURBS Wc(1) mismatch", &
group="getters"); ti=ti+1
call ut%test(ti)%check( &
name="get_knot(): matches input (NURBS)", &
res=nurbs%get_knot(), &
expected=knot, &
tol=TOL, &
msg="NURBS knot mismatch", &
group="getters"); ti=ti+1
call ut%test(ti)%check( &
name="get_knot(): matches input (B-spline)", &
res=bsp%get_knot(), &
expected=knot, &
tol=TOL, &
msg="B-spline knot mismatch", &
group="getters"); ti=ti+1
call ut%test(ti)%check( &
name="get_knot(1): matches input (NURBS)", &
res=nurbs%get_knot(1), &
expected=knot(1), &
tol=TOL, &
msg="NURBS knot(1) mismatch", &
group="getters"); ti=ti+1
call ut%test(ti)%check( &
name="get_knot(1): matches input (B-spline)", &
res=bsp%get_knot(1), &
expected=knot(1), &
tol=TOL, &
msg="B-spline knot(1) mismatch", &
group="getters"); ti=ti+1
call ut%test(ti)%check( &
name="get_ng(): matches Xg size (NURBS)", &
res=nurbs%get_ng(), &
expected=size(Xg,1), &
msg="NURBS ng mismatch", &
group="getters"); ti=ti+1
call ut%test(ti)%check( &
name="get_ng(): matches Xgb size (B-spline)", &
res=bsp%get_ng(), &
expected=size(Xgb,1), &
msg="B-spline ng mismatch", &
group="getters"); ti=ti+1
call ut%test(ti)%check( &
name="get_degree(): matches 2 (NURBS)", &
res=nurbs%get_degree(), &
expected=2, &
msg="NURBS degree not 2", &
group="getters"); ti=ti+1
call ut%test(ti)%check( &
name="get_degree(): matches 2 (B-spline)", &
res=bsp%get_degree(), &
expected=2, &
msg="B-spline degree not 2", &
group="getters"); ti=ti+1
call ut%test(ti)%check( &
name="get_multiplicity(): matches [3,3] (NURBS)", &
res=nurbs%get_multiplicity(), &
expected=[3,3], &
msg="NURBS multiplicity not [3,3]", &
group="knot-ops"); ti=ti+1
call ut%test(ti)%check( &
name="get_multiplicity(): matches [3,3] (B-spline)", &
res=bsp%get_multiplicity(), &
expected=[3,3], &
msg="B-spline multiplicity not [3,3]", &
group="knot-ops"); ti=ti+1
call ut%test(ti)%check( &
name="get_continuity(): matches [-1,-1] (NURBS)", &
res=nurbs%get_continuity(), &
expected=[-1,-1], &
msg="NURBS continuity not [-1,-1]", &
group="knot-ops"); ti=ti+1
call ut%test(ti)%check( &
name="get_continuity(): matches [-1,-1] (B-spline)", &
res=bsp%get_continuity(), &
expected=[-1,-1], &
msg="B-spline continuity not [-1,-1]", &
group="knot-ops"); ti=ti+1
call nurbs%cmp_nc()
call bsp%cmp_nc()
call ut%test(ti)%check( &
name="cmp_nc(): matches nc (NURBS)", &
res=nurbs%get_nc(), &
expected=size(Xc,1), &
msg="NURBS cmp_nc mismatch", &
group="knot-ops"); ti=ti+1
call ut%test(ti)%check( &
name="cmp_nc(): matches nc (B-spline)", &
res=bsp%get_nc(), &
expected=size(Xc,1), &
msg="B-spline cmp_nc mismatch", &
group="knot-ops"); ti=ti+1
! 11) Element connectivity tests
elemConn = nurbs%cmp_elem_Xc_vis(2)
call nurbs%set_elem_Xc_vis(elemConn)
call ut%test(ti)%check( &
name="set/get elem_Xc_vis(p=2): equality (NURBS)", &
res=nurbs%get_elem_Xc_vis(), &
expected=elemConn, &
msg="NURBS elem_Xc_vis(p=2) mismatch", &
group="elements"); ti=ti+1
deallocate(elemConn)
elemConn = nurbs%cmp_elem_Xc_vis()
call nurbs%set_elem_Xc_vis(elemConn)
call ut%test(ti)%check( &
name="set/get elem_Xc_vis(): equality (NURBS)", &
res=nurbs%get_elem_Xc_vis(), &
expected=elemConn, &
msg="NURBS elem_Xc_vis mismatch", &
group="elements"); ti=ti+1
deallocate(elemConn)
elemConn = bsp%cmp_elem_Xc_vis(2)
call bsp%set_elem_Xc_vis(elemConn)
call ut%test(ti)%check( &
name="set/get elem_Xc_vis(p=2): equality (B-spline)", &
res=bsp%get_elem_Xc_vis(), &
expected=elemConn, &
msg="B-spline elem_Xc_vis(p=2) mismatch", &
group="elements"); ti=ti+1
deallocate(elemConn)
elemConn = bsp%cmp_elem_Xc_vis()
call bsp%set_elem_Xc_vis(elemConn)
call ut%test(ti)%check( &
name="set/get elem_Xc_vis(): equality (B-spline)", &
res=bsp%get_elem_Xc_vis(), &
expected=elemConn, &
msg="B-spline elem_Xc_vis mismatch", &
group="elements"); ti=ti+1
deallocate(elemConn)
elemConn = nurbs%cmp_elem_Xg_vis(2)
call nurbs%set_elem_Xg_vis(elemConn)
call ut%test(ti)%check( &
name="set/get elem_Xg_vis(p=2): equality (NURBS)", &
res=nurbs%get_elem_Xg_vis(), &
expected=elemConn, &
msg="NURBS elem_Xg_vis(p=2) mismatch", &
group="elements"); ti=ti+1
deallocate(elemConn)
elemConn = nurbs%cmp_elem_Xg_vis()
call nurbs%set_elem_Xg_vis(elemConn)
call ut%test(ti)%check( &
name="set/get elem_Xg_vis(): equality (NURBS)", &
res=nurbs%get_elem_Xg_vis(), &
expected=elemConn, &
msg="NURBS elem_Xg_vis mismatch", &
group="elements"); ti=ti+1
deallocate(elemConn)
elemConn = bsp%cmp_elem_Xg_vis(2)
call bsp%set_elem_Xg_vis(elemConn)
call ut%test(ti)%check( &
name="set/get elem_Xg_vis(p=2): equality (B-spline)", &
res=bsp%get_elem_Xg_vis(), &
expected=elemConn, &
msg="B-spline elem_Xg_vis(p=2) mismatch", &
group="elements"); ti=ti+1
deallocate(elemConn)
elemConn = bsp%cmp_elem_Xg_vis()
call bsp%set_elem_Xg_vis(elemConn)
call ut%test(ti)%check( &
name="set/get elem_Xg_vis(): equality (B-spline)", &
res=bsp%get_elem_Xg_vis(), &
expected=elemConn, &
msg="B-spline elem_Xg_vis mismatch", &
group="elements"); ti=ti+1
deallocate(elemConn)
elemConn = nurbs%cmp_elem()
call nurbs%set_elem(elemConn)
call ut%test(ti)%check( &
name="set/get elem: equality (NURBS)", &
res=nurbs%get_elem(), &
expected=elemConn, &
msg="NURBS elem mismatch", &
group="elements"); ti=ti+1
deallocate(elemConn)
elemConn = bsp%cmp_elem()
call bsp%set_elem(elemConn)
call ut%test(ti)%check( &
name="set/get elem: equality (B-spline)", &
res=bsp%get_elem(), &
expected=elemConn, &
msg="B-spline elem mismatch", &
group="elements"); ti=ti+1
deallocate(elemConn)
! 12) Modify and export
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 nurbs%export_Xc(fXc_nurbs)
call bsp%export_Xc(fXc_bsp)
call nurbs%export_Xg(fXg_nurbs)
call bsp%export_Xg(fXg_bsp)
call ut%test(ti)%check( &
name="modify_Xc/Wc: geometry preserved (NURBS)", &
res=nurbs%get_Xg(), &
expected=Xg, &
tol=TOL, &
msg="NURBS geometry changed after modify_Xc/Wc", &
group="modify"); ti=ti+1
call ut%test(ti)%check( &
name="modify_Xc: geometry preserved (B-spline)", &
res=bsp%get_Xg(), &
expected=Xgb, &
tol=TOL, &
msg="B-spline geometry changed after modify_Xc", &
group="modify"); ti=ti+1
! 13) Basis and derivative tests
call nurbs%basis(res=23, Tgc=Tgc)
call bsp%basis(res=23, Tgc=Tgcb)
call ut%test(ti)%check( &
name="basis(res=23): sum(N)=1 rows (NURBS)", &
res=all(abs(sum(Tgc,dim=2)-1.0_rk) <= TOL), &
expected=.true., &
msg="NURBS basis partition of unity violated", &
group="basis"); ti=ti+1
call ut%test(ti)%check( &
name="basis(res=23): sum(N)=1 rows (B-spline)", &
res=all(abs(sum(Tgcb,dim=2)-1.0_rk) <= TOL), &
expected=.true., &
msg="B-spline basis partition of unity violated", &
group="basis"); ti=ti+1
call nurbs%basis(Xt=0.0_rk, Tgc=Tgc1)
call bsp%basis(Xt=0.0_rk, Tgc=Tgc1b)
call ut%test(ti)%check( &
name="basis(Xt=0): sum(N)=1 (NURBS)", &
res=sum(Tgc1), &
expected=1.0_rk, &
tol=TOL, &
msg="NURBS basis sum not 1 at t=0", &
group="basis"); ti=ti+1
call ut%test(ti)%check( &
name="basis(Xt=0): sum(N)=1 (B-spline)", &
res=sum(Tgc1b), &
expected=1.0_rk, &
tol=TOL, &
msg="B-spline basis sum not 1 at t=0", &
group="basis"); ti=ti+1
allocate(Xt(23))
do i = 1, 23
Xt(i) = real(i-1,rk) / real(23-1,rk)
end do
call nurbs%basis(Xt=Xt, Tgc=Tgc)
call bsp%basis(Xt=Xt, Tgc=Tgcb)
call ut%test(ti)%check( &
name="basis(Xt vector): sum(N)=1 rows (NURBS)", &
res=all(abs(sum(Tgc,dim=2)-1.0_rk) <= TOL), &
expected=.true., &
msg="NURBS basis partition of unity violated (vector)", &
group="basis"); ti=ti+1
call ut%test(ti)%check( &
name="basis(Xt vector): sum(N)=1 rows (B-spline)", &
res=all(abs(sum(Tgcb,dim=2)-1.0_rk) <= TOL), &
expected=.true., &
msg="B-spline basis partition of unity violated (vector)", &
group="basis"); ti=ti+1
call nurbs%derivative(res=23, dTgc=dTgc, Tgc=Tgc)
call bsp%derivative(res=23, dTgc=dTgcb, Tgc=Tgcb)
call ut%test(ti)%check( &
name="derivative(res=23): sum(dN)=0 rows (NURBS)", &
res=all(abs(sum(dTgc,dim=2)) <= TOL), &
expected=.true., &
msg="NURBS derivative sum not zero", &
group="derivatives"); ti=ti+1
call ut%test(ti)%check( &
name="derivative(res=23): sum(dN)=0 rows (B-spline)", &
res=all(abs(sum(dTgcb,dim=2)) <= TOL), &
expected=.true., &
msg="B-spline derivative sum not zero", &
group="derivatives"); ti=ti+1
call nurbs%derivative(Xt=Xt, dTgc=dTgc, Tgc=Tgc)
call bsp%derivative(Xt=Xt, dTgc=dTgcb, Tgc=Tgcb)
call ut%test(ti)%check( &
name="derivative(Xt vector): sum(dN)=0 rows (NURBS)", &
res=all(abs(sum(dTgc,dim=2)) <= TOL), &
expected=.true., &
msg="NURBS derivative sum not zero (vector)", &
group="derivatives"); ti=ti+1
call ut%test(ti)%check( &
name="derivative(Xt vector): sum(dN)=0 rows (B-spline)", &
res=all(abs(sum(dTgcb,dim=2)) <= TOL), &
expected=.true., &
msg="B-spline derivative sum not zero (vector)", &
group="derivatives"); ti=ti+1
call nurbs%derivative(Xt=0.0_rk, dTgc=dTgc1, Tgc=Tgc1)
call bsp%derivative(Xt=0.0_rk, dTgc=dTgc1b, Tgc=Tgc1b)
call ut%test(ti)%check( &
name="derivative(Xt=0): sum(dN)=0 (NURBS)", &
res=sum(dTgc1), &
expected=0.0_rk, &
tol=TOL, &
msg="NURBS derivative sum not zero at t=0", &
group="derivatives"); ti=ti+1
call ut%test(ti)%check( &
name="derivative(Xt=0): sum(dN)=0 (B-spline)", &
res=sum(dTgc1b), &
expected=0.0_rk, &
tol=TOL, &
msg="B-spline derivative sum not zero at t=0", &
group="derivatives"); ti=ti+1
call nurbs%derivative(Xt=0.0_rk, dTgc=dTgc1, Tgc=Tgc1, elem=[1,2,3])
call bsp%derivative(Xt=0.0_rk, dTgc=dTgc1b, Tgc=Tgc1b, elem=[1,2,3])
call ut%test(ti)%check( &
name="derivative(Xt=0, elem): sum(dN)=0 (NURBS)", &
res=sum(dTgc1), &
expected=0.0_rk, &
tol=TOL, &
msg="NURBS derivative sum not zero (elem)", &
group="derivatives"); ti=ti+1
call ut%test(ti)%check( &
name="derivative(Xt=0, elem): sum(dN)=0 (B-spline)", &
res=sum(dTgc1b), &
expected=0.0_rk, &
tol=TOL, &
msg="B-spline derivative sum not zero (elem)", &
group="derivatives"); ti=ti+1
call nurbs%derivative2(res=23, d2Tgc=d2Tgc, dTgc=dTgc, Tgc=Tgc)
call bsp%derivative2(res=23, d2Tgc=d2Tgcb, dTgc=dTgcb, Tgc=Tgcb)
call ut%test(ti)%check( &
name="derivative2(res=23): sum(d2N)=0 rows (NURBS)", &
res=all(abs(sum(d2Tgc,dim=2)) <= TOL), &
expected=.true., &
msg="NURBS second derivative sum not zero", &
group="derivatives"); ti=ti+1
call ut%test(ti)%check( &
name="derivative2(res=23): sum(d2N)=0 rows (B-spline)", &
res=all(abs(sum(d2Tgcb,dim=2)) <= TOL), &
expected=.true., &
msg="B-spline second derivative sum not zero", &
group="derivatives"); ti=ti+1
call nurbs%derivative2(Xt=Xt, d2Tgc=d2Tgc, dTgc=dTgc, Tgc=Tgc)
call bsp%derivative2(Xt=Xt, d2Tgc=d2Tgcb, dTgc=dTgcb, Tgc=Tgcb)
call ut%test(ti)%check( &
name="derivative2(Xt vector): sum(d2N)=0 rows (NURBS)", &
res=all(abs(sum(d2Tgc,dim=2)) <= TOL), &
expected=.true., &
msg="NURBS second derivative sum not zero (vector)", &
group="derivatives"); ti=ti+1
call ut%test(ti)%check( &
name="derivative2(Xt vector): sum(d2N)=0 rows (B-spline)", &
res=all(abs(sum(d2Tgcb,dim=2)) <= TOL), &
expected=.true., &
msg="B-spline second derivative sum not zero (vector)", &
group="derivatives"); ti=ti+1
call nurbs%derivative2(Xt=0.0_rk, d2Tgc=d2Tgc1, dTgc=dTgc1, Tgc=Tgc1)
call bsp%derivative2(Xt=0.0_rk, d2Tgc=d2Tgc1b, dTgc=dTgc1b, Tgc=Tgc1b)
call ut%test(ti)%check( &
name="derivative2(Xt=0): sum(d2N)=0 (NURBS)", &
res=sum(d2Tgc1), &
expected=0.0_rk, &
tol=TOL, &
msg="NURBS second derivative sum not zero at t=0", &
group="derivatives"); ti=ti+1
call ut%test(ti)%check( &
name="derivative2(Xt=0): sum(d2N)=0 (B-spline)", &
res=sum(d2Tgc1b), &
expected=0.0_rk, &
tol=TOL, &
msg="B-spline second derivative sum not zero at t=0", &
group="derivatives"); ti=ti+1
! 14) Rotation tests (Xc)
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%test(ti)%check( &
name="rotate_Xc(X-axis): geometry preserved (NURBS)", &
res=nurbs%get_Xc(), &
expected=Xc, &
tol=TOL, &
msg="NURBS Xc changed after X-axis rotation", &
group="transform"); ti=ti+1
call ut%test(ti)%check( &
name="rotate_Xc(X-axis): geometry preserved (B-spline)", &
res=bsp%get_Xc(), &
expected=Xc, &
tol=TOL, &
msg="B-spline Xc changed after X-axis rotation", &
group="transform"); ti=ti+1
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%test(ti)%check( &
name="rotate_Xc(Y-axis): geometry preserved (NURBS)", &
res=nurbs%get_Xc(), &
expected=Xc, &
tol=TOL, &
msg="NURBS Xc changed after Y-axis rotation", &
group="transform"); ti=ti+1
call ut%test(ti)%check( &
name="rotate_Xc(Y-axis): geometry preserved (B-spline)", &
res=bsp%get_Xc(), &
expected=Xc, &
tol=TOL, &
msg="B-spline Xc changed after Y-axis rotation", &
group="transform"); ti=ti+1
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%test(ti)%check( &
name="rotate_Xc(Z-axis): geometry preserved (NURBS)", &
res=nurbs%get_Xc(), &
expected=Xc, &
tol=TOL, &
msg="NURBS Xc changed after Z-axis rotation", &
group="transform"); ti=ti+1
call ut%test(ti)%check( &
name="rotate_Xc(Z-axis): geometry preserved (B-spline)", &
res=bsp%get_Xc(), &
expected=Xc, &
tol=TOL, &
msg="B-spline Xc changed after Z-axis rotation", &
group="transform"); ti=ti+1
! 15) Rotation tests (Xg)
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%test(ti)%check( &
name="rotate_Xg(X-axis): geometry preserved (NURBS)", &
res=nurbs%get_Xg(), &
expected=Xg, &
tol=TOL, &
msg="NURBS Xg changed after X-axis rotation", &
group="transform"); ti=ti+1
call ut%test(ti)%check( &
name="rotate_Xg(X-axis): geometry preserved (B-spline)", &
res=bsp%get_Xg(), &
expected=Xgb, &
tol=TOL, &
msg="B-spline Xg changed after X-axis rotation", &
group="transform"); ti=ti+1
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%test(ti)%check( &
name="rotate_Xg(Y-axis): geometry preserved (NURBS)", &
res=nurbs%get_Xg(), &
expected=Xg, &
tol=TOL, &
msg="NURBS Xg changed after Y-axis rotation", &
group="transform"); ti=ti+1
call ut%test(ti)%check( &
name="rotate_Xg(Y-axis): geometry preserved (B-spline)", &
res=bsp%get_Xg(), &
expected=Xgb, &
tol=TOL, &
msg="B-spline Xg changed after Y-axis rotation", &
group="transform"); ti=ti+1
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%test(ti)%check( &
name="rotate_Xg(Z-axis): geometry preserved (NURBS)", &
res=nurbs%get_Xg(), &
expected=Xg, &
tol=TOL, &
msg="NURBS Xg changed after Z-axis rotation", &
group="transform"); ti=ti+1
call ut%test(ti)%check( &
name="rotate_Xg(Z-axis): geometry preserved (B-spline)", &
res=bsp%get_Xg(), &
expected=Xgb, &
tol=TOL, &
msg="B-spline Xg changed after Z-axis rotation", &
group="transform"); ti=ti+1
! 16) Translation tests (Xc)
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%test(ti)%check( &
name="translate_Xc(): geometry preserved (NURBS)", &
res=nurbs%get_Xc(), &
expected=Xc, &
tol=TOL, &
msg="NURBS Xc changed after translation", &
group="transform"); ti=ti+1
call ut%test(ti)%check( &
name="translate_Xc(): geometry preserved (B-spline)", &
res=bsp%get_Xc(), &
expected=Xc, &
tol=TOL, &
msg="B-spline Xc changed after translation", &
group="transform"); ti=ti+1
! 17) Translation tests (Xg)
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%test(ti)%check( &
name="translate_Xg(): geometry preserved (NURBS)", &
res=nurbs%get_Xg(), &
expected=Xg, &
tol=TOL, &
msg="NURBS Xg changed after translation", &
group="transform"); ti=ti+1
call ut%test(ti)%check( &
name="translate_Xg(): geometry preserved (B-spline)", &
res=bsp%get_Xg(), &
expected=Xgb, &
tol=TOL, &
msg="B-spline Xg changed after translation", &
group="transform"); ti=ti+1
! 18) Knot insertion
call nurbs%insert_knots([0.25_rk, 0.75_rk], [2,1])
call bsp%insert_knots([0.25_rk, 0.75_rk], [2,1])
call nurbs%create()
call bsp%create()
call ut%test(ti)%check( &
name="insert_knots(): geometry preserved (NURBS)", &
res=nurbs%get_Xg(), &
expected=Xg, &
tol=TOL, &
msg="NURBS geometry changed after knot insertion", &
group="knot-ops"); ti=ti+1
call ut%test(ti)%check( &
name="insert_knots(): geometry preserved (B-spline)", &
res=bsp%get_Xg(), &
expected=Xgb, &
tol=TOL, &
msg="B-spline geometry changed after knot insertion", &
group="knot-ops"); ti=ti+1
! 19) Degree elevation
call nurbs%elevate_degree(2)
call bsp%elevate_degree(2)
call nurbs%create()
call bsp%create()
call ut%test(ti)%check( &
name="elevate_degree(): geometry preserved (NURBS)", &
res=nurbs%get_Xg(), &
expected=Xg, &
tol=TOL, &
msg="NURBS geometry changed after degree elevation", &
group="knot-ops"); ti=ti+1
call ut%test(ti)%check( &
name="elevate_degree(): geometry preserved (B-spline)", &
res=bsp%get_Xg(), &
expected=Xgb, &
tol=TOL, &
msg="B-spline geometry changed after degree elevation", &
group="knot-ops"); ti=ti+1
! 20) Knot removal
call nurbs%remove_knots([0.25_rk, 0.75_rk], [2,1])
call bsp%remove_knots([0.25_rk, 0.75_rk], [2,1])
call nurbs%create()
call bsp%create()
call ut%test(ti)%check( &
name="remove_knots(): geometry preserved (NURBS)", &
res=nurbs%get_Xg(), &
expected=Xg, &
tol=TOL, &
msg="NURBS geometry changed after knot removal", &
group="knot-ops"); ti=ti+1
call ut%test(ti)%check( &
name="remove_knots(): geometry preserved (B-spline)", &
res=bsp%get_Xg(), &
expected=Xgb, &
tol=TOL, &
msg="B-spline geometry changed after knot removal", &
group="knot-ops"); ti=ti+1
! 21) Shape functions (circle, half-circle)
call nurbs%set_circle([0.0_rk, 0.0_rk, 0.0_rk], 1.0_rk)
call nurbs%create(res=23)
call nurbs%nearest_point([1.0_rk, 0.0_rk, 0.0_rk], nearest_Xg, nearest_Xt, id)
call ut%test(ti)%check( &
name="set_circle(): nearest point param (NURBS)", &
res=nearest_Xt, &
expected=0.0_rk, &
tol=TOL, &
msg="NURBS circle nearest point param incorrect", &
group="shapes"); ti=ti+1
call nurbs%set_half_circle([0.0_rk, 0.0_rk, 0.0_rk], 1.0_rk)
call nurbs%create(res=23)
call nurbs%nearest_point([0.0_rk, 1.0_rk, 0.0_rk], nearest_Xg, nearest_Xt, id)
call ut%test(ti)%check( &
name="set_half_circle(): nearest point param (NURBS)", &
res=nearest_Xt, &
expected=0.5_rk, &
tol=TOL, &
msg="NURBS half-circle nearest point param incorrect", &
group="shapes"); ti=ti+1
! 22) Least-squares B-spline fitting
block
type(nurbs_curve) :: bsp_fit
integer :: n
real(rk), allocatable :: Xt_fit(:), Xdata(:,:), Xg_eval(:,:)
real(rk) :: err1, err2, err3, rms
n = 42
allocate(Xt_fit(n), Xdata(n, 3))
do i = 1, n
Xt_fit(i) = real(i-1, rk) / real(n-1, rk)
Xdata(i,1) = Xt_fit(i)
Xdata(i,2) = 0.3_rk * sin(4.0_rk * PI * Xt_fit(i))
Xdata(i,3) = 0.3_rk * cos(4.0_rk * PI * Xt_fit(i))
end do
call bsp_fit%set( &
degree=5, &
Xth_dir=[(real(i-1, rk)/10.0_rk, i=1,11)], &
continuity=[-1, 1, 1, 1, 1, 1, 1, 1, 1, 1, -1])
call bsp_fit%lsq_fit_bspline(Xt_fit, Xdata, n)
call bsp_fit%create(res=n)
Xg_eval = bsp_fit%get_Xg()
err1 = norm2(Xg_eval(:,1) - Xdata(:,1)) / norm2(Xdata(:,1))
err2 = norm2(Xg_eval(:,2) - Xdata(:,2)) / norm2(Xdata(:,2))
err3 = norm2(Xg_eval(:,3) - Xdata(:,3)) / norm2(Xdata(:,3))
rms = sqrt((err1**2 + err2**2 + err3**2) / 3.0_rk)
call ut%test(ti)%check( &
name="lsq_fit_bspline(): RMS error", &
res=rms, &
expected=0.0_rk, &
tol=1.0e-6_rk, &
msg="B-spline least-squares fit RMS error too high", &
group="basis"); ti=ti+1
call bsp_fit%finalize()
deallocate(Xt_fit, Xdata, Xg_eval)
end block
! Finalize
call nurbs%finalize()
call bsp%finalize()
deallocate(Xc, Wc, Xg, Xgb, Xt)
if (allocated(Tgc)) deallocate(Tgc)
if (allocated(Tgcb)) deallocate(Tgcb)
if (allocated(dTgc)) deallocate(dTgc)
if (allocated(dTgcb)) deallocate(dTgcb)
if (allocated(d2Tgc)) deallocate(d2Tgc)
if (allocated(d2Tgcb)) deallocate(d2Tgcb)
if (allocated(Tgc1)) deallocate(Tgc1)
if (allocated(Tgc1b)) deallocate(Tgc1b)
if (allocated(dTgc1)) deallocate(dTgc1)
if (allocated(dTgc1b)) deallocate(dTgc1b)
if (allocated(d2Tgc1)) deallocate(d2Tgc1)
if (allocated(d2Tgc1b)) deallocate(d2Tgc1b)
call ut%summary(verbose=3, required_score=100.0)
end program test_nurbs_curve