test: add 3d euler beam red fortran tests
This commit is contained in:
김경종
@@ -0,0 +1,127 @@
|
||||
program test_abi_static
|
||||
use uel_3d_euler_beam_test_support
|
||||
use uel_3d_euler_beam_kernel, only: &
|
||||
UEL3DEB_OK, uel3deb_global_stiffness
|
||||
use uel_3d_euler_beam_abi_adapter, only: uel3deb_abi_static
|
||||
implicit none
|
||||
|
||||
call test_static_full_contribution()
|
||||
call test_stiffness_only_and_residual_only()
|
||||
write(*, '(A)') 'PASS uel_3d_euler_beam_abi_static'
|
||||
|
||||
contains
|
||||
|
||||
subroutine valid_problem(coords, props, u)
|
||||
real(dp), intent(out) :: coords(3, 2)
|
||||
real(dp), intent(out) :: props(9)
|
||||
real(dp), intent(out) :: u(12)
|
||||
|
||||
coords(:, 1) = [0.0_dp, 0.0_dp, 0.0_dp]
|
||||
coords(:, 2) = [2.0_dp, 0.0_dp, 0.0_dp]
|
||||
props = [210.0e9_dp, 80.0e9_dp, 3.0e-3_dp, 4.0e-6_dp, 7.0e-6_dp, &
|
||||
2.5e-6_dp, 0.0_dp, 1.0_dp, 0.0_dp]
|
||||
u = [1.0e-4_dp, -2.0e-4_dp, 3.0e-4_dp, 4.0e-5_dp, -5.0e-5_dp, &
|
||||
6.0e-5_dp, -7.0e-4_dp, 8.0e-4_dp, -9.0e-4_dp, 1.0e-4_dp, &
|
||||
-1.1e-4_dp, 1.2e-4_dp]
|
||||
end subroutine valid_problem
|
||||
|
||||
subroutine expected_outputs(coords, props, u, expected_k, expected_rhs)
|
||||
real(dp), intent(in) :: coords(3, 2)
|
||||
real(dp), intent(in) :: props(9)
|
||||
real(dp), intent(in) :: u(12)
|
||||
real(dp), intent(out) :: expected_k(12, 12)
|
||||
real(dp), intent(out) :: expected_rhs(12)
|
||||
real(dp) :: transform(12, 12)
|
||||
integer :: status
|
||||
|
||||
call uel3deb_global_stiffness(coords, props, expected_k, transform, status)
|
||||
call assert_equal_int(status, UEL3DEB_OK, 'expected K status')
|
||||
expected_rhs = -matmul(expected_k, u)
|
||||
end subroutine expected_outputs
|
||||
|
||||
subroutine call_adapter(lflags3, rhs, amatrx, energy, pnewdt, status)
|
||||
integer, intent(in) :: lflags3
|
||||
real(dp), intent(out) :: rhs(12, 1)
|
||||
real(dp), intent(out) :: amatrx(12, 12)
|
||||
real(dp), intent(out) :: energy(8)
|
||||
real(dp), intent(inout) :: pnewdt
|
||||
integer, intent(out) :: status
|
||||
real(dp) :: coords(3, 2)
|
||||
real(dp) :: props(9)
|
||||
real(dp) :: u(12)
|
||||
integer :: lflags(5)
|
||||
integer :: jprops(1)
|
||||
|
||||
call valid_problem(coords, props, u)
|
||||
rhs = 999.0_dp
|
||||
amatrx = 999.0_dp
|
||||
energy = 999.0_dp
|
||||
lflags = 0
|
||||
lflags(2) = 0
|
||||
lflags(3) = lflags3
|
||||
jprops = 0
|
||||
|
||||
call uel3deb_abi_static(rhs, amatrx, energy, 12, 1, 0, props, 9, coords, &
|
||||
3, 2, u, lflags, 12, 0, pnewdt, jprops, 0, status)
|
||||
end subroutine call_adapter
|
||||
|
||||
subroutine test_static_full_contribution()
|
||||
real(dp) :: coords(3, 2)
|
||||
real(dp) :: props(9)
|
||||
real(dp) :: u(12)
|
||||
real(dp) :: rhs(12, 1)
|
||||
real(dp) :: amatrx(12, 12)
|
||||
real(dp) :: energy(8)
|
||||
real(dp) :: expected_k(12, 12)
|
||||
real(dp) :: expected_rhs(12)
|
||||
real(dp) :: pnewdt
|
||||
integer :: status
|
||||
|
||||
call valid_problem(coords, props, u)
|
||||
call expected_outputs(coords, props, u, expected_k, expected_rhs)
|
||||
|
||||
pnewdt = 0.25_dp
|
||||
call call_adapter(1, rhs, amatrx, energy, pnewdt, status)
|
||||
|
||||
call assert_equal_int(status, UEL3DEB_OK, 'NOA-A-AMATRX-001 status LFLAGS3=1')
|
||||
call assert_matrix_close(amatrx, expected_k, tol_stiff, 'NOA-A-AMATRX-001 K')
|
||||
call assert_vector_close(rhs(:, 1), expected_rhs, tol_vector, 'NOA-A-RHS-001 RHS')
|
||||
call assert_vector_close(energy, spread(0.0_dp, 1, 8), tol_vector, 'NOA-A-ENERGY-001 ENERGY')
|
||||
call assert_close(pnewdt, 0.25_dp, tol_vector, 'NOA-A-ENERGY-001 PNEWDT')
|
||||
end subroutine test_static_full_contribution
|
||||
|
||||
subroutine test_stiffness_only_and_residual_only()
|
||||
real(dp) :: coords(3, 2)
|
||||
real(dp) :: props(9)
|
||||
real(dp) :: u(12)
|
||||
real(dp) :: rhs(12, 1)
|
||||
real(dp) :: amatrx(12, 12)
|
||||
real(dp) :: energy(8)
|
||||
real(dp) :: expected_k(12, 12)
|
||||
real(dp) :: expected_rhs(12)
|
||||
real(dp) :: zeros12(12)
|
||||
real(dp) :: zeros_k(12, 12)
|
||||
real(dp) :: pnewdt
|
||||
integer :: status
|
||||
|
||||
call valid_problem(coords, props, u)
|
||||
call expected_outputs(coords, props, u, expected_k, expected_rhs)
|
||||
zeros12 = 0.0_dp
|
||||
zeros_k = 0.0_dp
|
||||
|
||||
pnewdt = 0.5_dp
|
||||
call call_adapter(2, rhs, amatrx, energy, pnewdt, status)
|
||||
call assert_equal_int(status, UEL3DEB_OK, 'NOA-A-AMATRX-001 status LFLAGS3=2')
|
||||
call assert_matrix_close(amatrx, expected_k, tol_stiff, 'NOA-A-AMATRX-001 K only')
|
||||
call assert_vector_close(rhs(:, 1), zeros12, tol_vector, 'NOA-A-AMATRX-001 zero RHS')
|
||||
call assert_close(pnewdt, 0.5_dp, tol_vector, 'NOA-A-ENERGY-001 PNEWDT K only')
|
||||
|
||||
pnewdt = 0.75_dp
|
||||
call call_adapter(5, rhs, amatrx, energy, pnewdt, status)
|
||||
call assert_equal_int(status, UEL3DEB_OK, 'NOA-A-RHS-001 status LFLAGS3=5')
|
||||
call assert_matrix_close(amatrx, zeros_k, tol_stiff, 'NOA-A-RHS-001 zero AMATRX')
|
||||
call assert_vector_close(rhs(:, 1), expected_rhs, tol_vector, 'NOA-A-RHS-001 RHS only')
|
||||
call assert_close(pnewdt, 0.75_dp, tol_vector, 'NOA-A-ENERGY-001 PNEWDT RHS only')
|
||||
end subroutine test_stiffness_only_and_residual_only
|
||||
|
||||
end program test_abi_static
|
||||
@@ -0,0 +1,131 @@
|
||||
program test_invalid_inputs
|
||||
use ieee_arithmetic, only: ieee_quiet_nan, ieee_value
|
||||
use uel_3d_euler_beam_test_support
|
||||
use uel_3d_euler_beam_kernel, only: &
|
||||
UEL3DEB_E001_NDOFEL, UEL3DEB_E002_NNODE, UEL3DEB_E003_MCRD, &
|
||||
UEL3DEB_E004_NPROPS, UEL3DEB_E005_NJPROP, UEL3DEB_E006_NSVARS, &
|
||||
UEL3DEB_E007_NRHS, UEL3DEB_E008_MLVARX, UEL3DEB_E009_NONFINITE, &
|
||||
UEL3DEB_E010_NONPOSITIVE_PROPERTY, UEL3DEB_E011_ZERO_LENGTH, &
|
||||
UEL3DEB_E012_ZERO_ORIENTATION, UEL3DEB_E013_PARALLEL_ORIENTATION, &
|
||||
UEL3DEB_E014_LFLAGS2, UEL3DEB_E015_LFLAGS3, UEL3DEB_E016_NDLOAD, &
|
||||
uel3deb_validate_kernel_inputs
|
||||
use uel_3d_euler_beam_abi_adapter, only: uel3deb_abi_static
|
||||
implicit none
|
||||
|
||||
call test_kernel_physical_diagnostics()
|
||||
call test_adapter_shape_and_request_diagnostics()
|
||||
write(*, '(A)') 'PASS uel_3d_euler_beam_invalid_inputs'
|
||||
|
||||
contains
|
||||
|
||||
subroutine valid_kernel_inputs(coords, props)
|
||||
real(dp), intent(out) :: coords(3, 2)
|
||||
real(dp), intent(out) :: props(9)
|
||||
|
||||
coords(:, 1) = [0.0_dp, 0.0_dp, 0.0_dp]
|
||||
coords(:, 2) = [2.0_dp, 0.0_dp, 0.0_dp]
|
||||
props = [210.0e9_dp, 80.0e9_dp, 3.0e-3_dp, 4.0e-6_dp, 7.0e-6_dp, &
|
||||
2.5e-6_dp, 0.0_dp, 1.0_dp, 0.0_dp]
|
||||
end subroutine valid_kernel_inputs
|
||||
|
||||
subroutine test_kernel_physical_diagnostics()
|
||||
real(dp) :: coords(3, 2)
|
||||
real(dp) :: props(9)
|
||||
integer :: status
|
||||
|
||||
call valid_kernel_inputs(coords, props)
|
||||
coords(1, 1) = ieee_value(coords(1, 1), ieee_quiet_nan)
|
||||
call uel3deb_validate_kernel_inputs(coords, props, status)
|
||||
call assert_equal_int(status, UEL3DEB_E009_NONFINITE, 'NOA-I-PHYS-001 nonfinite coordinate')
|
||||
|
||||
call valid_kernel_inputs(coords, props)
|
||||
props(1) = -1.0_dp
|
||||
call uel3deb_validate_kernel_inputs(coords, props, status)
|
||||
call assert_equal_int(status, UEL3DEB_E010_NONPOSITIVE_PROPERTY, 'NOA-I-PHYS-001 nonpositive E')
|
||||
|
||||
call valid_kernel_inputs(coords, props)
|
||||
coords(:, 2) = coords(:, 1)
|
||||
call uel3deb_validate_kernel_inputs(coords, props, status)
|
||||
call assert_equal_int(status, UEL3DEB_E011_ZERO_LENGTH, 'NOA-I-PHYS-001 zero length')
|
||||
|
||||
call valid_kernel_inputs(coords, props)
|
||||
props(7:9) = [0.0_dp, 0.0_dp, 0.0_dp]
|
||||
call uel3deb_validate_kernel_inputs(coords, props, status)
|
||||
call assert_equal_int(status, UEL3DEB_E012_ZERO_ORIENTATION, 'NOA-I-PHYS-001 zero orientation')
|
||||
|
||||
call valid_kernel_inputs(coords, props)
|
||||
props(7:9) = [1.0_dp, 0.0_dp, 0.0_dp]
|
||||
call uel3deb_validate_kernel_inputs(coords, props, status)
|
||||
call assert_equal_int(status, UEL3DEB_E013_PARALLEL_ORIENTATION, 'NOA-I-PHYS-001 parallel orientation')
|
||||
end subroutine test_kernel_physical_diagnostics
|
||||
|
||||
subroutine test_adapter_shape_and_request_diagnostics()
|
||||
call expect_adapter_status(11, 1, 0, 9, 3, 2, 12, 0, 0, 0, 1, &
|
||||
UEL3DEB_E001_NDOFEL, 'NOA-I-SHAPE-001 NDOFEL')
|
||||
call expect_adapter_status(12, 1, 0, 9, 3, 1, 12, 0, 0, 0, 1, &
|
||||
UEL3DEB_E002_NNODE, 'NOA-I-SHAPE-001 NNODE')
|
||||
call expect_adapter_status(12, 1, 0, 9, 2, 2, 12, 0, 0, 0, 1, &
|
||||
UEL3DEB_E003_MCRD, 'NOA-I-SHAPE-001 MCRD')
|
||||
call expect_adapter_status(12, 1, 0, 8, 3, 2, 12, 0, 0, 0, 1, &
|
||||
UEL3DEB_E004_NPROPS, 'NOA-I-SHAPE-001 NPROPS')
|
||||
call expect_adapter_status(12, 1, 0, 9, 3, 2, 12, 0, 1, 0, 1, &
|
||||
UEL3DEB_E005_NJPROP, 'NOA-I-SHAPE-001 NJPROP')
|
||||
call expect_adapter_status(12, 1, 1, 9, 3, 2, 12, 0, 0, 0, 1, &
|
||||
UEL3DEB_E006_NSVARS, 'NOA-I-SHAPE-001 NSVARS')
|
||||
call expect_adapter_status(12, 2, 0, 9, 3, 2, 12, 0, 0, 0, 1, &
|
||||
UEL3DEB_E007_NRHS, 'NOA-I-SHAPE-001 NRHS')
|
||||
call expect_adapter_status(12, 1, 0, 9, 3, 2, 11, 0, 0, 0, 1, &
|
||||
UEL3DEB_E008_MLVARX, 'NOA-I-SHAPE-001 MLVARX')
|
||||
call expect_adapter_status(12, 1, 0, 9, 3, 2, 12, 0, 0, 1, 1, &
|
||||
UEL3DEB_E014_LFLAGS2, 'NOA-I-LFLAGS-001 LFLAGS2')
|
||||
call expect_adapter_status(12, 1, 0, 9, 3, 2, 12, 0, 0, 0, 3, &
|
||||
UEL3DEB_E015_LFLAGS3, 'NOA-I-LFLAGS-001 LFLAGS3')
|
||||
call expect_adapter_status(12, 1, 0, 9, 3, 2, 12, 1, 0, 0, 1, &
|
||||
UEL3DEB_E016_NDLOAD, 'NOA-I-LFLAGS-001 NDLOAD')
|
||||
end subroutine test_adapter_shape_and_request_diagnostics
|
||||
|
||||
subroutine expect_adapter_status(ndofel, nrhs, nsvars, nprops, mcrd, nnode, &
|
||||
mlvarx, ndload, njprop, lflag2, lflag3, &
|
||||
expected_status, message)
|
||||
integer, intent(in) :: ndofel
|
||||
integer, intent(in) :: nrhs
|
||||
integer, intent(in) :: nsvars
|
||||
integer, intent(in) :: nprops
|
||||
integer, intent(in) :: mcrd
|
||||
integer, intent(in) :: nnode
|
||||
integer, intent(in) :: mlvarx
|
||||
integer, intent(in) :: ndload
|
||||
integer, intent(in) :: njprop
|
||||
integer, intent(in) :: lflag2
|
||||
integer, intent(in) :: lflag3
|
||||
integer, intent(in) :: expected_status
|
||||
character(len=*), intent(in) :: message
|
||||
real(dp) :: rhs(12, 2)
|
||||
real(dp) :: amatrx(12, 12)
|
||||
real(dp) :: energy(8)
|
||||
real(dp) :: coords(3, 2)
|
||||
real(dp) :: props(9)
|
||||
real(dp) :: u(12)
|
||||
real(dp) :: pnewdt
|
||||
integer :: lflags(5)
|
||||
integer :: jprops(1)
|
||||
integer :: status
|
||||
|
||||
call valid_kernel_inputs(coords, props)
|
||||
rhs = 0.0_dp
|
||||
amatrx = 0.0_dp
|
||||
energy = 0.0_dp
|
||||
u = 0.0_dp
|
||||
pnewdt = 1.0_dp
|
||||
lflags = 0
|
||||
lflags(2) = lflag2
|
||||
lflags(3) = lflag3
|
||||
jprops = 0
|
||||
|
||||
call uel3deb_abi_static(rhs, amatrx, energy, ndofel, nrhs, nsvars, props, &
|
||||
nprops, coords, mcrd, nnode, u, lflags, mlvarx, &
|
||||
ndload, pnewdt, jprops, njprop, status)
|
||||
call assert_equal_int(status, expected_status, message)
|
||||
end subroutine expect_adapter_status
|
||||
|
||||
end program test_invalid_inputs
|
||||
@@ -0,0 +1,134 @@
|
||||
program test_kernel_stiffness
|
||||
use uel_3d_euler_beam_test_support
|
||||
use uel_3d_euler_beam_kernel, only: &
|
||||
UEL3DEB_OK, uel3deb_local_stiffness
|
||||
implicit none
|
||||
|
||||
call test_local_stiffness_entries()
|
||||
call test_analytical_responses()
|
||||
write(*, '(A)') 'PASS uel_3d_euler_beam_kernel_stiffness'
|
||||
|
||||
contains
|
||||
|
||||
subroutine beam_values(length, young, shear, area, iy, iz, torsion_j)
|
||||
real(dp), intent(out) :: length
|
||||
real(dp), intent(out) :: young
|
||||
real(dp), intent(out) :: shear
|
||||
real(dp), intent(out) :: area
|
||||
real(dp), intent(out) :: iy
|
||||
real(dp), intent(out) :: iz
|
||||
real(dp), intent(out) :: torsion_j
|
||||
|
||||
length = 2.0_dp
|
||||
young = 210.0e9_dp
|
||||
shear = 80.0e9_dp
|
||||
area = 3.0e-3_dp
|
||||
iy = 4.0e-6_dp
|
||||
iz = 7.0e-6_dp
|
||||
torsion_j = 2.5e-6_dp
|
||||
end subroutine beam_values
|
||||
|
||||
subroutine build_local(k, length, young, shear, area, iy, iz, torsion_j)
|
||||
real(dp), intent(out) :: k(12, 12)
|
||||
real(dp), intent(out) :: length
|
||||
real(dp), intent(out) :: young
|
||||
real(dp), intent(out) :: shear
|
||||
real(dp), intent(out) :: area
|
||||
real(dp), intent(out) :: iy
|
||||
real(dp), intent(out) :: iz
|
||||
real(dp), intent(out) :: torsion_j
|
||||
integer :: status
|
||||
|
||||
call beam_values(length, young, shear, area, iy, iz, torsion_j)
|
||||
call uel3deb_local_stiffness(length, young, shear, area, iy, iz, torsion_j, k, status)
|
||||
call assert_equal_int(status, UEL3DEB_OK, 'local stiffness status')
|
||||
end subroutine build_local
|
||||
|
||||
subroutine test_local_stiffness_entries()
|
||||
real(dp) :: k(12, 12)
|
||||
real(dp) :: length, young, shear, area, iy, iz, torsion_j
|
||||
real(dp) :: axial, torsion, b2_12, b2_6, b2_4, b2_2
|
||||
real(dp) :: b3_12, b3_6, b3_4, b3_2
|
||||
|
||||
call build_local(k, length, young, shear, area, iy, iz, torsion_j)
|
||||
|
||||
axial = young * area / length
|
||||
torsion = shear * torsion_j / length
|
||||
b2_12 = 12.0_dp * young * iy / length**3
|
||||
b2_6 = 6.0_dp * young * iy / length**2
|
||||
b2_4 = 4.0_dp * young * iy / length
|
||||
b2_2 = 2.0_dp * young * iy / length
|
||||
b3_12 = 12.0_dp * young * iz / length**3
|
||||
b3_6 = 6.0_dp * young * iz / length**2
|
||||
b3_4 = 4.0_dp * young * iz / length
|
||||
b3_2 = 2.0_dp * young * iz / length
|
||||
|
||||
call assert_matrix_symmetric(k, tol_symmetry, 'NOA-K-SYM-001 local symmetry')
|
||||
|
||||
call assert_close(k(1, 1), axial, tol_stiff, 'NOA-K-STIFF-001 axial k11')
|
||||
call assert_close(k(1, 7), -axial, tol_stiff, 'NOA-K-STIFF-001 axial k17')
|
||||
call assert_close(k(7, 7), axial, tol_stiff, 'NOA-K-STIFF-001 axial k77')
|
||||
call assert_close(k(4, 4), torsion, tol_stiff, 'NOA-K-STIFF-001 torsion k44')
|
||||
call assert_close(k(4, 10), -torsion, tol_stiff, 'NOA-K-STIFF-001 torsion k4,10')
|
||||
call assert_close(k(10, 10), torsion, tol_stiff, 'NOA-K-STIFF-001 torsion k10,10')
|
||||
|
||||
call assert_close(k(3, 3), b2_12, tol_stiff, 'NOA-K-BEND2-001 k33')
|
||||
call assert_close(k(3, 5), -b2_6, tol_stiff, 'NOA-K-BEND2-001 k35')
|
||||
call assert_close(k(3, 9), -b2_12, tol_stiff, 'NOA-K-BEND2-001 k39')
|
||||
call assert_close(k(3, 11), -b2_6, tol_stiff, 'NOA-K-BEND2-001 k3,11')
|
||||
call assert_close(k(5, 5), b2_4, tol_stiff, 'NOA-K-BEND2-001 k55')
|
||||
call assert_close(k(5, 9), b2_6, tol_stiff, 'NOA-K-BEND2-001 k59')
|
||||
call assert_close(k(5, 11), b2_2, tol_stiff, 'NOA-K-BEND2-001 k5,11')
|
||||
call assert_close(k(11, 11), b2_4, tol_stiff, 'NOA-K-BEND2-001 k11,11')
|
||||
|
||||
call assert_close(k(2, 2), b3_12, tol_stiff, 'NOA-K-BEND3-001 k22')
|
||||
call assert_close(k(2, 6), b3_6, tol_stiff, 'NOA-K-BEND3-001 k26')
|
||||
call assert_close(k(2, 8), -b3_12, tol_stiff, 'NOA-K-BEND3-001 k28')
|
||||
call assert_close(k(2, 12), b3_6, tol_stiff, 'NOA-K-BEND3-001 k2,12')
|
||||
call assert_close(k(6, 6), b3_4, tol_stiff, 'NOA-K-BEND3-001 k66')
|
||||
call assert_close(k(6, 8), -b3_6, tol_stiff, 'NOA-K-BEND3-001 k68')
|
||||
call assert_close(k(6, 12), b3_2, tol_stiff, 'NOA-K-BEND3-001 k6,12')
|
||||
call assert_close(k(12, 12), b3_4, tol_stiff, 'NOA-K-BEND3-001 k12,12')
|
||||
call assert_close(k(2, 3), 0.0_dp, tol_stiff, 'NOA-K-STIFF-001 uncoupled k23')
|
||||
end subroutine test_local_stiffness_entries
|
||||
|
||||
subroutine test_analytical_responses()
|
||||
real(dp) :: k(12, 12)
|
||||
real(dp) :: force(12)
|
||||
real(dp) :: u(12)
|
||||
real(dp) :: length, young, shear, area, iy, iz, torsion_j
|
||||
real(dp) :: axial, torsion, b2_12, b3_12
|
||||
|
||||
call build_local(k, length, young, shear, area, iy, iz, torsion_j)
|
||||
|
||||
axial = young * area / length
|
||||
torsion = shear * torsion_j / length
|
||||
b2_12 = 12.0_dp * young * iy / length**3
|
||||
b3_12 = 12.0_dp * young * iz / length**3
|
||||
|
||||
u = 0.0_dp
|
||||
u(7) = 1.0e-2_dp
|
||||
force = matmul(k, u)
|
||||
call assert_close(force(1), -axial * u(7), tol_vector, 'NOA-K-STIFF-002 f1')
|
||||
call assert_close(force(7), axial * u(7), tol_vector, 'NOA-K-STIFF-002 f7')
|
||||
|
||||
u = 0.0_dp
|
||||
u(10) = 2.0e-2_dp
|
||||
force = matmul(k, u)
|
||||
call assert_close(force(4), -torsion * u(10), tol_vector, 'NOA-K-STIFF-003 m1')
|
||||
call assert_close(force(10), torsion * u(10), tol_vector, 'NOA-K-STIFF-003 m2')
|
||||
|
||||
u = 0.0_dp
|
||||
u(9) = 1.0e-2_dp
|
||||
force = matmul(k, u)
|
||||
call assert_close(force(3), -b2_12 * u(9), tol_vector, 'NOA-K-BEND2-001 f3')
|
||||
call assert_close(force(9), b2_12 * u(9), tol_vector, 'NOA-K-BEND2-001 f9')
|
||||
|
||||
u = 0.0_dp
|
||||
u(8) = 1.0e-2_dp
|
||||
force = matmul(k, u)
|
||||
call assert_close(force(2), -b3_12 * u(8), tol_vector, 'NOA-K-BEND3-001 f2')
|
||||
call assert_close(force(8), b3_12 * u(8), tol_vector, 'NOA-K-BEND3-001 f8')
|
||||
end subroutine test_analytical_responses
|
||||
|
||||
end program test_kernel_stiffness
|
||||
@@ -0,0 +1,168 @@
|
||||
program test_kernel_transform_modes
|
||||
use uel_3d_euler_beam_test_support
|
||||
use uel_3d_euler_beam_kernel, only: &
|
||||
UEL3DEB_OK, uel3deb_global_stiffness, uel3deb_local_stiffness
|
||||
implicit none
|
||||
|
||||
call test_identity_and_rotated_transform()
|
||||
call test_rigid_body_modes()
|
||||
call test_reversed_node_order()
|
||||
write(*, '(A)') 'PASS uel_3d_euler_beam_kernel_transform_modes'
|
||||
|
||||
contains
|
||||
|
||||
subroutine valid_props(props)
|
||||
real(dp), intent(out) :: props(9)
|
||||
|
||||
props = [210.0e9_dp, 80.0e9_dp, 3.0e-3_dp, 4.0e-6_dp, 7.0e-6_dp, &
|
||||
2.5e-6_dp, 0.0_dp, 1.0_dp, 0.0_dp]
|
||||
end subroutine valid_props
|
||||
|
||||
subroutine fill_block_transform(rotation, transform)
|
||||
real(dp), intent(in) :: rotation(3, 3)
|
||||
real(dp), intent(out) :: transform(12, 12)
|
||||
integer :: block
|
||||
integer :: i
|
||||
integer :: j
|
||||
|
||||
transform = 0.0_dp
|
||||
do block = 0, 3
|
||||
do i = 1, 3
|
||||
do j = 1, 3
|
||||
transform(3 * block + i, 3 * block + j) = rotation(i, j)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end subroutine fill_block_transform
|
||||
|
||||
subroutine test_identity_and_rotated_transform()
|
||||
real(dp) :: props(9)
|
||||
real(dp) :: coords(3, 2)
|
||||
real(dp) :: k_local(12, 12)
|
||||
real(dp) :: k_global(12, 12)
|
||||
real(dp) :: transform(12, 12)
|
||||
real(dp) :: expected(12, 12)
|
||||
real(dp) :: expected_transform(12, 12)
|
||||
real(dp) :: rotation(3, 3)
|
||||
integer :: status
|
||||
|
||||
call valid_props(props)
|
||||
coords(:, 1) = [0.0_dp, 0.0_dp, 0.0_dp]
|
||||
coords(:, 2) = [2.0_dp, 0.0_dp, 0.0_dp]
|
||||
|
||||
call uel3deb_local_stiffness(2.0_dp, props(1), props(2), props(3), props(4), &
|
||||
props(5), props(6), k_local, status)
|
||||
call assert_equal_int(status, UEL3DEB_OK, 'identity local stiffness status')
|
||||
call uel3deb_global_stiffness(coords, props, k_global, transform, status)
|
||||
call assert_equal_int(status, UEL3DEB_OK, 'identity global stiffness status')
|
||||
call set_identity(expected_transform)
|
||||
call assert_matrix_close(transform, expected_transform, tol_stiff, 'NOA-K-ROT-001 identity T')
|
||||
call assert_matrix_close(k_global, k_local, tol_stiff, 'NOA-K-ROT-001 identity K')
|
||||
|
||||
props(7:9) = [0.0_dp, 0.0_dp, 1.0_dp]
|
||||
coords(:, 2) = [0.0_dp, 2.0_dp, 0.0_dp]
|
||||
call uel3deb_global_stiffness(coords, props, k_global, transform, status)
|
||||
call assert_equal_int(status, UEL3DEB_OK, 'rotated global stiffness status')
|
||||
|
||||
rotation = 0.0_dp
|
||||
rotation(1, 2) = 1.0_dp
|
||||
rotation(2, 3) = 1.0_dp
|
||||
rotation(3, 1) = 1.0_dp
|
||||
call fill_block_transform(rotation, expected_transform)
|
||||
expected = matmul(transpose(expected_transform), matmul(k_local, expected_transform))
|
||||
call assert_matrix_close(transform, expected_transform, tol_stiff, 'NOA-K-ROT-001 rotated T')
|
||||
call assert_matrix_close(k_global, expected, tol_stiff, 'NOA-K-ROT-001 rotated K')
|
||||
call assert_matrix_symmetric(k_global, tol_symmetry, 'NOA-K-SYM-001 global symmetry')
|
||||
end subroutine test_identity_and_rotated_transform
|
||||
|
||||
subroutine test_rigid_body_modes()
|
||||
real(dp) :: props(9)
|
||||
real(dp) :: coords(3, 2)
|
||||
real(dp) :: k_global(12, 12)
|
||||
real(dp) :: transform(12, 12)
|
||||
real(dp) :: mode(12)
|
||||
integer :: status
|
||||
|
||||
call valid_props(props)
|
||||
coords(:, 1) = [0.0_dp, 0.0_dp, 0.0_dp]
|
||||
coords(:, 2) = [2.0_dp, 0.0_dp, 0.0_dp]
|
||||
call uel3deb_global_stiffness(coords, props, k_global, transform, status)
|
||||
call assert_equal_int(status, UEL3DEB_OK, 'rigid mode stiffness status')
|
||||
|
||||
mode = 0.0_dp
|
||||
mode([1, 7]) = 1.0_dp
|
||||
call assert_rigid_mode(k_global, mode, 'NOA-K-RBM-001 rigid U1')
|
||||
|
||||
mode = 0.0_dp
|
||||
mode([2, 8]) = 1.0_dp
|
||||
call assert_rigid_mode(k_global, mode, 'NOA-K-RBM-001 rigid U2')
|
||||
|
||||
mode = 0.0_dp
|
||||
mode([3, 9]) = 1.0_dp
|
||||
call assert_rigid_mode(k_global, mode, 'NOA-K-RBM-001 rigid U3')
|
||||
|
||||
mode = 0.0_dp
|
||||
mode([4, 10]) = 1.0_dp
|
||||
call assert_rigid_mode(k_global, mode, 'NOA-K-RBM-001 rigid UR1')
|
||||
|
||||
mode = 0.0_dp
|
||||
mode(3) = 0.0_dp
|
||||
mode(5) = 1.0_dp
|
||||
mode(9) = -2.0_dp
|
||||
mode(11) = 1.0_dp
|
||||
call assert_rigid_mode(k_global, mode, 'NOA-K-RBM-001 rigid UR2')
|
||||
|
||||
mode = 0.0_dp
|
||||
mode(2) = 0.0_dp
|
||||
mode(6) = 1.0_dp
|
||||
mode(8) = 2.0_dp
|
||||
mode(12) = 1.0_dp
|
||||
call assert_rigid_mode(k_global, mode, 'NOA-K-RBM-001 rigid UR3')
|
||||
end subroutine test_rigid_body_modes
|
||||
|
||||
subroutine assert_rigid_mode(k_global, mode, message)
|
||||
real(dp), intent(in) :: k_global(12, 12)
|
||||
real(dp), intent(in) :: mode(12)
|
||||
character(len=*), intent(in) :: message
|
||||
real(dp) :: residual(12)
|
||||
real(dp) :: limit
|
||||
|
||||
residual = matmul(k_global, mode)
|
||||
limit = tol_vector * max(1.0_dp, max_abs_matrix(k_global) * max_abs_vector(mode))
|
||||
call assert_vector_norm_le(residual, limit, message)
|
||||
end subroutine assert_rigid_mode
|
||||
|
||||
subroutine test_reversed_node_order()
|
||||
real(dp) :: props(9)
|
||||
real(dp) :: coords_forward(3, 2)
|
||||
real(dp) :: coords_reversed(3, 2)
|
||||
real(dp) :: k_forward(12, 12)
|
||||
real(dp) :: k_reversed(12, 12)
|
||||
real(dp) :: transform(12, 12)
|
||||
real(dp) :: expected_reversed(12, 12)
|
||||
integer :: status
|
||||
integer :: i
|
||||
integer :: j
|
||||
integer :: permutation(12)
|
||||
|
||||
call valid_props(props)
|
||||
coords_forward(:, 1) = [0.0_dp, 0.0_dp, 0.0_dp]
|
||||
coords_forward(:, 2) = [2.0_dp, 0.0_dp, 0.0_dp]
|
||||
coords_reversed(:, 1) = coords_forward(:, 2)
|
||||
coords_reversed(:, 2) = coords_forward(:, 1)
|
||||
|
||||
call uel3deb_global_stiffness(coords_forward, props, k_forward, transform, status)
|
||||
call assert_equal_int(status, UEL3DEB_OK, 'forward node stiffness status')
|
||||
call uel3deb_global_stiffness(coords_reversed, props, k_reversed, transform, status)
|
||||
call assert_equal_int(status, UEL3DEB_OK, 'reversed node stiffness status')
|
||||
|
||||
permutation = [7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6]
|
||||
do i = 1, 12
|
||||
do j = 1, 12
|
||||
expected_reversed(i, j) = k_forward(permutation(i), permutation(j))
|
||||
end do
|
||||
end do
|
||||
call assert_matrix_close(k_reversed, expected_reversed, tol_stiff, 'NOA-K-REVNODE-001')
|
||||
end subroutine test_reversed_node_order
|
||||
|
||||
end program test_kernel_transform_modes
|
||||
@@ -0,0 +1,138 @@
|
||||
module uel_3d_euler_beam_test_support
|
||||
use iso_fortran_env, only: real64
|
||||
implicit none
|
||||
|
||||
integer, parameter :: dp = real64
|
||||
real(dp), parameter :: tol_stiff = 1.0e-10_dp
|
||||
real(dp), parameter :: tol_vector = 1.0e-10_dp
|
||||
real(dp), parameter :: tol_symmetry = 1.0e-12_dp
|
||||
|
||||
contains
|
||||
|
||||
subroutine fail_test(message)
|
||||
character(len=*), intent(in) :: message
|
||||
|
||||
write(*, '(A)') 'FAIL: ' // trim(message)
|
||||
error stop 1
|
||||
end subroutine fail_test
|
||||
|
||||
subroutine assert_true(condition, message)
|
||||
logical, intent(in) :: condition
|
||||
character(len=*), intent(in) :: message
|
||||
|
||||
if (.not. condition) call fail_test(message)
|
||||
end subroutine assert_true
|
||||
|
||||
subroutine assert_equal_int(actual, expected, message)
|
||||
integer, intent(in) :: actual
|
||||
integer, intent(in) :: expected
|
||||
character(len=*), intent(in) :: message
|
||||
|
||||
if (actual /= expected) then
|
||||
write(*, *) 'actual=', actual, ' expected=', expected
|
||||
call fail_test(message)
|
||||
end if
|
||||
end subroutine assert_equal_int
|
||||
|
||||
subroutine assert_close(actual, expected, tolerance, message)
|
||||
real(dp), intent(in) :: actual
|
||||
real(dp), intent(in) :: expected
|
||||
real(dp), intent(in) :: tolerance
|
||||
character(len=*), intent(in) :: message
|
||||
real(dp) :: limit
|
||||
|
||||
limit = tolerance * max(1.0_dp, abs(expected))
|
||||
if (abs(actual - expected) > limit) then
|
||||
write(*, *) 'actual=', actual, ' expected=', expected, ' limit=', limit
|
||||
call fail_test(message)
|
||||
end if
|
||||
end subroutine assert_close
|
||||
|
||||
function max_abs_matrix(values) result(max_abs)
|
||||
real(dp), intent(in) :: values(:, :)
|
||||
real(dp) :: max_abs
|
||||
|
||||
max_abs = maxval(abs(values))
|
||||
end function max_abs_matrix
|
||||
|
||||
function max_abs_vector(values) result(max_abs)
|
||||
real(dp), intent(in) :: values(:)
|
||||
real(dp) :: max_abs
|
||||
|
||||
max_abs = maxval(abs(values))
|
||||
end function max_abs_vector
|
||||
|
||||
subroutine assert_matrix_close(actual, expected, tolerance, message)
|
||||
real(dp), intent(in) :: actual(:, :)
|
||||
real(dp), intent(in) :: expected(:, :)
|
||||
real(dp), intent(in) :: tolerance
|
||||
character(len=*), intent(in) :: message
|
||||
real(dp) :: diff
|
||||
real(dp) :: limit
|
||||
|
||||
call assert_true(all(shape(actual) == shape(expected)), trim(message) // ' shape')
|
||||
diff = maxval(abs(actual - expected))
|
||||
limit = tolerance * max(1.0_dp, max_abs_matrix(expected))
|
||||
if (diff > limit) then
|
||||
write(*, *) 'max matrix diff=', diff, ' limit=', limit
|
||||
call fail_test(message)
|
||||
end if
|
||||
end subroutine assert_matrix_close
|
||||
|
||||
subroutine assert_vector_close(actual, expected, tolerance, message)
|
||||
real(dp), intent(in) :: actual(:)
|
||||
real(dp), intent(in) :: expected(:)
|
||||
real(dp), intent(in) :: tolerance
|
||||
character(len=*), intent(in) :: message
|
||||
real(dp) :: diff
|
||||
real(dp) :: limit
|
||||
|
||||
call assert_true(size(actual) == size(expected), trim(message) // ' size')
|
||||
diff = maxval(abs(actual - expected))
|
||||
limit = tolerance * max(1.0_dp, max_abs_vector(expected))
|
||||
if (diff > limit) then
|
||||
write(*, *) 'max vector diff=', diff, ' limit=', limit
|
||||
call fail_test(message)
|
||||
end if
|
||||
end subroutine assert_vector_close
|
||||
|
||||
subroutine assert_vector_norm_le(actual, limit, message)
|
||||
real(dp), intent(in) :: actual(:)
|
||||
real(dp), intent(in) :: limit
|
||||
character(len=*), intent(in) :: message
|
||||
real(dp) :: norm_inf
|
||||
|
||||
norm_inf = max_abs_vector(actual)
|
||||
if (norm_inf > limit) then
|
||||
write(*, *) 'max vector norm=', norm_inf, ' limit=', limit
|
||||
call fail_test(message)
|
||||
end if
|
||||
end subroutine assert_vector_norm_le
|
||||
|
||||
subroutine assert_matrix_symmetric(values, tolerance, message)
|
||||
real(dp), intent(in) :: values(:, :)
|
||||
real(dp), intent(in) :: tolerance
|
||||
character(len=*), intent(in) :: message
|
||||
real(dp) :: diff
|
||||
real(dp) :: limit
|
||||
|
||||
call assert_true(size(values, 1) == size(values, 2), trim(message) // ' square')
|
||||
diff = maxval(abs(values - transpose(values)))
|
||||
limit = tolerance * max(1.0_dp, max_abs_matrix(values))
|
||||
if (diff > limit) then
|
||||
write(*, *) 'symmetry diff=', diff, ' limit=', limit
|
||||
call fail_test(message)
|
||||
end if
|
||||
end subroutine assert_matrix_symmetric
|
||||
|
||||
subroutine set_identity(values)
|
||||
real(dp), intent(out) :: values(:, :)
|
||||
integer :: i
|
||||
|
||||
values = 0.0_dp
|
||||
do i = 1, min(size(values, 1), size(values, 2))
|
||||
values(i, i) = 1.0_dp
|
||||
end do
|
||||
end subroutine set_identity
|
||||
|
||||
end module uel_3d_euler_beam_test_support
|
||||
Reference in New Issue
Block a user