!> Author: Seyed Ali Ghasemi !> License: BSD 3-Clause !> Unit test program for `forcad_utils`. !> using ForUnitTest: https://github.com/gha3mi/forunittest program test_forcad_utils use forcad_kinds, only: rk use forcad_utils, only: basis_bspline, basis_bspline_der, basis_bspline_2der, & basis_bernstein, compute_multiplicity, ndgrid, dyad, kron, unique, findspan, & compute_knot_vector, insert_knot_A_5_1, remove_knots_A_5_8, elevate_degree_A_5_9, & hexahedron_Xc, tetragon_Xc, elemConn_C0, elemConn_Cn, rotation, det, inv, gauss_leg, & export_vtk_legacy, solve, repelem, linspace, eye, kron_eye use forunittest, only: unit_tests implicit none type(unit_tests) :: ut real(rk) :: Xt integer :: nc, degree real(rk) :: knot(7) real(rk) :: B4(4), dB(4), d2B(4), A4(2,2) real(rk) :: B_ref(4), dB_ref(4), d2B_ref(4) real(rk) :: u(2), v(2), w(4) real(rk), allocatable :: u2(:), v2(:), w2(:) real(rk) :: A2x2(2,2), Bk(4,2) real(rk), allocatable :: A(:), vec(:), M(:,:) real(rk), allocatable :: X1(:), X2(:), X3(:) real(rk), allocatable :: Xt2(:,:), Xt3(:,:) real(rk), allocatable :: R(:,:), R_expected(:,:) real(rk), allocatable :: knot_in(:), knot_out(:), Pw(:,:), Qw(:,:) real(rk), allocatable :: Xksi(:,:), Wksi(:) real(rk), allocatable :: K3(:), K2(:), K1(:), out(:) integer, allocatable :: conn1D(:,:), conn2D(:,:), conn3D(:,:) integer :: nq, p, rr, s, k character(len=*), parameter :: vtk_file = "vtk/test_output.vtk" real(rk), allocatable :: A2(:,:), A_inv(:,:) ! Initialize unit tests call ut%initialize(n=57) ! ---------------------------- ! Test: basis_bspline ! ---------------------------- degree = 2 nc = 4 knot = [0.0_rk, 0.0_rk, 0.0_rk, 0.5_rk, 1.0_rk, 1.0_rk, 1.0_rk] Xt = 0.5_rk B4 = basis_bspline(Xt, knot, nc, degree) B_ref = [0.0_rk, 0.5_rk, 0.5_rk, 0.0_rk] call ut%test(1)%check( & name = "basis_bspline", & res = B4, & expected = B_ref, & msg = "Partition of unity and shape check", & group = "basis") call ut%test(2)%check( & name = "basis_sum", & res = sum(B4), & expected = 1.0_rk, & msg = "Partition of unity", & group = "basis") call basis_bspline_der(Xt, knot, nc, degree, dB, B4) dB_ref = [0.0_rk, -2.0_rk, 2.0_rk, 0.0_rk] call ut%test(3)%check( & name = "basis_bspline_der", & res = dB, & expected = dB_ref, & msg = "1st derivative shape check", & group = "basis_der") call basis_bspline_der(Xt, knot, nc, degree, dB) call ut%test(4)%check( & name = "basis_bspline_der_B", & res = dB, & expected = dB_ref, & msg = "1st derivative alternate", & group = "basis_der") call basis_bspline_2der(Xt, knot, nc, degree, d2B, dB, B4) d2B_ref = [0.0_rk, 4.0_rk, -12.0_rk, 8.0_rk] call ut%test(5)%check( & name = "basis_bspline_2der_A", & res = d2B, & expected = d2B_ref, & msg = "2nd derivative shape check A", & group = "basis_2der") call basis_bspline_2der(Xt, knot, nc, degree, d2B, dB) call ut%test(6)%check( & name = "basis_bspline_2der_B", & res = d2B, & expected = d2B_ref, & msg = "2nd derivative shape check B", & group = "basis_2der") call basis_bspline_2der(Xt, knot, nc, degree, d2B) call ut%test(7)%check( & name = "basis_bspline_2der_C", & res = d2B, & expected = d2B_ref, & msg = "2nd derivative shape check C", & group = "basis_2der") Xt = -0.1_rk B4 = basis_bspline(Xt, knot, nc, degree) call ut%test(8)%check( & name = "basis_out_of_bounds_low", & res = sum(B4), & expected = 0.0_rk, & msg = "Out-of-bounds left", & group = "bounds") Xt = 1.1_rk B4 = basis_bspline(Xt, knot, nc, degree) call ut%test(9)%check( & name = "basis_out_of_bounds_high", & res = sum(B4), & expected = 0.0_rk, & msg = "Out-of-bounds right", & group = "bounds") ! ---------------------------- ! Test: basis_bernstein ! ---------------------------- call ut%test(10)%check( & name = "basis_bernstein_sum", & res = sum(basis_bernstein(0.3_rk, 3)), & expected = 1.0_rk, & msg = "Bernstein basis partition of unity", & group = "bernstein") ! ---------------------------- ! Test: compute_multiplicity ! ---------------------------- call ut%test(11)%check( & name = "compute_multiplicity1", & res = compute_multiplicity([0.0_rk, 0.0_rk, 1.0_rk, 1.0_rk]), & expected = [2, 2], & msg = "Multiplicity vector", & group = "multiplicity") call ut%test(12)%check( & name = "compute_multiplicity2", & res = compute_multiplicity([0.0_rk, 0.0_rk, 0.0_rk, 1.0_rk, 1.0_rk, 1.0_rk], 1.0_rk), & expected = 3, & msg = "Multiplicity at value", & group = "multiplicity") ! ---------------------------- ! Test: ndgrid ! ---------------------------- X1 = [1.0_rk, 2.0_rk] X2 = [10.0_rk, 20.0_rk] X3 = [100.0_rk] call ndgrid(X1, X2, Xt2) call ndgrid(X1, X2, X3, Xt3) call ut%test(13)%check( & name = "ndgrid2_shape", & res = shape(Xt2), & expected = [4,2], & msg = "ndgrid2 shape", & group = "ndgrid") call ut%test(14)%check( & name = "ndgrid3_value", & res = Xt3(:,3), & expected = [100.0_rk, 100.0_rk, 100.0_rk, 100.0_rk], & msg = "ndgrid3 constant Z", & group = "ndgrid") ! ---------------------------- ! Test: dyad ! ---------------------------- A = [1.0_rk, 2.0_rk] vec = [3.0_rk, 4.0_rk, 5.0_rk] M = dyad(A, vec) call ut%test(15)%check( & name = "dyad_check", & res = M, & expected = reshape([3.0_rk, 6.0_rk, 4.0_rk, 8.0_rk, 5.0_rk, 10.0_rk], [2,3]), & msg = "Outer product a .dyad. b", & group = "dyad") ! ---------------------------- ! Test: kron ! ---------------------------- u = [1.0_rk, 2.0_rk] v = [3.0_rk, 4.0_rk] w = kron(u, v) call ut%test(16)%check( & name = "kron_vector", & res = w, & expected = [3.0_rk, 4.0_rk, 6.0_rk, 8.0_rk], & msg = "u .kron. v", & group = "kron") A2x2 = reshape([1.0_rk, 2.0_rk, 3.0_rk, 4.0_rk], [2,2]) Bk = kron(u, A2x2) call ut%test(17)%check( & name = "kron_matrix", & res = Bk, & expected = reshape([1.0_rk, 2.0_rk, 2.0_rk, 4.0_rk, 3.0_rk, 4.0_rk, 6.0_rk, 8.0_rk], [4,2]), & msg = "u .kron. A", & group = "kron") ! ---------------------------- ! Test: unique ! ---------------------------- call ut%test(18)%check( & name = "unique_integer", & res = unique([1,2,2,3,1,4]), & expected = [1,2,3,4], & msg = "Unique integers", & group = "unique") call ut%test(19)%check( & name = "unique_real", & res = unique([1.0_rk, 1.0_rk + 1e-16_rk, 2.0_rk, 1.0_rk, 3.0_rk]), & expected = [1.0_rk, 2.0_rk, 3.0_rk], & msg = "Unique real with tolerance", & group = "unique") ! ---------------------------- ! Test: findspan ! ---------------------------- call ut%test(20)%check( & name = "findspan_middle", & res = findspan(4, 2, 0.5_rk, knot), & expected = 3, & msg = "Find span index at Xt=0.5", & group = "findspan") ! ---------------------------- ! Test: compute_knot_vector ! ---------------------------- call ut%test(21)%check( & name = "compute_knot_vector", & res = compute_knot_vector([0.0_rk, 1.0_rk, 2.0_rk], 2, [-1, 1, -1]), & expected = [0.0_rk, 0.0_rk, 0.0_rk, 1.0_rk, 2.0_rk, 2.0_rk, 2.0_rk], & msg = "Knot vector construction", & group = "knot") ! ---------------------------- ! Test: repelem ! ---------------------------- call ut%test(22)%check( & name = "repelem", & res = repelem([1.0_rk, 2.0_rk, 3.0_rk], [2, 1, 3]), & expected = [1.0_rk, 1.0_rk, 2.0_rk, 3.0_rk, 3.0_rk, 3.0_rk], & msg = "Repeat vector elements", & group = "repelem") ! ---------------------------- ! Test: hexahedron_Xc ! ---------------------------- call ut%test(23)%check( & name = "hexahedron_Xc_shape", & res = shape(hexahedron_Xc([1.0_rk, 1.0_rk, 1.0_rk], [2,2,2])), & expected = [8,3], & msg = "3D grid shape", & group = "geometry") ! ---------------------------- ! Test: tetragon_Xc ! ---------------------------- call ut%test(24)%check( & name = "tetragon_Xc_shape", & res = shape(tetragon_Xc([1.0_rk, 1.0_rk], [2,2])), & expected = [4,3], & msg = "2D grid shape", & group = "geometry") ! ---------------------------- ! Test: rotation ! ---------------------------- call ut%test(25)%check( & name = "rotation_identity", & res = rotation(0.0_rk, 0.0_rk, 0.0_rk), & expected = reshape([1.0_rk,0.0_rk,0.0_rk, 0.0_rk,1.0_rk,0.0_rk, 0.0_rk,0.0_rk,1.0_rk], [3,3]), & msg = "Rotation(0,0,0) = I", & group = "rotation") ! ---------------------------- ! Test: det ! ---------------------------- call ut%test(26)%check( & name = "det_2x2", & res = det(reshape([1.0_rk, 2.0_rk, 3.0_rk, 4.0_rk], [2,2])), & expected = -2.0_rk, & msg = "Determinant 2x2", & group = "matrix") ! ---------------------------- ! Test: inv ! ---------------------------- call ut%test(27)%check( & name = "inv_2x2", & res = inv(reshape([4.0_rk, 7.0_rk, 2.0_rk, 6.0_rk], [2,2])), & expected = reshape([0.6_rk, -0.7_rk, -0.2_rk, 0.4_rk], [2,2]), & msg = "Inverse of 2x2", & group = "matrix") ! ---------------------------- ! Test: solve ! ---------------------------- A4 = reshape([4.0_rk, 1.0_rk, 1.0_rk, 3.0_rk], [2,2]) R_expected = reshape([1.0_rk/11.0_rk, 7.0_rk/11.0_rk], [2,1]) R = solve(A4, reshape([1.0_rk, 2.0_rk], [2,1])) call ut%test(28)%check( & name = "solve_linear_system", & res = R, & expected = R_expected, & msg = "Solving A.X = B", & group = "matrix") ! ---------------------------- ! Test: insert_knot_A_5_1 ! ---------------------------- p = 2 rr = 1 s = 1 k = 3 allocate(knot_in(0:6)) knot_in = [0.0_rk, 0.0_rk, 0.0_rk, 0.5_rk, 1.0_rk, 1.0_rk, 1.0_rk] allocate(Pw(0:2,1:2)) Pw(0,:) = [0.0_rk, 0.0_rk] Pw(1,:) = [0.5_rk, 0.5_rk] Pw(2,:) = [1.0_rk, 1.0_rk] call insert_knot_A_5_1(p, knot_in, Pw, 0.5_rk, k, s, rr, nq, knot_out, Qw) call ut%test(29)%check( & name = "insert_knot_A_5_1_nq", & res = nq, & expected = 3, & msg = "Inserted knot, new number of control points", & group = "insert_knot") ! ---------------------------- ! Test: elevate_degree_A_5_9 ! ---------------------------- deallocate(knot_in, knot_out, Pw, Qw) knot_in = [0.0_rk, 0.0_rk, 1.0_rk, 1.0_rk] allocate(Pw(1:2,1:2)) Pw(1,:) = [0.0_rk, 0.0_rk] Pw(2,:) = [0.5_rk, 0.5_rk] call elevate_degree_A_5_9(t=1, knot=knot_in, degree=1, Xcw=Pw, nc_new=nc, knot_new=knot_out, Xcw_new=Qw) call ut%test(30)%check( & name = "elevate_degree_nc", & res = nc, & expected = 3, & msg = "New number of control points after elevation", & group = "elevate_degree") ! ---------------------------- ! Test: gauss_legendre_1D ! ---------------------------- call gauss_leg([0.0_rk, 1.0_rk], 2, Xksi=vec, Wksi=A) call ut%test(31)%check( & name = "gauss_legendre_1D_pts", & res = size(vec), & expected = 3, & msg = "Number of Gauss points (1D)", & group = "quadrature") ! ---------------------------- ! Test: gauss_legendre_2D ! ---------------------------- call gauss_leg([0.0_rk, 1.0_rk], [0.0_rk, 1.0_rk], [2, 2], Xksi=Xksi, Wksi=Wksi) call ut%test(32)%check( & name = "gauss_legendre_2D_shape", & res = shape(Xksi), & expected = [9,2], & msg = "Gauss points shape (2D)", & group = "quadrature") ! ---------------------------- ! Test: gauss_legendre_3D ! ---------------------------- call gauss_leg([0.0_rk, 1.0_rk], [0.0_rk, 1.0_rk], [0.0_rk, 1.0_rk], [1, 1, 1], Xksi=Xksi, Wksi=Wksi) call ut%test(33)%check( & name = "gauss_legendre_3D_size", & res = size(Xksi,1), & expected = 8, & msg = "Number of Gauss points (3D)", & group = "quadrature") ! ---------------------------- ! Test: export_vtk_legacy ! ---------------------------- call export_vtk_legacy(filename=vtk_file, points=reshape([0.0_rk, 0.0_rk, 0.0_rk, 1.0_rk, 0.0_rk, 0.0_rk], [2,3]), & elemConn=reshape([1,2], [1,2]), vtkCellType=3, encoding="binary") call ut%test(34)%check( & name = "export_vtk_legacy", & res = .true., & expected = .true., & msg = "Export to VTK (not crashing)", & group = "export") call export_vtk_legacy(filename=vtk_file, points=reshape([0.0_rk, 0.0_rk, 0.0_rk, 1.0_rk, 0.0_rk, 0.0_rk], [2,3]), & elemConn=reshape([1,2], [1,2]), vtkCellType=3, encoding="ascii") call ut%test(35)%check( & name = "export_vtk_legacy", & res = .true., & expected = .true., & msg = "Export to VTK (not crashing)", & group = "export") ! ---------------------------- ! Test: linspace ! ---------------------------- call ut%test(36)%check( & name = "linspace_uniform", & res = linspace(0.0_rk, 1.0_rk, 5), & expected = [0.0_rk, 0.25_rk, 0.5_rk, 0.75_rk, 1.0_rk], & msg = "Uniform linspace from 0 to 1", & group = "linspace") ! ---------------------------- ! Test: kron3 ! ---------------------------- K1 = [1.0_rk, 2.0_rk] K2 = [3.0_rk] K3 = [4.0_rk, 5.0_rk] out = kron(K1, K2, K3) call ut%test(37)%check( & name = "kron3_product", & res = out, & expected = [12.0_rk, 15.0_rk, 24.0_rk, 30.0_rk], & msg = "kron3 result values", & group = "kron") ! ---------------------------- ! Test: elemConn_C0 ! ---------------------------- conn1D = elemConn_C0(5, 2) call ut%test(38)%check( & name = "elemConn_C0_L", & res = conn1D, & expected = reshape([1,3,2,4,3,5], [2,3]), & msg = "Linear C0 connectivity", & group = "connectivity") ! ---------------------------- ! Test: elemConn_Cn (L) ! ---------------------------- call elemConn_Cn(5, 2, [0.0_rk, 0.5_rk, 1.0_rk], [3,1,3], conn1D) call ut%test(39)%check( & name = "elemConn_Cn_L", & res = conn1D, & expected = reshape([1,2,2,3,3,4], [2,3]), & msg = "Linear Cn connectivity", & group = "connectivity") ! ---------------------------- ! Test: elemConn_C0 (2D) ! ---------------------------- conn2D = elemConn_C0(5, 5, 2, 2) call ut%test(40)%check( & name = "elemConn_C0_S", & res = shape(conn2D), & expected = [4,9], & msg = "Surface C0 connectivity shape", & group = "connectivity") ! ---------------------------- ! Test: elemConn_Cn (2D) ! ---------------------------- call elemConn_Cn(5, 5, 2, 2, [0.0_rk, 0.5_rk, 1.0_rk], [0.0_rk, 0.5_rk, 1.0_rk], & [3,1,3], [3,1,3], conn2D) call ut%test(41)%check( & name = "elemConn_Cn_S", & res = shape(conn2D), & expected = [4,9], & msg = "Surface Cn connectivity shape", & group = "connectivity") ! ---------------------------- ! Test: elemConn_C0 (3D) ! ---------------------------- conn3D = elemConn_C0(5, 5, 5, 2, 2, 2) call ut%test(42)%check( & name = "elemConn_C0_V", & res = shape(conn3D), & expected = [8,27], & msg = "Volume C0 connectivity shape", & group = "connectivity") ! ---------------------------- ! Test: elemConn_Cn (3D) ! ---------------------------- call elemConn_Cn(5, 5, 5, 2, 2, 2, & [0.0_rk, 0.5_rk, 1.0_rk], & [0.0_rk, 0.5_rk, 1.0_rk], & [0.0_rk, 0.5_rk, 1.0_rk], & [3,1,3], [3,1,3], [3,1,3], conn3D) call ut%test(43)%check( & name = "elemConn_Cn_V", & res = shape(conn3D), & expected = [8,27], & msg = "Volume Cn connectivity shape", & group = "connectivity") ! ---------------------------- ! Test: inv (3x3 matrix) ! ---------------------------- allocate(A2(3,3)) A2 = reshape([1.0_rk, 2.0_rk, 3.0_rk, & 0.0_rk, 1.0_rk, 4.0_rk, & 5.0_rk, 6.0_rk, 0.0_rk], [3,3]) A_inv = inv(A2) call ut%test(44)%check( & name = "inv_3x3", & res = matmul(A2, A_inv), & expected = eye(3), & msg = "A . inv(A) = I for 3x3", & group = "matrix") ! ---------------------------- ! Test: inv (rectangular 3x2 matrix) ! ---------------------------- deallocate(A2, A_inv) allocate(A2(3,2)) A2 = reshape([1.0_rk, 2.0_rk, 3.0_rk, 4.0_rk, 5.0_rk, 6.0_rk], [3,2]) A_inv = inv(A2) call ut%test(45)%check( & name = "inv_rectangular_3x2", & res = matmul(A_inv, A2), & expected = eye(2), & msg = "inv(A) . A = I for 3x2", & group = "matrix") ! ---------------------------- ! Test: inv (rectangular 2x3 matrix) ! ---------------------------- deallocate(A2, A_inv) allocate(A2(2,3)) A2 = reshape([1.0_rk, 4.0_rk, 2.0_rk, 5.0_rk, 3.0_rk, 6.0_rk], [2,3]) A_inv = inv(A2) call ut%test(46)%check( & name = "inv_rectangular_2x3", & res = matmul(A2, A_inv), & expected = eye(2), & msg = "A . inv(A) = I for 2x3", & group = "matrix") ! ---------------------------- ! Test: inv (identity matrix) ! ---------------------------- deallocate(A2, A_inv) allocate(A2(4,4)) A2 = eye(4) A_inv = inv(A2) call ut%test(47)%check( & name = "inv_identity", & res = A_inv, & expected = A2, & msg = "inv(I) = I", & group = "matrix") call ut%test(48)%check( & name = "inv_rectangular_2x3_proj", & res = matmul(A_inv, A2), & expected = transpose(matmul(A_inv, A2)), & msg = "inv(A) . A is symmetric projection (2x3)", & group = "matrix") ! ---------------------------- ! Test: kron_eye ! ---------------------------- if (allocated(A2)) deallocate(A2) allocate(A2(2,2)) A2 = reshape([1.0_rk, 2.0_rk, 3.0_rk, 4.0_rk],[2,2]) R = kron_eye(A2, 2) call ut%test(49)%check( & name = "kron_eye_block_diag", & res = R, & expected = reshape([ & 1.0_rk,0.0_rk,2.0_rk,0.0_rk, & 0.0_rk,1.0_rk,0.0_rk,2.0_rk, & 3.0_rk,0.0_rk,4.0_rk,0.0_rk, & 0.0_rk,3.0_rk,0.0_rk,4.0_rk], [4,4]), & msg = "Kronecker with identity blocks", & group = "kron") ! ---------------------------- ! Test: kron_t1_t1 (vector .kron. vector) ! ---------------------------- ! different sizes + negatives/zeros u2 = [-1.0_rk, 0.0_rk, 2.0_rk] v2 = [ 5.0_rk, -3.0_rk] w2 = kron(u2, v2) call ut%test(50)%check( & name = "kron_t1_t1_values", & res = w2, & expected = [-5.0_rk, 3.0_rk, 0.0_rk, 0.0_rk, 10.0_rk, -6.0_rk], & msg = "kron(u,v) concatenates u(i)*v blocks", & group = "kron") ! length matches size(u)*size(v) call ut%test(51)%check( & name = "kron_t1_t1_size", & res = size(w2), & expected = size(u2)*size(v2), & msg = "Length of kron(u,v)", & group = "kron") ! non-commutativity (order matters) call ut%test(52)%check( & name = "kron_t1_t1_noncommutative", & res = all(w2 == kron(v2, u2)), & expected = .false., & msg = "kron(u,v) /= kron(v,u) for vectors", & group = "kron") ! ---------------------------- ! Test: kron3(u, v, w) ! ---------------------------- ! values with zeros/negatives and |v|=2, |w|=1 u2 = [-1.0_rk, 2.0_rk] v2 = [0.0_rk, 3.0_rk] K3 = [-2.0_rk] out = kron(u2, v2, K3) ! kron3 call ut%test(53)%check( & name = "kron3_values_mixed", & res = out, & expected = [0.0_rk, 6.0_rk, 0.0_rk, -12.0_rk], & msg = "ordering: (u1*v1*w1, u1*v2*w1, u2*v1*w1, u2*v2*w1)", & group = "kron") ! size = |u|*|v|*|w| call ut%test(54)%check( & name = "kron3_size", & res = size(out), & expected = size(u2) * size(v2) * size(K3), & msg = "length of kron3 output", & group = "kron") ! associativity check: kron3(u,v,w) == kron(u, kron(v,w)) u2 = [1.0_rk, 2.0_rk] v2 = [3.0_rk, 4.0_rk] K3 = [5.0_rk] out = kron(u2, v2, K3) w2 = kron(v2, K3) A = kron(u2, w2) call ut%test(55)%check( & name = "kron3_associativity_vec", & res = all(abs(out - A) <= epsilon(0.0_rk)), & expected = .true., & msg = "kron3 equals kron(u, kron(v,w))", & group = "kron") ! another ordering/values case (|u|=2, |v|=1, |w|=3) u2 = [2.0_rk, -1.0_rk] v2 = [7.0_rk] K3 = [1.0_rk, 0.0_rk, -2.0_rk] out = kron(u2, v2, K3) call ut%test(56)%check( & name = "kron3_values_ordering", & res = out, & expected = [14.0_rk, 0.0_rk, -28.0_rk, -7.0_rk, 0.0_rk, 14.0_rk], & msg = "iterate i (u), then j (v), then k (w)", & group = "kron") ! non-commutativity across arguments u2 = [1.0_rk, 2.0_rk] v2 = [3.0_rk, 4.0_rk] K3 = [5.0_rk] ! length-1 w is fine; non-commutativity comes from u vs v call ut%test(57)%check( & name = "kron3_noncommutative", & res = all(kron(u2, v2, K3) == kron(v2, u2, K3)), & expected = .false., & msg = "kron3(u,v,w) /= kron3(v,u,w)", & group = "kron") ! summary of tests call ut%summary( & required_score = 100.0, & verbose = 3, & stop_fail = .false.) end program test_forcad_utils