module test_polygon_in_mesh_m !! Unit tests for polygon_on_mesh_t use funit use helpers_m, only : almost_equal implicit none contains @test subroutine test_polygon_in_mesh_cerfons() !! Tests polygon on mesh in CERFONS geometry use precision_grillix_m, only : GP use descriptors_m, only : DISTRICT_CORE, DISTRICT_WALL, DISTRICT_DOME, DISTRICT_OUT use equilibrium_m, only : equilibrium_t use mesh_cart_m, only : mesh_cart_t use helpers_mesh_sample_m, only : mesh_sample_cerfons use polygon_in_mesh_m, only : polygon_in_mesh_t class(equilibrium_t), allocatable :: equi type(mesh_cart_t) :: mesh integer, parameter :: npoly = 5, gnseg = 6 logical :: lret, fine integer :: it, ki, l, district, ipoly, iret, iseg, gseg, ls, ln, k, nshft integer :: ncore, nwall, ndome, nout integer, allocatable, dimension(:) :: inds_core, inds_wall, inds_dome, inds_out, inds_comb, inds_corner type(polygon_in_mesh_t), dimension(npoly) :: polys type(polygon_in_mesh_t) :: poly_empty, poly_single real(GP) :: rnd, xs, ys, xn, yn, lseg, aseg ! Reference data integer, dimension(npoly), parameter :: ref_nsegs = [1, 1, 1, 2, 1] logical, dimension(npoly), parameter :: ref_is_periodic = [.true., .false., .false., .false., .true.] integer, dimension(npoly), parameter :: ref_np_total = [129, 375, 30, 9, 414] integer, dimension(gnseg), parameter :: ref_nps = [129, 375, 30, 5, 4, 414] real(GP), dimension(gnseg), parameter :: ref_lseg = [1.509533188058E+000_GP, 4.299188309204E+000_GP, 3.769848480983E-001_GP, & 4.000000000000E-002_GP, 3.000000000000E-002_GP, 4.786173157302E+000_GP] real(GP), dimension(gnseg), parameter :: ref_aseg = [1.513500000000E-001_GP, 8.813000000000E-001_GP, 1.108500000000E-001_GP, & 2.000000000000E-002_GP, 1.500000000000E-002_GP, -1.038900000000E+000_GP] write(*,*)'' write(*,'(A80)')'test_polygon_on_mesh_cerfons ----------------------------------------------' call mesh_sample_cerfons(equi, 0.0_GP, mesh) ! Create some indices for polygons do it = 1, 2 ncore = 0 nwall = 0 ndome = 0 nout = 0 do ki = 1, mesh%get_n_points_boundary() l = mesh%boundary_indices(ki) district = mesh%get_district(l) if (district == DISTRICT_CORE) then ncore = ncore + 1 if (it == 2) then inds_core(ncore) = l endif elseif (district == DISTRICT_WALL) then nwall = nwall + 1 if (it == 2) then inds_wall(nwall) = l inds_comb(nwall+ndome+nout) = l endif elseif (district == DISTRICT_DOME) then ndome = ndome + 1 if (it == 2) then inds_dome(ndome) = l inds_comb(nwall+ndome+nout) = l endif elseif (district == DISTRICT_OUT) then nout = nout + 1 if (it == 2) then inds_out(nout) = l inds_comb(nwall+ndome+nout) = l endif endif enddo if (it == 1) then allocate(inds_core(ncore)) allocate(inds_wall(nwall)) allocate(inds_dome(ndome)) allocate(inds_out(nout)) allocate(inds_comb(nwall+ndome+nout)) endif enddo ! Create polygons call polys(1)%init(equi, mesh, inds_core, dbgout=1) call polys(2)%init(equi, mesh, inds_wall, dbgout=1) call polys(3)%init(equi, mesh, inds_dome, dbgout=1) call polys(4)%init(equi, mesh, inds_out, dbgout=1) call polys(5)%init(equi, mesh, inds_comb, dir_counterclockwise=.false., dbgout=1) ! Test polygons gseg = 0 do ipoly = 1, npoly lret = polys(ipoly)%is_periodic() @assertEqual(lret, ref_is_periodic(ipoly)) iret = polys(ipoly)%get_np_total() @assertEqual(iret, ref_np_total(ipoly)) if (polys(ipoly)%is_periodic()) then ! Create a random shift index between [-1000, 1000] for periodic polygons ! to check in subsequent tests if periodicity is accounted for correctly call random_number(rnd) nshft = -1000 + FLOOR((1000+1+1000)*rnd) else nshft = 0 endif iret = polys(ipoly)%get_nsegs() @assertEqual(iret, ref_nsegs(ipoly)) ! Check individual segments do iseg = 1, polys(ipoly)%get_nsegs() gseg = gseg + 1 iret = polys(ipoly)%get_nps(iseg) @assertEqual(iret, ref_nps(gseg)) ! Check length and area (actually only for periodic) of segment lseg = 0.0_GP aseg = 0.0_GP do k = 1, polys(ipoly)%get_nps(iseg) ls = polys(ipoly)%get_ind_on_mesh(k+nshft, iseg) ln = polys(ipoly)%get_ind_on_mesh(k+nshft+1, iseg) if (ln /= 0) then xs = mesh%get_x(ls) ys = mesh%get_y(ls) xn = mesh%get_x(ln) yn = mesh%get_y(ln) lseg = lseg + sqrt((xn - xs)**2 + (yn - ys)**2) aseg = aseg + 0.5_GP * (xs * yn - xn * ys) endif enddo fine = almost_equal(lseg, ref_lseg(gseg), rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) fine = almost_equal(aseg, ref_aseg(gseg), rtol=0.0_GP, atol=1.0E-8_GP) @assertTrue(fine) enddo enddo ! Check corner case: empty list allocate(inds_corner(0)) call poly_empty%init(equi, mesh, inds_corner, dbgout=1) lret = poly_empty%is_periodic() @assertEqual(lret, .false.) iret = poly_empty%get_np_total() @assertEqual(iret, 0) iret = poly_empty%get_nsegs() @assertEqual(iret, 0) iret = poly_empty%get_ind_on_mesh(1) @assertEqual(iret, 0) deallocate(inds_corner) ! Check corner case: list with single entry allocate(inds_corner(1)) inds_corner(1) = inds_core(1) call poly_single%init(equi, mesh, inds_corner, dbgout=1) lret = poly_single%is_periodic() @assertEqual(lret, .false.) iret = poly_single%get_np_total() @assertEqual(iret, 0) iret = poly_single%get_nsegs() @assertEqual(iret, 0) iret = poly_single%get_ind_on_mesh(1) @assertEqual(iret, 0) deallocate(inds_corner) deallocate(inds_comb) deallocate(inds_out) deallocate(inds_dome) deallocate(inds_wall) deallocate(inds_core) write(*,'(A80)')'---------------------------------------------------------------------------' write(*,*)'' write(*,*)'' end subroutine end module