From 67de22512d3d9033467b13a86b92e570a70b2d57 Mon Sep 17 00:00:00 2001 From: FlorisBuwaldaDeltares Date: Mon, 1 Dec 2025 13:46:33 +0100 Subject: [PATCH 01/32] initial implementation of new, elemental routine --- .../prepost/pol_to_cellmask.f90 | 256 ++++++++++++++++-- 1 file changed, 227 insertions(+), 29 deletions(-) diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.f90 index a909b61287..3fbf644f4c 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.f90 +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.f90 @@ -43,51 +43,249 @@ module m_pol_to_cellmask public :: pol_to_cellmask + ! Module-level variables for cellmask polygon checking + real(kind=dp), allocatable :: xpmin_cellmask(:), ypmin_cellmask(:) + real(kind=dp), allocatable :: xpmax_cellmask(:), ypmax_cellmask(:) + real(kind=dp), allocatable :: zpl_start_cellmask(:) + real(kind=dp), allocatable :: xpl_cellmask(:), ypl_cellmask(:) + integer, allocatable :: iistart_cellmask(:), iiend_cellmask(:) + integer :: Npoly_cellmask = 0 + logical :: cellmask_initialized = .false. + contains subroutine pol_to_cellmask() - use network_data, only: cellmask, nump1d2d, npl, nump, xzw, yzw, xpl, ypl, zpl - use m_readyy, only: readyy - use m_missing, only: dmiss, JINS - use geometry_module, only: dbpinpol_optinside_perpol + use network_data, only: cellmask, nump1d2d, npl, nump, xzw, yzw, xpl, ypl, zpl + + if (allocated(cellmask)) deallocate(cellmask) + allocate(cellmask(nump1d2d)) + cellmask = 0 + + if (NPL == 0) return + + ! Initialize once + call dbpinpol_cellmask_init(NPL, xpl, ypl, zpl) + + ! Vectorized assignment - beautiful and simple! + cellmask(1:nump) = dbpinpol_cellmask(xzw(1:nump), yzw(1:nump)) + + ! Cleanup + call dbpinpol_cellmask_cleanup() - integer :: in, k, KMOD - integer :: num +end subroutine pol_to_cellmask - if (allocated(cellmask)) then - deallocate (cellmask) + !> Initialize polygon data structures for cellmask checking. + !! Must be called once before using dbpinpol_cellmask. + subroutine dbpinpol_cellmask_init(NPL, xpl, ypl, zpl) + use m_alloc + use m_missing, only: dp, dmiss + + implicit none + + integer, intent(in) :: NPL + real(kind=dp), intent(in) :: xpl(NPL), ypl(NPL), zpl(NPL) + + integer :: MAXPOLY + integer :: ipoint, istart, iend, ipoly + + if (cellmask_initialized) then + call dbpinpol_cellmask_cleanup() end if - allocate (cellmask(nump1d2d)) - cellmask = 0 - if (NPL == 0) return ! no polygon + if (NPL == 0) then + cellmask_initialized = .true. + return + end if - call READYY('Applying polygon cellmask', 0.0_dp) - KMOD = max(1, NUMP / 100) + MAXPOLY = 1000 + call realloc(xpmin_cellmask, maxpoly, keepExisting=.false.) + call realloc(xpmax_cellmask, maxpoly, keepExisting=.false.) + call realloc(ypmin_cellmask, maxpoly, keepExisting=.false.) + call realloc(ypmax_cellmask, maxpoly, keepExisting=.false.) + call realloc(iistart_cellmask, maxpoly, keepExisting=.false.) + call realloc(iiend_cellmask, maxpoly, keepExisting=.false.) + call realloc(zpl_start_cellmask, maxpoly, keepExisting=.false.) - ! generate cell mask - in = -1 - do k = 1, nump -! check if cell is in any "zpl>0" polygons - call dbpinpol_optinside_perpol(xzw(k), yzw(k), 0, 1, in, num, dmiss, JINS, NPL, xpl, ypl, zpl) + ipoint = 1 + ipoly = 0 - if (in == 0) then -! check if cell is outside all "zpl<0" polygons (enclosure) - call dbpinpol_optinside_perpol(xzw(k), yzw(k), 0, -1, in, num, dmiss, JINS, NPL, xpl, ypl, zpl) ! in=0: outside all "-1"-pol - if (num > 0) in = 1 - in ! only if at least one "-1"-type polygon was encountered + do while (ipoint < NPL) + ipoly = ipoly + 1 + if (ipoly > maxpoly) then + maxpoly = ceiling(maxpoly * 1.1) + call realloc(xpmin_cellmask, maxpoly, keepExisting=.true.) + call realloc(xpmax_cellmask, maxpoly, keepExisting=.true.) + call realloc(ypmin_cellmask, maxpoly, keepExisting=.true.) + call realloc(ypmax_cellmask, maxpoly, keepExisting=.true.) + call realloc(iistart_cellmask, maxpoly, keepExisting=.true.) + call realloc(iiend_cellmask, maxpoly, keepExisting=.true.) + call realloc(zpl_start_cellmask, maxpoly, keepExisting=.true.) end if - if (in > 0) then -! mask cell - cellmask(k) = 1 + ! Get polygon start and end pointer + call get_startend(NPL - ipoint + 1, xpl(ipoint:NPL), ypl(ipoint:NPL), istart, iend, dmiss) + istart = istart + ipoint - 1 + iend = iend + ipoint - 1 + + if (istart >= iend .or. iend > NPL) exit + + ! Store bounding box + xpmin_cellmask(ipoly) = minval(xpl(istart:iend)) + xpmax_cellmask(ipoly) = maxval(xpl(istart:iend)) + ypmin_cellmask(ipoly) = minval(ypl(istart:iend)) + ypmax_cellmask(ipoly) = maxval(ypl(istart:iend)) + + iistart_cellmask(ipoly) = istart + iiend_cellmask(ipoly) = iend + + ! Store polygon type (sign of zpl) + if (zpl(istart) /= dmiss) then + zpl_start_cellmask(ipoly) = zpl(istart) + else + zpl_start_cellmask(ipoly) = dmiss end if - if (mod(k, KMOD) == 0) then - call READYY(' ', min(1.0_dp, real(k, kind=dp) / nump)) + + ! Advance pointer + ipoint = iend + 2 + end do + + Npoly_cellmask = ipoly + cellmask_initialized = .true. + + end subroutine dbpinpol_cellmask_init + + !> Check if point should be masked. + elemental function dbpinpol_cellmask(xp, yp) result(mask) + use m_missing, only: dmiss !< Use globals directly + implicit none + + integer :: mask + real(kind=dp), intent(in) :: xp, yp + + integer :: istart, iend, ipoly, in_test + logical :: found_inside_enclosure + integer :: num_enclosures + + mask = 0 + if (.not. cellmask_initialized) return + + found_inside_enclosure = .false. + num_enclosures = 0 + + do ipoly = 1, Npoly_cellmask + istart = iistart_cellmask(ipoly) + iend = iiend_cellmask(ipoly) + + if (zpl_start_cellmask(ipoly) == dmiss) cycle + + ! Bounding box check + if (xp < xpmin_cellmask(ipoly) .or. xp > xpmax_cellmask(ipoly) .or. & + yp < ypmin_cellmask(ipoly) .or. yp > ypmax_cellmask(ipoly)) cycle + + ! Point-in-polygon test using globals dmiss and JINS + !call pinpok(xp, yp, iend - istart + 1, & + ! xpl_cellmask(istart:iend), ypl_cellmask(istart:iend), & + ! in_test, JINS, dmiss) + in_test = pinpok_elemental(xp, yp, iend - istart + 1) + + if (zpl_start_cellmask(ipoly) > 0.0_dp) then + ! Dry point polygon + if (in_test == 1) then + mask = 1 + return + end if + else if (zpl_start_cellmask(ipoly) < 0.0_dp) then + ! Enclosure polygon + num_enclosures = num_enclosures + 1 + if (in_test == 1) found_inside_enclosure = .true. end if end do - call READYY(' ', -1.0_dp) + + ! Outside all enclosures = dry + if (num_enclosures > 0 .and. .not. found_inside_enclosure) mask = 1 + + end function dbpinpol_cellmask + + !> Clean up module-level cellmask polygon data structures. + subroutine dbpinpol_cellmask_cleanup() + implicit none + + if (allocated(xpmin_cellmask)) deallocate(xpmin_cellmask) + if (allocated(xpmax_cellmask)) deallocate(xpmax_cellmask) + if (allocated(ypmin_cellmask)) deallocate(ypmin_cellmask) + if (allocated(ypmax_cellmask)) deallocate(ypmax_cellmask) + if (allocated(zpl_start_cellmask)) deallocate(zpl_start_cellmask) + if (allocated(iistart_cellmask)) deallocate(iistart_cellmask) + if (allocated(iiend_cellmask)) deallocate(iiend_cellmask) + + Npoly_cellmask = 0 + cellmask_initialized = .false. return - end subroutine pol_to_cellmask + end subroutine dbpinpol_cellmask_cleanup + + !> Optimized elemental point-in-polygon test using ray casting algorithm. + !! Accesses polygon data via module arrays. + elemental function pinpok_elemental(xl, yl, ipoly) result(inside) + use m_missing, only: dmiss, JINS !< Use globals directly + + implicit none + + real(kind=dp), intent(in) :: xl, yl !< Point coordinates (scalar) + integer, intent(in) :: ipoly !< Polygon index + integer :: inside !< Result: 1=inside, 0=outside + + integer :: i, j, istart, iend, crossings + real(kind=dp) :: x1, x2, y1, y2, xinters + + inside = 0 + + ! Get polygon bounds from module variables + istart = iistart_cellmask(ipoly) + iend = iiend_cellmask(ipoly) + + if (iend - istart + 1 <= 2) then + inside = 1 + goto 999 + end if + + ! Ray-casting algorithm + crossings = 0 + j = iend + + do i = istart, iend + if (xpl_cellmask(i) == dmiss) exit + + x1 = xpl_cellmask(j) + y1 = ypl_cellmask(j) + x2 = xpl_cellmask(i) + y2 = ypl_cellmask(i) + + ! Check if point is on vertex + if (xl == x1 .and. yl == y1) then + inside = 1 + goto 999 + end if + + ! Check if ray crosses edge + if ((y1 > yl) .neqv. (y2 > yl)) then + xinters = x1 + (yl - y1) * (x2 - x1) / (y2 - y1) + + if (xl < xinters) then + crossings = crossings + 1 + else if (xl == xinters) then + inside = 1 + goto 999 + end if + end if + + j = i + end do + + inside = mod(crossings, 2) + 999 continue + if (jins == 0) inside = 1 - inside + + end function pinpok_elemental end module m_pol_to_cellmask From 750cd9eaf4cf8e439d336b6dcc95e158589f8d61 Mon Sep 17 00:00:00 2001 From: FlorisBuwaldaDeltares Date: Mon, 1 Dec 2025 14:16:18 +0100 Subject: [PATCH 02/32] missing use statement --- .../src/dflowfm_kernel/prepost/pol_to_cellmask.f90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.f90 index 3fbf644f4c..99f6087829 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.f90 +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.f90 @@ -79,6 +79,7 @@ end subroutine pol_to_cellmask subroutine dbpinpol_cellmask_init(NPL, xpl, ypl, zpl) use m_alloc use m_missing, only: dp, dmiss + use geometry_module, only: get_startend implicit none From fad6a555ab8e7fa63e801810fe1f000b2139f13d Mon Sep 17 00:00:00 2001 From: FlorisBuwaldaDeltares Date: Mon, 1 Dec 2025 15:24:31 +0100 Subject: [PATCH 03/32] parallellize dbpinpol_cellmask --- .../prepost/pol_to_cellmask.f90 | 29 +++++++++++++------ 1 file changed, 20 insertions(+), 9 deletions(-) diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.f90 index 99f6087829..850f7eba28 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.f90 +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.f90 @@ -56,7 +56,12 @@ module m_pol_to_cellmask subroutine pol_to_cellmask() use network_data, only: cellmask, nump1d2d, npl, nump, xzw, yzw, xpl, ypl, zpl - + use m_partitioninfo, only: jampi +#ifdef _OPENMP + use omp_lib + integer :: temp_threads +#endif + integer :: k if (allocated(cellmask)) deallocate(cellmask) allocate(cellmask(nump1d2d)) cellmask = 0 @@ -65,9 +70,20 @@ subroutine pol_to_cellmask() ! Initialize once call dbpinpol_cellmask_init(NPL, xpl, ypl, zpl) - - ! Vectorized assignment - beautiful and simple! - cellmask(1:nump) = dbpinpol_cellmask(xzw(1:nump), yzw(1:nump)) +#ifdef _OPENMP + temp_threads = omp_get_max_threads() !> Save old number of threads + if (jampi == 0) then + call omp_set_num_threads(OMP_GET_NUM_PROCS()) !> Set number of threads to max for this O(N^2) operation + end if !> no else, in MPI mode omp num threads is already set to 1 +#endif + !$OMP PARALLEL DO SCHEDULE(DYNAMIC, 100) + do k = 1, nump + cellmask(k) = dbpinpol_cellmask(xzw(k), yzw(k)) + end do + !$OMP END PARALLEL DO +#ifdef _OPENMP + call omp_set_num_threads(temp_threads) +#endif ! Cleanup call dbpinpol_cellmask_cleanup() @@ -184,9 +200,6 @@ elemental function dbpinpol_cellmask(xp, yp) result(mask) yp < ypmin_cellmask(ipoly) .or. yp > ypmax_cellmask(ipoly)) cycle ! Point-in-polygon test using globals dmiss and JINS - !call pinpok(xp, yp, iend - istart + 1, & - ! xpl_cellmask(istart:iend), ypl_cellmask(istart:iend), & - ! in_test, JINS, dmiss) in_test = pinpok_elemental(xp, yp, iend - istart + 1) if (zpl_start_cellmask(ipoly) > 0.0_dp) then @@ -222,7 +235,6 @@ subroutine dbpinpol_cellmask_cleanup() Npoly_cellmask = 0 cellmask_initialized = .false. - return end subroutine dbpinpol_cellmask_cleanup !> Optimized elemental point-in-polygon test using ray casting algorithm. @@ -279,7 +291,6 @@ elemental function pinpok_elemental(xl, yl, ipoly) result(inside) goto 999 end if end if - j = i end do From b8517279702b438628903bcd8410101b3c541f12 Mon Sep 17 00:00:00 2001 From: FlorisBuwaldaDeltares Date: Mon, 1 Dec 2025 15:51:11 +0100 Subject: [PATCH 04/32] rename 1 --- .../prepost/{pol_to_cellmask.f90 => pol_to_cellmask_rename.f90} | 0 1 file changed, 0 insertions(+), 0 deletions(-) rename src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/{pol_to_cellmask.f90 => pol_to_cellmask_rename.f90} (100%) diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask_rename.f90 similarity index 100% rename from src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.f90 rename to src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask_rename.f90 From 2a797cbf1fb0fa2d4d069d58ac804654a34cf847 Mon Sep 17 00:00:00 2001 From: FlorisBuwaldaDeltares Date: Mon, 1 Dec 2025 15:53:59 +0100 Subject: [PATCH 05/32] rename 2, to F90. Add preprocessor flag in cmake --- src/engines_gpl/dflowfm/packages/dflowfm_kernel/CMakeLists.txt | 1 + .../prepost/{pol_to_cellmask_rename.f90 => pol_to_cellmask.F90} | 0 2 files changed, 1 insertion(+) rename src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/{pol_to_cellmask_rename.f90 => pol_to_cellmask.F90} (100%) diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/CMakeLists.txt b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/CMakeLists.txt index 28b352f490..b4285546a8 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/CMakeLists.txt +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/CMakeLists.txt @@ -216,6 +216,7 @@ set_source_files_properties(${fortran_version_files} ${CMAKE_CURRENT_SOURCE_DIR}/${dflowfm_kernel_compute_path}/filter.F90 ${CMAKE_CURRENT_SOURCE_DIR}/${dflowfm_kernel_compute_waves_surfbeat_path}/xbeach_filefunctions.F90 ${CMAKE_CURRENT_SOURCE_DIR}/${dflowfm_kernel_prepost_path}/find1dcells.f90 + ${CMAKE_CURRENT_SOURCE_DIR}/${dflowfm_kernel_prepost_path}/pol_to_cellmask.F90 ${CMAKE_CURRENT_SOURCE_DIR}/${dflowfm_kernel_prepost_path}/flow_modelinit.F90 ${CMAKE_CURRENT_SOURCE_DIR}/${dflowfm_kernel_prepost_path}/inidat.F90 ${CMAKE_CURRENT_SOURCE_DIR}/${dflowfm_kernel_prepost_path}/print_help_commandline.F90 diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask_rename.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.F90 similarity index 100% rename from src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask_rename.f90 rename to src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.F90 From d36ed146a0b607f0bbfca8415c6d66f97dc1580b Mon Sep 17 00:00:00 2001 From: FlorisBuwaldaDeltares Date: Mon, 1 Dec 2025 15:55:14 +0100 Subject: [PATCH 06/32] add dynamic scheduling for incells check --- .../dflowfm_kernel/src/dflowfm_kernel/prepost/find1dcells.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/find1dcells.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/find1dcells.f90 index 4d4d500add..1ccf470714 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/find1dcells.f90 +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/find1dcells.f90 @@ -70,7 +70,7 @@ subroutine find1dcells() temp_threads = omp_get_max_threads() !> Save old number of threads call omp_set_num_threads(OMP_GET_NUM_PROCS()) !> Set number of threads to max for this O(N^2) operation #endif - !$OMP PARALLEL DO + !$OMP PARALLEL DO SCHEDULE(DYNAMIC, 100) do L = 1, NUML1D if (KN(1, L) /= 0 .and. kn(3, L) /= LINK_1D .and. kn(3, L) /= LINK_1D_MAINBRANCH) then call INCELLS(Xk(KN(1, L)), Yk(KN(1, L)), left_2D_cells(L)) From cdd426930d07a9b3b807b800e3cc727ec8a6d011 Mon Sep 17 00:00:00 2001 From: FlorisBuwaldaDeltares Date: Tue, 2 Dec 2025 16:11:23 +0100 Subject: [PATCH 07/32] small fixes --- .../prepost/pol_to_cellmask.F90 | 126 +++++++++--------- 1 file changed, 64 insertions(+), 62 deletions(-) diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.F90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.F90 index 850f7eba28..f802261b17 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.F90 +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.F90 @@ -46,8 +46,7 @@ module m_pol_to_cellmask ! Module-level variables for cellmask polygon checking real(kind=dp), allocatable :: xpmin_cellmask(:), ypmin_cellmask(:) real(kind=dp), allocatable :: xpmax_cellmask(:), ypmax_cellmask(:) - real(kind=dp), allocatable :: zpl_start_cellmask(:) - real(kind=dp), allocatable :: xpl_cellmask(:), ypl_cellmask(:) + real(kind=dp), allocatable :: xpl_cellmask(:), ypl_cellmask(:), zpl_cellmask(:) integer, allocatable :: iistart_cellmask(:), iiend_cellmask(:) integer :: Npoly_cellmask = 0 logical :: cellmask_initialized = .false. @@ -55,21 +54,21 @@ module m_pol_to_cellmask contains subroutine pol_to_cellmask() - use network_data, only: cellmask, nump1d2d, npl, nump, xzw, yzw, xpl, ypl, zpl - use m_partitioninfo, only: jampi + use network_data, only: cellmask, nump1d2d, npl, nump, xzw, yzw, xpl, ypl, zpl + use m_partitioninfo, only: jampi #ifdef _OPENMP use omp_lib integer :: temp_threads #endif - integer :: k - if (allocated(cellmask)) deallocate(cellmask) - allocate(cellmask(nump1d2d)) - cellmask = 0 + integer :: k + if (allocated(cellmask)) deallocate (cellmask) + allocate (cellmask(nump1d2d)) + cellmask = 0 - if (NPL == 0) return + if (NPL == 0) return - ! Initialize once - call dbpinpol_cellmask_init(NPL, xpl, ypl, zpl) + ! Initialize once + call dbpinpol_cellmask_init(NPL, xpl, ypl, zpl) #ifdef _OPENMP temp_threads = omp_get_max_threads() !> Save old number of threads if (jampi == 0) then @@ -85,12 +84,12 @@ subroutine pol_to_cellmask() call omp_set_num_threads(temp_threads) #endif - ! Cleanup - call dbpinpol_cellmask_cleanup() + ! Cleanup + call dbpinpol_cellmask_cleanup() -end subroutine pol_to_cellmask + end subroutine pol_to_cellmask - !> Initialize polygon data structures for cellmask checking. + !> Initialize polygon data structures for cellmask checking. !! Must be called once before using dbpinpol_cellmask. subroutine dbpinpol_cellmask_init(NPL, xpl, ypl, zpl) use m_alloc @@ -114,6 +113,10 @@ subroutine dbpinpol_cellmask_init(NPL, xpl, ypl, zpl) return end if + ! Allocate and copy full coordinate arrays + xpl_cellmask = xpl + ypl_cellmask = ypl + MAXPOLY = 1000 call realloc(xpmin_cellmask, maxpoly, keepExisting=.false.) call realloc(xpmax_cellmask, maxpoly, keepExisting=.false.) @@ -121,7 +124,7 @@ subroutine dbpinpol_cellmask_init(NPL, xpl, ypl, zpl) call realloc(ypmax_cellmask, maxpoly, keepExisting=.false.) call realloc(iistart_cellmask, maxpoly, keepExisting=.false.) call realloc(iiend_cellmask, maxpoly, keepExisting=.false.) - call realloc(zpl_start_cellmask, maxpoly, keepExisting=.false.) + call realloc(zpl_cellmask, maxpoly, keepExisting=.false.) ipoint = 1 ipoly = 0 @@ -136,7 +139,7 @@ subroutine dbpinpol_cellmask_init(NPL, xpl, ypl, zpl) call realloc(ypmax_cellmask, maxpoly, keepExisting=.true.) call realloc(iistart_cellmask, maxpoly, keepExisting=.true.) call realloc(iiend_cellmask, maxpoly, keepExisting=.true.) - call realloc(zpl_start_cellmask, maxpoly, keepExisting=.true.) + call realloc(zpl_cellmask, maxpoly, keepExisting=.true.) end if ! Get polygon start and end pointer @@ -154,34 +157,31 @@ subroutine dbpinpol_cellmask_init(NPL, xpl, ypl, zpl) iistart_cellmask(ipoly) = istart iiend_cellmask(ipoly) = iend - + ! Store polygon type (sign of zpl) - if (zpl(istart) /= dmiss) then - zpl_start_cellmask(ipoly) = zpl(istart) - else - zpl_start_cellmask(ipoly) = dmiss - end if + zpl_cellmask(ipoly) = zpl(istart) ! Advance pointer ipoint = iend + 2 end do - + Npoly_cellmask = ipoly cellmask_initialized = .true. end subroutine dbpinpol_cellmask_init - !> Check if point should be masked. + !> Check if point should be masked. elemental function dbpinpol_cellmask(xp, yp) result(mask) - use m_missing, only: dmiss !< Use globals directly + use m_missing, only: dmiss implicit none integer :: mask real(kind=dp), intent(in) :: xp, yp - integer :: istart, iend, ipoly, in_test + integer :: ipoly, in_test logical :: found_inside_enclosure integer :: num_enclosures + real(kind=dp) :: zpl_val mask = 0 if (.not. cellmask_initialized) return @@ -189,34 +189,36 @@ elemental function dbpinpol_cellmask(xp, yp) result(mask) found_inside_enclosure = .false. num_enclosures = 0 + ! Single loop over all polygons do ipoly = 1, Npoly_cellmask - istart = iistart_cellmask(ipoly) - iend = iiend_cellmask(ipoly) + zpl_val = zpl_cellmask(ipoly) - if (zpl_start_cellmask(ipoly) == dmiss) cycle - - ! Bounding box check + ! Bounding box check (applies to all polygon types) if (xp < xpmin_cellmask(ipoly) .or. xp > xpmax_cellmask(ipoly) .or. & yp < ypmin_cellmask(ipoly) .or. yp > ypmax_cellmask(ipoly)) cycle - ! Point-in-polygon test using globals dmiss and JINS - in_test = pinpok_elemental(xp, yp, iend - istart + 1) + ! Point-in-polygon test + in_test = pinpok_elemental(xp, yp, ipoly) - if (zpl_start_cellmask(ipoly) > 0.0_dp) then - ! Dry point polygon + if (zpl_val == dmiss .or. zpl_val > 0.0_dp) then + ! Dry point polygon (including dmiss): inside = masked if (in_test == 1) then mask = 1 - return + return ! Early exit end if - else if (zpl_start_cellmask(ipoly) < 0.0_dp) then - ! Enclosure polygon + else if (zpl_val < 0.0_dp) then + ! Enclosure polygon: track if we're inside any enclosure num_enclosures = num_enclosures + 1 - if (in_test == 1) found_inside_enclosure = .true. + if (in_test == 1) then + found_inside_enclosure = .true. + end if end if end do - ! Outside all enclosures = dry - if (num_enclosures > 0 .and. .not. found_inside_enclosure) mask = 1 + ! After checking all polygons: outside all enclosures = dry + if (num_enclosures > 0 .and. .not. found_inside_enclosure) then + mask = 1 + end if end function dbpinpol_cellmask @@ -224,14 +226,14 @@ end function dbpinpol_cellmask subroutine dbpinpol_cellmask_cleanup() implicit none - if (allocated(xpmin_cellmask)) deallocate(xpmin_cellmask) - if (allocated(xpmax_cellmask)) deallocate(xpmax_cellmask) - if (allocated(ypmin_cellmask)) deallocate(ypmin_cellmask) - if (allocated(ypmax_cellmask)) deallocate(ypmax_cellmask) - if (allocated(zpl_start_cellmask)) deallocate(zpl_start_cellmask) - if (allocated(iistart_cellmask)) deallocate(iistart_cellmask) - if (allocated(iiend_cellmask)) deallocate(iiend_cellmask) - + if (allocated(xpmin_cellmask)) deallocate (xpmin_cellmask) + if (allocated(xpmax_cellmask)) deallocate (xpmax_cellmask) + if (allocated(ypmin_cellmask)) deallocate (ypmin_cellmask) + if (allocated(ypmax_cellmask)) deallocate (ypmax_cellmask) + if (allocated(zpl_cellmask)) deallocate (zpl_cellmask) + if (allocated(iistart_cellmask)) deallocate (iistart_cellmask) + if (allocated(iiend_cellmask)) deallocate (iiend_cellmask) + Npoly_cellmask = 0 cellmask_initialized = .false. @@ -240,23 +242,23 @@ end subroutine dbpinpol_cellmask_cleanup !> Optimized elemental point-in-polygon test using ray casting algorithm. !! Accesses polygon data via module arrays. elemental function pinpok_elemental(xl, yl, ipoly) result(inside) - use m_missing, only: dmiss, JINS !< Use globals directly + use m_missing, only: dmiss, JINS !< Use globals directly implicit none - real(kind=dp), intent(in) :: xl, yl !< Point coordinates (scalar) - integer, intent(in) :: ipoly !< Polygon index - integer :: inside !< Result: 1=inside, 0=outside + real(kind=dp), intent(in) :: xl, yl !< Point coordinates (scalar) + integer, intent(in) :: ipoly !< Polygon index + integer :: inside !< Result: 1=inside, 0=outside integer :: i, j, istart, iend, crossings real(kind=dp) :: x1, x2, y1, y2, xinters inside = 0 - + ! Get polygon bounds from module variables istart = iistart_cellmask(ipoly) iend = iiend_cellmask(ipoly) - + if (iend - istart + 1 <= 2) then inside = 1 goto 999 @@ -265,25 +267,25 @@ elemental function pinpok_elemental(xl, yl, ipoly) result(inside) ! Ray-casting algorithm crossings = 0 j = iend - + do i = istart, iend if (xpl_cellmask(i) == dmiss) exit - + x1 = xpl_cellmask(j) y1 = ypl_cellmask(j) x2 = xpl_cellmask(i) y2 = ypl_cellmask(i) - + ! Check if point is on vertex if (xl == x1 .and. yl == y1) then inside = 1 goto 999 end if - + ! Check if ray crosses edge if ((y1 > yl) .neqv. (y2 > yl)) then xinters = x1 + (yl - y1) * (x2 - x1) / (y2 - y1) - + if (xl < xinters) then crossings = crossings + 1 else if (xl == xinters) then @@ -295,7 +297,7 @@ elemental function pinpok_elemental(xl, yl, ipoly) result(inside) end do inside = mod(crossings, 2) - 999 continue +999 continue if (jins == 0) inside = 1 - inside end function pinpok_elemental From 93f95fd70a83665f38e92907c4576b00e6b54c39 Mon Sep 17 00:00:00 2001 From: FlorisBuwaldaDeltares Date: Tue, 2 Dec 2025 17:19:49 +0100 Subject: [PATCH 08/32] restore nested polygon check --- .../prepost/pol_to_cellmask.F90 | 39 +++++++++++-------- 1 file changed, 22 insertions(+), 17 deletions(-) diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.F90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.F90 index f802261b17..2a1c53b2a5 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.F90 +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.F90 @@ -1,4 +1,4 @@ -!----- AGPL -------------------------------------------------------------------- + !----- AGPL -------------------------------------------------------------------- ! ! Copyright (C) Stichting Deltares, 2017-2025. ! @@ -89,8 +89,6 @@ subroutine pol_to_cellmask() end subroutine pol_to_cellmask - !> Initialize polygon data structures for cellmask checking. - !! Must be called once before using dbpinpol_cellmask. subroutine dbpinpol_cellmask_init(NPL, xpl, ypl, zpl) use m_alloc use m_missing, only: dp, dmiss @@ -142,14 +140,12 @@ subroutine dbpinpol_cellmask_init(NPL, xpl, ypl, zpl) call realloc(zpl_cellmask, maxpoly, keepExisting=.true.) end if - ! Get polygon start and end pointer call get_startend(NPL - ipoint + 1, xpl(ipoint:NPL), ypl(ipoint:NPL), istart, iend, dmiss) istart = istart + ipoint - 1 iend = iend + ipoint - 1 if (istart >= iend .or. iend > NPL) exit - ! Store bounding box xpmin_cellmask(ipoly) = minval(xpl(istart:iend)) xpmax_cellmask(ipoly) = maxval(xpl(istart:iend)) ypmin_cellmask(ipoly) = minval(ypl(istart:iend)) @@ -157,28 +153,27 @@ subroutine dbpinpol_cellmask_init(NPL, xpl, ypl, zpl) iistart_cellmask(ipoly) = istart iiend_cellmask(ipoly) = iend - - ! Store polygon type (sign of zpl) zpl_cellmask(ipoly) = zpl(istart) - ! Advance pointer ipoint = iend + 2 end do Npoly_cellmask = ipoly + cellmask_initialized = .true. end subroutine dbpinpol_cellmask_init - !> Check if point should be masked. + !> Check if point should be masked - OPTIMIZED VERSION elemental function dbpinpol_cellmask(xp, yp) result(mask) - use m_missing, only: dmiss + use m_missing, only: dmiss, JINS implicit none integer :: mask real(kind=dp), intent(in) :: xp, yp integer :: ipoly, in_test + integer :: count_drypoint logical :: found_inside_enclosure integer :: num_enclosures real(kind=dp) :: zpl_val @@ -186,14 +181,14 @@ elemental function dbpinpol_cellmask(xp, yp) result(mask) mask = 0 if (.not. cellmask_initialized) return - found_inside_enclosure = .false. num_enclosures = 0 + count_drypoint = 0 ! Single loop over all polygons do ipoly = 1, Npoly_cellmask zpl_val = zpl_cellmask(ipoly) - ! Bounding box check (applies to all polygon types) + ! Bounding box check if (xp < xpmin_cellmask(ipoly) .or. xp > xpmax_cellmask(ipoly) .or. & yp < ypmin_cellmask(ipoly) .or. yp > ypmax_cellmask(ipoly)) cycle @@ -201,13 +196,12 @@ elemental function dbpinpol_cellmask(xp, yp) result(mask) in_test = pinpok_elemental(xp, yp, ipoly) if (zpl_val == dmiss .or. zpl_val > 0.0_dp) then - ! Dry point polygon (including dmiss): inside = masked + ! Dry point polygon if (in_test == 1) then - mask = 1 - return ! Early exit + count_drypoint = count_drypoint + 1 end if else if (zpl_val < 0.0_dp) then - ! Enclosure polygon: track if we're inside any enclosure + ! Enclosure polygon num_enclosures = num_enclosures + 1 if (in_test == 1) then found_inside_enclosure = .true. @@ -215,7 +209,18 @@ elemental function dbpinpol_cellmask(xp, yp) result(mask) end if end do - ! After checking all polygons: outside all enclosures = dry + ! Apply odd-even rule only if counting was needed + if (JINS == 1) then + if (mod(count_drypoint, 2) == 1) then + mask = 1 + end if + else + if (mod(count_drypoint, 2) == 0) then + mask = 1 + end if + end if + + ! Override with enclosure logic if applicable if (num_enclosures > 0 .and. .not. found_inside_enclosure) then mask = 1 end if From 9b12dffed80ec70f3367c22db05824b5287b4622 Mon Sep 17 00:00:00 2001 From: FlorisBuwaldaDeltares Date: Wed, 3 Dec 2025 10:31:08 +0100 Subject: [PATCH 09/32] added pol_to_cellmask unit tests --- .../packages/dflowfm_kernel/test/CMakeLists.txt | 9 +++++++++ .../dflowfm_kernel/test/test_pol_to_cellmask.cpp | 4 ++++ .../dflowfm_kernel/test/test_pol_to_cellmask.f90 | 16 ++++++++++++++++ 3 files changed, 29 insertions(+) create mode 100644 src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/test_pol_to_cellmask.cpp create mode 100644 src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/test_pol_to_cellmask.f90 diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/CMakeLists.txt b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/CMakeLists.txt index 08986fcc3e..767e3fa31e 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/CMakeLists.txt +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/CMakeLists.txt @@ -15,3 +15,12 @@ f90twtest(test_3d_layer_administration LIBRARIES dflowfm_kernel VISUAL_STUDIO_PATH engines_gpl/dflowfm/test ) + +f90twtest(test_pol_to_cellmask + CFILES test_pol_to_cellmask.cpp + F90FILES test_pol_to_cellmask.f90 + F2HFILES ${CMAKE_CURRENT_SOURCE_DIR}/test_pol_to_cellmask.f90 + OUTDIR "${CMAKE_CURRENT_BINARY_DIR}" + LIBRARIES dflowfm_kernel + VISUAL_STUDIO_PATH engines_gpl/dflowfm/test +) diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/test_pol_to_cellmask.cpp b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/test_pol_to_cellmask.cpp new file mode 100644 index 0000000000..ae203f825d --- /dev/null +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/test_pol_to_cellmask.cpp @@ -0,0 +1,4 @@ +/* set the fpp where the tests are implemented using f90tw "directives" in the fortran. */ +#define CURRENTTESTFILE "test_pol_to_cellmask.f90.h" + +#include "f90tw_gtest.h" diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/test_pol_to_cellmask.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/test_pol_to_cellmask.f90 new file mode 100644 index 0000000000..f420a5059c --- /dev/null +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/test_pol_to_cellmask.f90 @@ -0,0 +1,16 @@ +module test_pol_to_cellmask + use assertions_gtest + implicit none + +contains + + !$f90tw TESTCODE(TEST, test_pol_to_cellmask, test_layer_height_water_level_consistency, + subroutine test_mixed_polygon() bind(C) + + call pol_to_cellmask() + ! Check results + call f90_expect_near(mask, expected_mask, tolerance, "mask does not match expected value") + + end subroutine test_mixed_polygon + !$f90tw) +end module test_pol_to_cellmask From 9ce81e69ed628d9c1d28b6fa02083d81df5a673c Mon Sep 17 00:00:00 2001 From: FlorisBuwaldaDeltares Date: Wed, 3 Dec 2025 11:26:09 +0100 Subject: [PATCH 10/32] new tests --- .../test/test_pol_to_cellmask.f90 | 169 +++++++++++++++++- 1 file changed, 164 insertions(+), 5 deletions(-) diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/test_pol_to_cellmask.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/test_pol_to_cellmask.f90 index f420a5059c..d839476f2a 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/test_pol_to_cellmask.f90 +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/test_pol_to_cellmask.f90 @@ -1,16 +1,175 @@ module test_pol_to_cellmask use assertions_gtest + use precision, only: dp + use m_missing, only: dmiss + use network_data, only: cellmask, nump1d2d, npl, nump, xzw, yzw, xpl, ypl, zpl + use m_pol_to_cellmask, only: dbpinpol_cellmask_init, dbpinpol_cellmask, dbpinpol_cellmask_cleanup implicit none contains - !$f90tw TESTCODE(TEST, test_pol_to_cellmask, test_layer_height_water_level_consistency, + !$f90tw TESTCODE(TEST, test_pol_to_cellmask, test_mixed_polygon, test_nested_drypoint_polygons, subroutine test_mixed_polygon() bind(C) + ! Test with both enclosure (-1) and dry point (1) polygons + integer :: i, k - call pol_to_cellmask() - ! Check results - call f90_expect_near(mask, expected_mask, tolerance, "mask does not match expected value") + ! Setup test grid: 5x5 grid of cells + nump = 25 + + if (allocated(xzw)) deallocate(xzw) + if (allocated(yzw)) deallocate(yzw) + if (allocated(cellmask)) deallocate(cellmask) + allocate(xzw(nump), yzw(nump), cellmask(nump)) + + ! Create 5x5 grid of cell centers + do i = 1, 25 + xzw(i) = mod(i-1, 5) * 10.0_dp + 5.0_dp ! x: 5, 15, 25, 35, 45 + yzw(i) = ((i-1) / 5) * 10.0_dp + 5.0_dp ! y: 5, 15, 25, 35, 45 + end do + + ! Setup polygons: + ! 1. Enclosure polygon (zpl=-1): covers cells at x=[0,30], y=[0,30] + ! 2. Dry point polygon (zpl=1): inside enclosure at x=[10,20], y=[10,20] + + npl = 15 ! 5 points for enclosure + separator + 5 points for dry point + separators + + if (allocated(xpl)) deallocate(xpl) + if (allocated(ypl)) deallocate(ypl) + if (allocated(zpl)) deallocate(zpl) + allocate(xpl(npl), ypl(npl), zpl(npl)) + + ! Enclosure polygon (rectangle from 0,0 to 30,30) with zpl=-1 + xpl(1) = 0.0_dp; ypl(1) = 0.0_dp; zpl(1) = -1.0_dp + xpl(2) = 30.0_dp; ypl(2) = 0.0_dp; zpl(2) = -1.0_dp + xpl(3) = 30.0_dp; ypl(3) = 30.0_dp; zpl(3) = -1.0_dp + xpl(4) = 0.0_dp; ypl(4) = 30.0_dp; zpl(4) = -1.0_dp + xpl(5) = 0.0_dp; ypl(5) = 0.0_dp; zpl(5) = -1.0_dp + + ! Separator + xpl(6) = dmiss; ypl(6) = dmiss; zpl(6) = dmiss + + ! Dry point polygon (rectangle from 10,10 to 20,20) with zpl=1 + xpl(7) = 10.0_dp; ypl(7) = 10.0_dp; zpl(7) = 1.0_dp + xpl(8) = 20.0_dp; ypl(8) = 10.0_dp; zpl(8) = 1.0_dp + xpl(9) = 20.0_dp; ypl(9) = 20.0_dp; zpl(9) = 1.0_dp + xpl(10) = 10.0_dp; ypl(10) = 20.0_dp; zpl(10) = 1.0_dp + xpl(11) = 10.0_dp; ypl(11) = 10.0_dp; zpl(11) = 1.0_dp + + ! Separators + xpl(12) = dmiss; ypl(12) = dmiss; zpl(12) = dmiss + xpl(13) = dmiss; ypl(13) = dmiss; zpl(13) = dmiss + xpl(14) = dmiss; ypl(14) = dmiss; zpl(14) = dmiss + xpl(15) = dmiss; ypl(15) = dmiss; zpl(15) = dmiss + + ! Initialize polygon data structures + call dbpinpol_cellmask_init(NPL, xpl, ypl, zpl) + + ! Process all cells + cellmask = 0 + cellmask = dbpinpol_cellmask(xzw, yzw) + + ! Cleanup + call dbpinpol_cellmask_cleanup() + + ! Check results: + ! Cell at (5,5) - inside enclosure, outside dry point -> mask=0 + call f90_expect_eq(cellmask(1), 0, "Cell (5,5) should not be masked") + + ! Cell at (15,15) - inside enclosure AND inside dry point -> mask=1 + call f90_expect_eq(cellmask(7), 1, "Cell (15,15) should be masked (dry point)") + + ! Cell at (35,5) - outside enclosure -> mask=1 + call f90_expect_eq(cellmask(4), 1, "Cell (35,5) should be masked (outside enclosure)") + + ! Cell at (25,25) - inside enclosure, outside dry point -> mask=0 + call f90_expect_eq(cellmask(13), 0, "Cell (25,25) should not be masked") + + ! Cleanup + deallocate(xzw, yzw, xpl, ypl, zpl, cellmask) end subroutine test_mixed_polygon + + subroutine test_nested_drypoint_polygons() bind(C) + ! Test nested dry point polygons (odd-even rule) + integer :: i, k + + ! Setup test grid: 5x5 grid of cells + nump = 25 + + if (allocated(xzw)) deallocate(xzw) + if (allocated(yzw)) deallocate(yzw) + if (allocated(cellmask)) deallocate(cellmask) + allocate(xzw(nump), yzw(nump), cellmask(nump)) + + ! Create 5x5 grid of cell centers + do i = 1, 25 + xzw(i) = mod(i-1, 5) * 10.0_dp + 5.0_dp + yzw(i) = ((i-1) / 5) * 10.0_dp + 5.0_dp + end do + + ! Setup nested dry point polygons: + ! 1. Outer dry point polygon: x=[0,40], y=[0,40] with zpl=1 + ! 2. Inner dry point polygon: x=[10,30], y=[10,30] with zpl=1 + ! Cells inside odd number of polygons are masked + + npl = 13 + + if (allocated(xpl)) deallocate(xpl) + if (allocated(ypl)) deallocate(ypl) + if (allocated(zpl)) deallocate(zpl) + allocate(xpl(npl), ypl(npl), zpl(npl)) + + ! Outer dry point polygon (rectangle from 0,0 to 40,40) + xpl(1) = 0.0_dp; ypl(1) = 0.0_dp; zpl(1) = 1.0_dp + xpl(2) = 40.0_dp; ypl(2) = 0.0_dp; zpl(2) = 1.0_dp + xpl(3) = 40.0_dp; ypl(3) = 40.0_dp; zpl(3) = 1.0_dp + xpl(4) = 0.0_dp; ypl(4) = 40.0_dp; zpl(4) = 1.0_dp + xpl(5) = 0.0_dp; ypl(5) = 0.0_dp; zpl(5) = 1.0_dp + + ! Separator + xpl(6) = dmiss; ypl(6) = dmiss; zpl(6) = dmiss + + ! Inner dry point polygon (rectangle from 10,10 to 30,30) + xpl(7) = 10.0_dp; ypl(7) = 10.0_dp; zpl(7) = 1.0_dp + xpl(8) = 30.0_dp; ypl(8) = 10.0_dp; zpl(8) = 1.0_dp + xpl(9) = 30.0_dp; ypl(9) = 30.0_dp; zpl(9) = 1.0_dp + xpl(10) = 10.0_dp; ypl(10) = 30.0_dp; zpl(10) = 1.0_dp + xpl(11) = 10.0_dp; ypl(11) = 10.0_dp; zpl(11) = 1.0_dp + + ! Separators + xpl(12) = dmiss; ypl(12) = dmiss; zpl(12) = dmiss + xpl(13) = dmiss; ypl(13) = dmiss; zpl(13) = dmiss + + ! Initialize polygon data structures + call dbpinpol_cellmask_init(NPL, xpl, ypl, zpl) + + ! Process all cells + cellmask = 0 + cellmask = dbpinpol_cellmask(xzw, yzw) + + ! Cleanup + call dbpinpol_cellmask_cleanup() + + ! Check results (odd-even rule): + ! Cell at (5,5) - inside 1 polygon (outer) -> mask=1 + call f90_expect_eq(cellmask(1), 1, "Cell (5,5) should be masked (inside 1 polygon)") + + ! Cell at (15,15) - inside 2 polygons (outer+inner) -> mask=0 + call f90_expect_eq(cellmask(7), 0, "Cell (15,15) should not be masked (inside 2 polygons)") + + ! Cell at (25,25) - inside 2 polygons (outer+inner) -> mask=0 + call f90_expect_eq(cellmask(13), 0, "Cell (25,25) should not be masked (inside 2 polygons)") + + ! Cell at (35,35) - inside 1 polygon (outer) -> mask=1 + call f90_expect_eq(cellmask(19), 1, "Cell (35,35) should be masked (inside 1 polygon)") + + ! Cell at (45,45) - outside all polygons -> mask=0 + call f90_expect_eq(cellmask(25), 0, "Cell (45,45) should not be masked (outside all)") + + ! Cleanup + deallocate(xzw, yzw, xpl, ypl, zpl, cellmask) + + end subroutine test_nested_drypoint_polygons !$f90tw) -end module test_pol_to_cellmask + +end module test_pol_to_cellmask \ No newline at end of file From de6f60065f4da5bab927a28ec591249fc9683ac2 Mon Sep 17 00:00:00 2001 From: FlorisBuwaldaDeltares Date: Wed, 3 Dec 2025 11:45:14 +0100 Subject: [PATCH 11/32] separate module --- .../prepost/pol_to_cellmask.F90 | 110 ++++++++++-------- .../test/test_pol_to_cellmask.f90 | 2 +- 2 files changed, 60 insertions(+), 52 deletions(-) diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.F90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.F90 index 2a1c53b2a5..16d9810b05 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.F90 +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.F90 @@ -27,21 +27,16 @@ ! !------------------------------------------------------------------------------- -! -! - -!> update cellmask from samples -!> a cell is dry when it is: -!> 1) inside ANY "1"-polygon (drypnts), OR -!> 2) outside ALL "-1"-polygons (enclosures) -module m_pol_to_cellmask - +module m_dbpinpol_cellmask + use m_missing, only: jins, dmiss use precision, only: dp + implicit none private - public :: pol_to_cellmask + !> dbpinpol routines are public to avoid PetSC dependency in unit tests + public :: dbpinpol_cellmask_init, dbpinpol_cellmask_cleanup, dbpinpol_cellmask ! Module-level variables for cellmask polygon checking real(kind=dp), allocatable :: xpmin_cellmask(:), ypmin_cellmask(:) @@ -53,45 +48,8 @@ module m_pol_to_cellmask contains - subroutine pol_to_cellmask() - use network_data, only: cellmask, nump1d2d, npl, nump, xzw, yzw, xpl, ypl, zpl - use m_partitioninfo, only: jampi -#ifdef _OPENMP - use omp_lib - integer :: temp_threads -#endif - integer :: k - if (allocated(cellmask)) deallocate (cellmask) - allocate (cellmask(nump1d2d)) - cellmask = 0 - - if (NPL == 0) return - - ! Initialize once - call dbpinpol_cellmask_init(NPL, xpl, ypl, zpl) -#ifdef _OPENMP - temp_threads = omp_get_max_threads() !> Save old number of threads - if (jampi == 0) then - call omp_set_num_threads(OMP_GET_NUM_PROCS()) !> Set number of threads to max for this O(N^2) operation - end if !> no else, in MPI mode omp num threads is already set to 1 -#endif - !$OMP PARALLEL DO SCHEDULE(DYNAMIC, 100) - do k = 1, nump - cellmask(k) = dbpinpol_cellmask(xzw(k), yzw(k)) - end do - !$OMP END PARALLEL DO -#ifdef _OPENMP - call omp_set_num_threads(temp_threads) -#endif - - ! Cleanup - call dbpinpol_cellmask_cleanup() - - end subroutine pol_to_cellmask - subroutine dbpinpol_cellmask_init(NPL, xpl, ypl, zpl) use m_alloc - use m_missing, only: dp, dmiss use geometry_module, only: get_startend implicit none @@ -166,8 +124,6 @@ end subroutine dbpinpol_cellmask_init !> Check if point should be masked - OPTIMIZED VERSION elemental function dbpinpol_cellmask(xp, yp) result(mask) - use m_missing, only: dmiss, JINS - implicit none integer :: mask real(kind=dp), intent(in) :: xp, yp @@ -229,7 +185,6 @@ end function dbpinpol_cellmask !> Clean up module-level cellmask polygon data structures. subroutine dbpinpol_cellmask_cleanup() - implicit none if (allocated(xpmin_cellmask)) deallocate (xpmin_cellmask) if (allocated(xpmax_cellmask)) deallocate (xpmax_cellmask) @@ -247,7 +202,6 @@ end subroutine dbpinpol_cellmask_cleanup !> Optimized elemental point-in-polygon test using ray casting algorithm. !! Accesses polygon data via module arrays. elemental function pinpok_elemental(xl, yl, ipoly) result(inside) - use m_missing, only: dmiss, JINS !< Use globals directly implicit none @@ -307,4 +261,58 @@ elemental function pinpok_elemental(xl, yl, ipoly) result(inside) end function pinpok_elemental +end module m_dbpinpol_cellmask + +!> update cellmask from samples +!> a cell is dry when it is: +!> 1) inside ANY "1"-polygon (drypnts), OR +!> 2) outside ALL "-1"-polygons (enclosures) +module m_pol_to_cellmask + + use precision, only: dp + use m_dbpinpol_cellmask, only: dbpinpol_cellmask_init, dbpinpol_cellmask_cleanup, dbpinpol_cellmask + implicit none + + private + + public :: pol_to_cellmask + +contains + + subroutine pol_to_cellmask() + use network_data, only: cellmask, nump1d2d, npl, nump, xzw, yzw, xpl, ypl, zpl + use m_partitioninfo, only: jampi +#ifdef _OPENMP + use omp_lib + integer :: temp_threads +#endif + integer :: k + if (allocated(cellmask)) deallocate (cellmask) + allocate (cellmask(nump1d2d)) + cellmask = 0 + + if (NPL == 0) return + + ! Initialize once + call dbpinpol_cellmask_init(NPL, xpl, ypl, zpl) +#ifdef _OPENMP + temp_threads = omp_get_max_threads() !> Save old number of threads + if (jampi == 0) then + call omp_set_num_threads(OMP_GET_NUM_PROCS()) !> Set number of threads to max for this O(N^2) operation + end if !> no else, in MPI mode omp num threads is already set to 1 +#endif + !$OMP PARALLEL DO SCHEDULE(DYNAMIC, 100) + do k = 1, nump + cellmask(k) = dbpinpol_cellmask(xzw(k), yzw(k)) + end do + !$OMP END PARALLEL DO +#ifdef _OPENMP + call omp_set_num_threads(temp_threads) +#endif + + ! Cleanup + call dbpinpol_cellmask_cleanup() + + end subroutine pol_to_cellmask + end module m_pol_to_cellmask diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/test_pol_to_cellmask.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/test_pol_to_cellmask.f90 index d839476f2a..8216d4b1c0 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/test_pol_to_cellmask.f90 +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/test_pol_to_cellmask.f90 @@ -3,7 +3,7 @@ module test_pol_to_cellmask use precision, only: dp use m_missing, only: dmiss use network_data, only: cellmask, nump1d2d, npl, nump, xzw, yzw, xpl, ypl, zpl - use m_pol_to_cellmask, only: dbpinpol_cellmask_init, dbpinpol_cellmask, dbpinpol_cellmask_cleanup + use m_dbpinpol_cellmask, only: dbpinpol_cellmask_init, dbpinpol_cellmask, dbpinpol_cellmask_cleanup implicit none contains From e3013bf91ebbfeafe66fa0de3f8051e52801b825 Mon Sep 17 00:00:00 2001 From: FlorisBuwaldaDeltares Date: Wed, 3 Dec 2025 12:06:21 +0100 Subject: [PATCH 12/32] separate module to remove PetSC dependency attempt 2 --- .../prepost/dbpinpol_cellmask.f90 | 235 +++++++++++++++++ .../prepost/delete_drypoints_from_netgeom.f90 | 4 +- .../prepost/pol_to_cellmask.F90 | 236 ------------------ 3 files changed, 237 insertions(+), 238 deletions(-) create mode 100644 src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 new file mode 100644 index 0000000000..1a62482632 --- /dev/null +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 @@ -0,0 +1,235 @@ +module m_dbpinpol_cellmask + use m_missing, only: jins, dmiss + use precision, only: dp + + implicit none + + private + + !> dbpinpol routines are public to avoid PetSC dependency in unit tests + public :: dbpinpol_cellmask_init, dbpinpol_cellmask_cleanup, dbpinpol_cellmask + + ! Module-level variables for cellmask polygon checking + real(kind=dp), allocatable :: xpmin_cellmask(:), ypmin_cellmask(:) + real(kind=dp), allocatable :: xpmax_cellmask(:), ypmax_cellmask(:) + real(kind=dp), allocatable :: xpl_cellmask(:), ypl_cellmask(:), zpl_cellmask(:) + integer, allocatable :: iistart_cellmask(:), iiend_cellmask(:) + integer :: Npoly_cellmask = 0 + logical :: cellmask_initialized = .false. + +contains + + subroutine dbpinpol_cellmask_init(NPL, xpl, ypl, zpl) + use m_alloc + use geometry_module, only: get_startend + + implicit none + + integer, intent(in) :: NPL + real(kind=dp), intent(in) :: xpl(NPL), ypl(NPL), zpl(NPL) + + integer :: MAXPOLY + integer :: ipoint, istart, iend, ipoly + + if (cellmask_initialized) then + call dbpinpol_cellmask_cleanup() + end if + + if (NPL == 0) then + cellmask_initialized = .true. + return + end if + + ! Allocate and copy full coordinate arrays + xpl_cellmask = xpl + ypl_cellmask = ypl + + MAXPOLY = 1000 + call realloc(xpmin_cellmask, maxpoly, keepExisting=.false.) + call realloc(xpmax_cellmask, maxpoly, keepExisting=.false.) + call realloc(ypmin_cellmask, maxpoly, keepExisting=.false.) + call realloc(ypmax_cellmask, maxpoly, keepExisting=.false.) + call realloc(iistart_cellmask, maxpoly, keepExisting=.false.) + call realloc(iiend_cellmask, maxpoly, keepExisting=.false.) + call realloc(zpl_cellmask, maxpoly, keepExisting=.false.) + + ipoint = 1 + ipoly = 0 + + do while (ipoint < NPL) + ipoly = ipoly + 1 + if (ipoly > maxpoly) then + maxpoly = ceiling(maxpoly * 1.1) + call realloc(xpmin_cellmask, maxpoly, keepExisting=.true.) + call realloc(xpmax_cellmask, maxpoly, keepExisting=.true.) + call realloc(ypmin_cellmask, maxpoly, keepExisting=.true.) + call realloc(ypmax_cellmask, maxpoly, keepExisting=.true.) + call realloc(iistart_cellmask, maxpoly, keepExisting=.true.) + call realloc(iiend_cellmask, maxpoly, keepExisting=.true.) + call realloc(zpl_cellmask, maxpoly, keepExisting=.true.) + end if + + call get_startend(NPL - ipoint + 1, xpl(ipoint:NPL), ypl(ipoint:NPL), istart, iend, dmiss) + istart = istart + ipoint - 1 + iend = iend + ipoint - 1 + + if (istart >= iend .or. iend > NPL) exit + + xpmin_cellmask(ipoly) = minval(xpl(istart:iend)) + xpmax_cellmask(ipoly) = maxval(xpl(istart:iend)) + ypmin_cellmask(ipoly) = minval(ypl(istart:iend)) + ypmax_cellmask(ipoly) = maxval(ypl(istart:iend)) + + iistart_cellmask(ipoly) = istart + iiend_cellmask(ipoly) = iend + zpl_cellmask(ipoly) = zpl(istart) + + ipoint = iend + 2 + end do + + Npoly_cellmask = ipoly + + cellmask_initialized = .true. + + end subroutine dbpinpol_cellmask_init + + !> Check if point should be masked - OPTIMIZED VERSION + elemental function dbpinpol_cellmask(xp, yp) result(mask) + + integer :: mask + real(kind=dp), intent(in) :: xp, yp + + integer :: ipoly, in_test + integer :: count_drypoint + logical :: found_inside_enclosure + integer :: num_enclosures + real(kind=dp) :: zpl_val + + mask = 0 + if (.not. cellmask_initialized) return + + num_enclosures = 0 + count_drypoint = 0 + + ! Single loop over all polygons + do ipoly = 1, Npoly_cellmask + zpl_val = zpl_cellmask(ipoly) + + ! Bounding box check + if (xp < xpmin_cellmask(ipoly) .or. xp > xpmax_cellmask(ipoly) .or. & + yp < ypmin_cellmask(ipoly) .or. yp > ypmax_cellmask(ipoly)) cycle + + ! Point-in-polygon test + in_test = pinpok_elemental(xp, yp, ipoly) + + if (zpl_val == dmiss .or. zpl_val > 0.0_dp) then + ! Dry point polygon + if (in_test == 1) then + count_drypoint = count_drypoint + 1 + end if + else if (zpl_val < 0.0_dp) then + ! Enclosure polygon + num_enclosures = num_enclosures + 1 + if (in_test == 1) then + found_inside_enclosure = .true. + end if + end if + end do + + ! Apply odd-even rule only if counting was needed + if (JINS == 1) then + if (mod(count_drypoint, 2) == 1) then + mask = 1 + end if + else + if (mod(count_drypoint, 2) == 0) then + mask = 1 + end if + end if + + ! Override with enclosure logic if applicable + if (num_enclosures > 0 .and. .not. found_inside_enclosure) then + mask = 1 + end if + + end function dbpinpol_cellmask + + !> Clean up module-level cellmask polygon data structures. + subroutine dbpinpol_cellmask_cleanup() + + if (allocated(xpmin_cellmask)) deallocate (xpmin_cellmask) + if (allocated(xpmax_cellmask)) deallocate (xpmax_cellmask) + if (allocated(ypmin_cellmask)) deallocate (ypmin_cellmask) + if (allocated(ypmax_cellmask)) deallocate (ypmax_cellmask) + if (allocated(zpl_cellmask)) deallocate (zpl_cellmask) + if (allocated(iistart_cellmask)) deallocate (iistart_cellmask) + if (allocated(iiend_cellmask)) deallocate (iiend_cellmask) + + Npoly_cellmask = 0 + cellmask_initialized = .false. + + end subroutine dbpinpol_cellmask_cleanup + + !> Optimized elemental point-in-polygon test using ray casting algorithm. + !! Accesses polygon data via module arrays. + elemental function pinpok_elemental(xl, yl, ipoly) result(inside) + + implicit none + + real(kind=dp), intent(in) :: xl, yl !< Point coordinates (scalar) + integer, intent(in) :: ipoly !< Polygon index + integer :: inside !< Result: 1=inside, 0=outside + + integer :: i, j, istart, iend, crossings + real(kind=dp) :: x1, x2, y1, y2, xinters + + inside = 0 + + ! Get polygon bounds from module variables + istart = iistart_cellmask(ipoly) + iend = iiend_cellmask(ipoly) + + if (iend - istart + 1 <= 2) then + inside = 1 + goto 999 + end if + + ! Ray-casting algorithm + crossings = 0 + j = iend + + do i = istart, iend + if (xpl_cellmask(i) == dmiss) exit + + x1 = xpl_cellmask(j) + y1 = ypl_cellmask(j) + x2 = xpl_cellmask(i) + y2 = ypl_cellmask(i) + + ! Check if point is on vertex + if (xl == x1 .and. yl == y1) then + inside = 1 + goto 999 + end if + + ! Check if ray crosses edge + if ((y1 > yl) .neqv. (y2 > yl)) then + xinters = x1 + (yl - y1) * (x2 - x1) / (y2 - y1) + + if (xl < xinters) then + crossings = crossings + 1 + else if (xl == xinters) then + inside = 1 + goto 999 + end if + end if + j = i + end do + + inside = mod(crossings, 2) +999 continue + if (jins == 0) inside = 1 - inside + + end function pinpok_elemental + +end module m_dbpinpol_cellmask \ No newline at end of file diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/delete_drypoints_from_netgeom.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/delete_drypoints_from_netgeom.f90 index f50d4a6a4d..696210dd7d 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/delete_drypoints_from_netgeom.f90 +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/delete_drypoints_from_netgeom.f90 @@ -34,7 +34,7 @@ !! Grid enclosures are handled via the jinside=-1 option. module m_delete_drypoints_from_netgeom use m_remove_masked_netcells, only: remove_masked_netcells - use m_pol_to_cellmask, only: pol_to_cellmask + !use m_pol_to_cellmask, only: pol_to_cellmask use m_fix_global_polygons, only: fix_global_polygons implicit none @@ -154,7 +154,7 @@ subroutine delete_drypoints_from_netgeom(dryptsfilelist, jaconfirm, jinside, upd call fix_global_polygons(1, 0) end if - call pol_to_cellmask() ! third column in pol-file may be used to specify inside (1), or outside (0) mode, only 0 or 1 allowed. + !call pol_to_cellmask() ! third column in pol-file may be used to specify inside (1), or outside (0) mode, only 0 or 1 allowed. call delpol() call restorepol() diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.F90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.F90 index 16d9810b05..ef09524ffc 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.F90 +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.F90 @@ -27,242 +27,6 @@ ! !------------------------------------------------------------------------------- -module m_dbpinpol_cellmask - use m_missing, only: jins, dmiss - use precision, only: dp - - implicit none - - private - - !> dbpinpol routines are public to avoid PetSC dependency in unit tests - public :: dbpinpol_cellmask_init, dbpinpol_cellmask_cleanup, dbpinpol_cellmask - - ! Module-level variables for cellmask polygon checking - real(kind=dp), allocatable :: xpmin_cellmask(:), ypmin_cellmask(:) - real(kind=dp), allocatable :: xpmax_cellmask(:), ypmax_cellmask(:) - real(kind=dp), allocatable :: xpl_cellmask(:), ypl_cellmask(:), zpl_cellmask(:) - integer, allocatable :: iistart_cellmask(:), iiend_cellmask(:) - integer :: Npoly_cellmask = 0 - logical :: cellmask_initialized = .false. - -contains - - subroutine dbpinpol_cellmask_init(NPL, xpl, ypl, zpl) - use m_alloc - use geometry_module, only: get_startend - - implicit none - - integer, intent(in) :: NPL - real(kind=dp), intent(in) :: xpl(NPL), ypl(NPL), zpl(NPL) - - integer :: MAXPOLY - integer :: ipoint, istart, iend, ipoly - - if (cellmask_initialized) then - call dbpinpol_cellmask_cleanup() - end if - - if (NPL == 0) then - cellmask_initialized = .true. - return - end if - - ! Allocate and copy full coordinate arrays - xpl_cellmask = xpl - ypl_cellmask = ypl - - MAXPOLY = 1000 - call realloc(xpmin_cellmask, maxpoly, keepExisting=.false.) - call realloc(xpmax_cellmask, maxpoly, keepExisting=.false.) - call realloc(ypmin_cellmask, maxpoly, keepExisting=.false.) - call realloc(ypmax_cellmask, maxpoly, keepExisting=.false.) - call realloc(iistart_cellmask, maxpoly, keepExisting=.false.) - call realloc(iiend_cellmask, maxpoly, keepExisting=.false.) - call realloc(zpl_cellmask, maxpoly, keepExisting=.false.) - - ipoint = 1 - ipoly = 0 - - do while (ipoint < NPL) - ipoly = ipoly + 1 - if (ipoly > maxpoly) then - maxpoly = ceiling(maxpoly * 1.1) - call realloc(xpmin_cellmask, maxpoly, keepExisting=.true.) - call realloc(xpmax_cellmask, maxpoly, keepExisting=.true.) - call realloc(ypmin_cellmask, maxpoly, keepExisting=.true.) - call realloc(ypmax_cellmask, maxpoly, keepExisting=.true.) - call realloc(iistart_cellmask, maxpoly, keepExisting=.true.) - call realloc(iiend_cellmask, maxpoly, keepExisting=.true.) - call realloc(zpl_cellmask, maxpoly, keepExisting=.true.) - end if - - call get_startend(NPL - ipoint + 1, xpl(ipoint:NPL), ypl(ipoint:NPL), istart, iend, dmiss) - istart = istart + ipoint - 1 - iend = iend + ipoint - 1 - - if (istart >= iend .or. iend > NPL) exit - - xpmin_cellmask(ipoly) = minval(xpl(istart:iend)) - xpmax_cellmask(ipoly) = maxval(xpl(istart:iend)) - ypmin_cellmask(ipoly) = minval(ypl(istart:iend)) - ypmax_cellmask(ipoly) = maxval(ypl(istart:iend)) - - iistart_cellmask(ipoly) = istart - iiend_cellmask(ipoly) = iend - zpl_cellmask(ipoly) = zpl(istart) - - ipoint = iend + 2 - end do - - Npoly_cellmask = ipoly - - cellmask_initialized = .true. - - end subroutine dbpinpol_cellmask_init - - !> Check if point should be masked - OPTIMIZED VERSION - elemental function dbpinpol_cellmask(xp, yp) result(mask) - - integer :: mask - real(kind=dp), intent(in) :: xp, yp - - integer :: ipoly, in_test - integer :: count_drypoint - logical :: found_inside_enclosure - integer :: num_enclosures - real(kind=dp) :: zpl_val - - mask = 0 - if (.not. cellmask_initialized) return - - num_enclosures = 0 - count_drypoint = 0 - - ! Single loop over all polygons - do ipoly = 1, Npoly_cellmask - zpl_val = zpl_cellmask(ipoly) - - ! Bounding box check - if (xp < xpmin_cellmask(ipoly) .or. xp > xpmax_cellmask(ipoly) .or. & - yp < ypmin_cellmask(ipoly) .or. yp > ypmax_cellmask(ipoly)) cycle - - ! Point-in-polygon test - in_test = pinpok_elemental(xp, yp, ipoly) - - if (zpl_val == dmiss .or. zpl_val > 0.0_dp) then - ! Dry point polygon - if (in_test == 1) then - count_drypoint = count_drypoint + 1 - end if - else if (zpl_val < 0.0_dp) then - ! Enclosure polygon - num_enclosures = num_enclosures + 1 - if (in_test == 1) then - found_inside_enclosure = .true. - end if - end if - end do - - ! Apply odd-even rule only if counting was needed - if (JINS == 1) then - if (mod(count_drypoint, 2) == 1) then - mask = 1 - end if - else - if (mod(count_drypoint, 2) == 0) then - mask = 1 - end if - end if - - ! Override with enclosure logic if applicable - if (num_enclosures > 0 .and. .not. found_inside_enclosure) then - mask = 1 - end if - - end function dbpinpol_cellmask - - !> Clean up module-level cellmask polygon data structures. - subroutine dbpinpol_cellmask_cleanup() - - if (allocated(xpmin_cellmask)) deallocate (xpmin_cellmask) - if (allocated(xpmax_cellmask)) deallocate (xpmax_cellmask) - if (allocated(ypmin_cellmask)) deallocate (ypmin_cellmask) - if (allocated(ypmax_cellmask)) deallocate (ypmax_cellmask) - if (allocated(zpl_cellmask)) deallocate (zpl_cellmask) - if (allocated(iistart_cellmask)) deallocate (iistart_cellmask) - if (allocated(iiend_cellmask)) deallocate (iiend_cellmask) - - Npoly_cellmask = 0 - cellmask_initialized = .false. - - end subroutine dbpinpol_cellmask_cleanup - - !> Optimized elemental point-in-polygon test using ray casting algorithm. - !! Accesses polygon data via module arrays. - elemental function pinpok_elemental(xl, yl, ipoly) result(inside) - - implicit none - - real(kind=dp), intent(in) :: xl, yl !< Point coordinates (scalar) - integer, intent(in) :: ipoly !< Polygon index - integer :: inside !< Result: 1=inside, 0=outside - - integer :: i, j, istart, iend, crossings - real(kind=dp) :: x1, x2, y1, y2, xinters - - inside = 0 - - ! Get polygon bounds from module variables - istart = iistart_cellmask(ipoly) - iend = iiend_cellmask(ipoly) - - if (iend - istart + 1 <= 2) then - inside = 1 - goto 999 - end if - - ! Ray-casting algorithm - crossings = 0 - j = iend - - do i = istart, iend - if (xpl_cellmask(i) == dmiss) exit - - x1 = xpl_cellmask(j) - y1 = ypl_cellmask(j) - x2 = xpl_cellmask(i) - y2 = ypl_cellmask(i) - - ! Check if point is on vertex - if (xl == x1 .and. yl == y1) then - inside = 1 - goto 999 - end if - - ! Check if ray crosses edge - if ((y1 > yl) .neqv. (y2 > yl)) then - xinters = x1 + (yl - y1) * (x2 - x1) / (y2 - y1) - - if (xl < xinters) then - crossings = crossings + 1 - else if (xl == xinters) then - inside = 1 - goto 999 - end if - end if - j = i - end do - - inside = mod(crossings, 2) -999 continue - if (jins == 0) inside = 1 - inside - - end function pinpok_elemental - -end module m_dbpinpol_cellmask - !> update cellmask from samples !> a cell is dry when it is: !> 1) inside ANY "1"-polygon (drypnts), OR From 4240b9d1f4f4ade252d25eac413ec84ece424a08 Mon Sep 17 00:00:00 2001 From: FlorisBuwaldaDeltares Date: Wed, 3 Dec 2025 12:55:35 +0100 Subject: [PATCH 13/32] fix test, new global "enclosures_present" --- .../prepost/dbpinpol_cellmask.f90 | 19 ++++++++++--------- .../test/test_pol_to_cellmask.f90 | 4 +++- 2 files changed, 13 insertions(+), 10 deletions(-) diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 index 1a62482632..772bfde968 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 @@ -16,6 +16,7 @@ module m_dbpinpol_cellmask integer, allocatable :: iistart_cellmask(:), iiend_cellmask(:) integer :: Npoly_cellmask = 0 logical :: cellmask_initialized = .false. + logical :: enclosures_present = .false. contains @@ -89,6 +90,8 @@ subroutine dbpinpol_cellmask_init(NPL, xpl, ypl, zpl) Npoly_cellmask = ipoly + ! check if there are any enclosure polygons + enclosures_present = any(zpl_cellmask(1:Npoly_cellmask) < 0.0_dp) cellmask_initialized = .true. end subroutine dbpinpol_cellmask_init @@ -110,6 +113,7 @@ elemental function dbpinpol_cellmask(xp, yp) result(mask) num_enclosures = 0 count_drypoint = 0 + found_inside_enclosure = .false. ! Single loop over all polygons do ipoly = 1, Npoly_cellmask @@ -127,12 +131,8 @@ elemental function dbpinpol_cellmask(xp, yp) result(mask) if (in_test == 1) then count_drypoint = count_drypoint + 1 end if - else if (zpl_val < 0.0_dp) then - ! Enclosure polygon - num_enclosures = num_enclosures + 1 - if (in_test == 1) then - found_inside_enclosure = .true. - end if + else if (zpl_val < 0.0_dp .and. in_test == 1) then + found_inside_enclosure = .true. end if end do @@ -147,8 +147,9 @@ elemental function dbpinpol_cellmask(xp, yp) result(mask) end if end if - ! Override with enclosure logic if applicable - if (num_enclosures > 0 .and. .not. found_inside_enclosure) then + ! if an enclosure is present, the point must lie inside at least one + ! NOTE: this means we do not handle nested enclosure polygons. + if (enclosures_present .and. .not. found_inside_enclosure) then mask = 1 end if @@ -232,4 +233,4 @@ elemental function pinpok_elemental(xl, yl, ipoly) result(inside) end function pinpok_elemental -end module m_dbpinpol_cellmask \ No newline at end of file +end module m_dbpinpol_cellmask diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/test_pol_to_cellmask.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/test_pol_to_cellmask.f90 index 8216d4b1c0..73f4c6898f 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/test_pol_to_cellmask.f90 +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/test_pol_to_cellmask.f90 @@ -8,7 +8,7 @@ module test_pol_to_cellmask contains - !$f90tw TESTCODE(TEST, test_pol_to_cellmask, test_mixed_polygon, test_nested_drypoint_polygons, + !$f90tw TESTCODE(TEST, test_pol_to_cellmask, test_mixed_polygon, test_mixed_polygon, subroutine test_mixed_polygon() bind(C) ! Test with both enclosure (-1) and dry point (1) polygons integer :: i, k @@ -88,7 +88,9 @@ subroutine test_mixed_polygon() bind(C) deallocate(xzw, yzw, xpl, ypl, zpl, cellmask) end subroutine test_mixed_polygon + !$f90tw) + !$f90tw TESTCODE(TEST, test_pol_to_cellmask, test_nested_drypoint_polygons, test_nested_drypoint_polygons, subroutine test_nested_drypoint_polygons() bind(C) ! Test nested dry point polygons (odd-even rule) integer :: i, k From ad76c0d83cbbc03f3136d77a8c221fbae4a7a5d6 Mon Sep 17 00:00:00 2001 From: FlorisBuwaldaDeltares Date: Wed, 3 Dec 2025 13:30:40 +0100 Subject: [PATCH 14/32] revert temporary out-commenting --- .../dflowfm_kernel/prepost/delete_drypoints_from_netgeom.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/delete_drypoints_from_netgeom.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/delete_drypoints_from_netgeom.f90 index 696210dd7d..f50d4a6a4d 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/delete_drypoints_from_netgeom.f90 +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/delete_drypoints_from_netgeom.f90 @@ -34,7 +34,7 @@ !! Grid enclosures are handled via the jinside=-1 option. module m_delete_drypoints_from_netgeom use m_remove_masked_netcells, only: remove_masked_netcells - !use m_pol_to_cellmask, only: pol_to_cellmask + use m_pol_to_cellmask, only: pol_to_cellmask use m_fix_global_polygons, only: fix_global_polygons implicit none @@ -154,7 +154,7 @@ subroutine delete_drypoints_from_netgeom(dryptsfilelist, jaconfirm, jinside, upd call fix_global_polygons(1, 0) end if - !call pol_to_cellmask() ! third column in pol-file may be used to specify inside (1), or outside (0) mode, only 0 or 1 allowed. + call pol_to_cellmask() ! third column in pol-file may be used to specify inside (1), or outside (0) mode, only 0 or 1 allowed. call delpol() call restorepol() From 8083983de2e840e4dea71c95403aa127725fa5ef Mon Sep 17 00:00:00 2001 From: FlorisBuwaldaDeltares Date: Wed, 3 Dec 2025 13:40:32 +0100 Subject: [PATCH 15/32] forgot to include dmiss in dry area boolean check --- .../src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 index 772bfde968..fb1b580271 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 @@ -91,7 +91,7 @@ subroutine dbpinpol_cellmask_init(NPL, xpl, ypl, zpl) Npoly_cellmask = ipoly ! check if there are any enclosure polygons - enclosures_present = any(zpl_cellmask(1:Npoly_cellmask) < 0.0_dp) + enclosures_present = any(zpl_cellmask(1:Npoly_cellmask) < 0.0_dp .and. zpl_cellmask(1:Npoly_cellmask) /= dmiss) cellmask_initialized = .true. end subroutine dbpinpol_cellmask_init From 012cf87e2c79b1c70e7dd89de31eee121242317f Mon Sep 17 00:00:00 2001 From: FlorisBuwaldaDeltares Date: Wed, 3 Dec 2025 13:56:06 +0100 Subject: [PATCH 16/32] doxygen comments --- .../prepost/dbpinpol_cellmask.f90 | 57 ++++++++++++++----- 1 file changed, 42 insertions(+), 15 deletions(-) diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 index fb1b580271..fb5b82752a 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 @@ -1,3 +1,32 @@ +!----- AGPL -------------------------------------------------------------------- +! +! Copyright (C) Stichting Deltares, 2017-2025. +! +! This file is part of Delft3D (D-Flow Flexible Mesh component). +! +! Delft3D is free software: you can redistribute it and/or modify +! it under the terms of the GNU Affero General Public License as +! published by the Free Software Foundation version 3. +! +! Delft3D is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU Affero General Public License for more details. +! +! You should have received a copy of the GNU Affero General Public License +! along with Delft3D. If not, see . +! +! contact: delft3d.support@deltares.nl +! Stichting Deltares +! P.O. Box 177 +! 2600 MH Delft, The Netherlands +! +! All indications and logos of, and references to, "Delft3D", +! "D-Flow Flexible Mesh" and "Deltares" are registered trademarks of Stichting +! Deltares, and remain the property of Stichting Deltares. All rights reserved. +! +!------------------------------------------------------------------------------- + module m_dbpinpol_cellmask use m_missing, only: jins, dmiss use precision, only: dp @@ -9,25 +38,24 @@ module m_dbpinpol_cellmask !> dbpinpol routines are public to avoid PetSC dependency in unit tests public :: dbpinpol_cellmask_init, dbpinpol_cellmask_cleanup, dbpinpol_cellmask - ! Module-level variables for cellmask polygon checking - real(kind=dp), allocatable :: xpmin_cellmask(:), ypmin_cellmask(:) - real(kind=dp), allocatable :: xpmax_cellmask(:), ypmax_cellmask(:) - real(kind=dp), allocatable :: xpl_cellmask(:), ypl_cellmask(:), zpl_cellmask(:) - integer, allocatable :: iistart_cellmask(:), iiend_cellmask(:) - integer :: Npoly_cellmask = 0 - logical :: cellmask_initialized = .false. - logical :: enclosures_present = .false. + real(kind=dp), allocatable :: xpmin_cellmask(:), ypmin_cellmask(:) !< Polygon bounding box min coordinates + real(kind=dp), allocatable :: xpmax_cellmask(:), ypmax_cellmask(:) !< Polygon bounding box max coordinates + real(kind=dp), allocatable :: xpl_cellmask(:), ypl_cellmask(:), zpl_cellmask(:) !< Polygon coordinate arrays + integer, allocatable :: iistart_cellmask(:), iiend_cellmask(:) !< Polygon start and end indices in coordinate arrays (dim = number of polygons) + integer :: Npoly_cellmask = 0 !< Number of polygons stored in module arrays + logical :: cellmask_initialized = .false. !< Flag indicating if cellmask data structures have been initialized for safety + logical :: enclosures_present = .false. !< Flag indicating if any enclosures are present in the polygon dataset contains + !> Initialize module-level cellmask polygon data structures, such as the bounding boxes and iistart/iiend + ! this keeps the actual calculation routines elemental. subroutine dbpinpol_cellmask_init(NPL, xpl, ypl, zpl) use m_alloc use geometry_module, only: get_startend - implicit none - - integer, intent(in) :: NPL - real(kind=dp), intent(in) :: xpl(NPL), ypl(NPL), zpl(NPL) + integer, intent(in) :: NPL !< Number of polygon points + real(kind=dp), intent(in) :: xpl(NPL), ypl(NPL), zpl(NPL) !< Polygon coordinate arrays integer :: MAXPOLY integer :: ipoint, istart, iend, ipoly @@ -41,7 +69,6 @@ subroutine dbpinpol_cellmask_init(NPL, xpl, ypl, zpl) return end if - ! Allocate and copy full coordinate arrays xpl_cellmask = xpl ypl_cellmask = ypl @@ -96,11 +123,11 @@ subroutine dbpinpol_cellmask_init(NPL, xpl, ypl, zpl) end subroutine dbpinpol_cellmask_init - !> Check if point should be masked - OPTIMIZED VERSION + !> Check if a point should be masked, either inside a dry-area polygon or outside an enclosure polygon. elemental function dbpinpol_cellmask(xp, yp) result(mask) integer :: mask - real(kind=dp), intent(in) :: xp, yp + real(kind=dp), intent(in) :: xp, yp !< Point coordinates integer :: ipoly, in_test integer :: count_drypoint From b815f6e3371591be14b8261092ebe0d9e117dfb5 Mon Sep 17 00:00:00 2001 From: FlorisBuwaldaDeltares Date: Wed, 3 Dec 2025 14:22:49 +0100 Subject: [PATCH 17/32] removed needless copies of m_polygon variables --- .../prepost/dbpinpol_cellmask.f90 | 24 +++++++------------ 1 file changed, 9 insertions(+), 15 deletions(-) diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 index fb5b82752a..bc0cf10c91 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 @@ -30,6 +30,7 @@ module m_dbpinpol_cellmask use m_missing, only: jins, dmiss use precision, only: dp + use m_polygon, only: xpl, ypl, npl implicit none @@ -40,9 +41,8 @@ module m_dbpinpol_cellmask real(kind=dp), allocatable :: xpmin_cellmask(:), ypmin_cellmask(:) !< Polygon bounding box min coordinates real(kind=dp), allocatable :: xpmax_cellmask(:), ypmax_cellmask(:) !< Polygon bounding box max coordinates - real(kind=dp), allocatable :: xpl_cellmask(:), ypl_cellmask(:), zpl_cellmask(:) !< Polygon coordinate arrays + real(kind=dp), allocatable :: zpl_cellmask(:) !< Polygon coordinate arrays integer, allocatable :: iistart_cellmask(:), iiend_cellmask(:) !< Polygon start and end indices in coordinate arrays (dim = number of polygons) - integer :: Npoly_cellmask = 0 !< Number of polygons stored in module arrays logical :: cellmask_initialized = .false. !< Flag indicating if cellmask data structures have been initialized for safety logical :: enclosures_present = .false. !< Flag indicating if any enclosures are present in the polygon dataset @@ -69,9 +69,6 @@ subroutine dbpinpol_cellmask_init(NPL, xpl, ypl, zpl) return end if - xpl_cellmask = xpl - ypl_cellmask = ypl - MAXPOLY = 1000 call realloc(xpmin_cellmask, maxpoly, keepExisting=.false.) call realloc(xpmax_cellmask, maxpoly, keepExisting=.false.) @@ -115,10 +112,8 @@ subroutine dbpinpol_cellmask_init(NPL, xpl, ypl, zpl) ipoint = iend + 2 end do - Npoly_cellmask = ipoly - ! check if there are any enclosure polygons - enclosures_present = any(zpl_cellmask(1:Npoly_cellmask) < 0.0_dp .and. zpl_cellmask(1:Npoly_cellmask) /= dmiss) + enclosures_present = any(zpl_cellmask(1:npl) < 0.0_dp .and. zpl_cellmask(1:npl) /= dmiss) cellmask_initialized = .true. end subroutine dbpinpol_cellmask_init @@ -143,7 +138,7 @@ elemental function dbpinpol_cellmask(xp, yp) result(mask) found_inside_enclosure = .false. ! Single loop over all polygons - do ipoly = 1, Npoly_cellmask + do ipoly = 1, npl zpl_val = zpl_cellmask(ipoly) ! Bounding box check @@ -193,7 +188,6 @@ subroutine dbpinpol_cellmask_cleanup() if (allocated(iistart_cellmask)) deallocate (iistart_cellmask) if (allocated(iiend_cellmask)) deallocate (iiend_cellmask) - Npoly_cellmask = 0 cellmask_initialized = .false. end subroutine dbpinpol_cellmask_cleanup @@ -227,12 +221,12 @@ elemental function pinpok_elemental(xl, yl, ipoly) result(inside) j = iend do i = istart, iend - if (xpl_cellmask(i) == dmiss) exit + if (xpl(i) == dmiss) exit - x1 = xpl_cellmask(j) - y1 = ypl_cellmask(j) - x2 = xpl_cellmask(i) - y2 = ypl_cellmask(i) + x1 = xpl(j) + y1 = ypl(j) + x2 = xpl(i) + y2 = ypl(i) ! Check if point is on vertex if (xl == x1 .and. yl == y1) then From cf60c8adfb2dd50b51411bf9787143a76f071560 Mon Sep 17 00:00:00 2001 From: FlorisBuwaldaDeltares Date: Wed, 3 Dec 2025 16:37:46 +0100 Subject: [PATCH 18/32] restore Npoly_cellmask as it was quite distinct from npl --- .../src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 index bc0cf10c91..2923163a9a 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 @@ -43,6 +43,7 @@ module m_dbpinpol_cellmask real(kind=dp), allocatable :: xpmax_cellmask(:), ypmax_cellmask(:) !< Polygon bounding box max coordinates real(kind=dp), allocatable :: zpl_cellmask(:) !< Polygon coordinate arrays integer, allocatable :: iistart_cellmask(:), iiend_cellmask(:) !< Polygon start and end indices in coordinate arrays (dim = number of polygons) + integer :: Npoly_cellmask = 0 !< Number of polygons stored in module arrays logical :: cellmask_initialized = .false. !< Flag indicating if cellmask data structures have been initialized for safety logical :: enclosures_present = .false. !< Flag indicating if any enclosures are present in the polygon dataset @@ -112,8 +113,10 @@ subroutine dbpinpol_cellmask_init(NPL, xpl, ypl, zpl) ipoint = iend + 2 end do + Npoly_cellmask = ipoly + ! check if there are any enclosure polygons - enclosures_present = any(zpl_cellmask(1:npl) < 0.0_dp .and. zpl_cellmask(1:npl) /= dmiss) + enclosures_present = any(zpl_cellmask(1:Npoly_cellmask) < 0.0_dp .and. zpl_cellmask(1:Npoly_cellmask) /= dmiss) cellmask_initialized = .true. end subroutine dbpinpol_cellmask_init @@ -138,7 +141,7 @@ elemental function dbpinpol_cellmask(xp, yp) result(mask) found_inside_enclosure = .false. ! Single loop over all polygons - do ipoly = 1, npl + do ipoly = 1, Npoly_cellmask zpl_val = zpl_cellmask(ipoly) ! Bounding box check @@ -188,6 +191,7 @@ subroutine dbpinpol_cellmask_cleanup() if (allocated(iistart_cellmask)) deallocate (iistart_cellmask) if (allocated(iiend_cellmask)) deallocate (iiend_cellmask) + Npoly_cellmask = 0 cellmask_initialized = .false. end subroutine dbpinpol_cellmask_cleanup From a476d81ebcc9fd524945cb427f17e282fc67da3c Mon Sep 17 00:00:00 2001 From: FlorisBuwaldaDeltares Date: Thu, 4 Dec 2025 11:54:19 +0100 Subject: [PATCH 19/32] add sequential dependency for unit tests as a workaround --- .../dflowfm/packages/dflowfm_kernel/test/CMakeLists.txt | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/CMakeLists.txt b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/CMakeLists.txt index 767e3fa31e..8d2439b028 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/CMakeLists.txt +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/CMakeLists.txt @@ -24,3 +24,7 @@ f90twtest(test_pol_to_cellmask LIBRARIES dflowfm_kernel VISUAL_STUDIO_PATH engines_gpl/dflowfm/test ) + +# Force sequential build order for the test projects (workaround, to be fixed properly) +add_dependencies(test_3d_layer_administration test_dflowfm_kernel_example) +add_dependencies(test_pol_to_cellmask test_3d_layer_administration) \ No newline at end of file From 417b9acfaaea101f98b363bae55c62643ee86c08 Mon Sep 17 00:00:00 2001 From: FlorisBuwaldaDeltares Date: Mon, 8 Dec 2025 14:18:53 +0100 Subject: [PATCH 20/32] dbpinpol_cellmask => cellmask_from_polygon_set --- ...sk.f90 => m_cellmask_from_polygon_set.f90} | 21 +++++++++---------- .../prepost/pol_to_cellmask.F90 | 8 +++---- .../test/test_pol_to_cellmask.f90 | 14 ++++++------- 3 files changed, 21 insertions(+), 22 deletions(-) rename src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/{dbpinpol_cellmask.f90 => m_cellmask_from_polygon_set.f90} (93%) diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/m_cellmask_from_polygon_set.f90 similarity index 93% rename from src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 rename to src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/m_cellmask_from_polygon_set.f90 index 2923163a9a..2ad317616a 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/m_cellmask_from_polygon_set.f90 @@ -27,7 +27,7 @@ ! !------------------------------------------------------------------------------- -module m_dbpinpol_cellmask +module m_cellmask_from_polygon_set use m_missing, only: jins, dmiss use precision, only: dp use m_polygon, only: xpl, ypl, npl @@ -36,8 +36,7 @@ module m_dbpinpol_cellmask private - !> dbpinpol routines are public to avoid PetSC dependency in unit tests - public :: dbpinpol_cellmask_init, dbpinpol_cellmask_cleanup, dbpinpol_cellmask + public :: cellmask_from_polygon_set_init, cellmask_from_polygon_set_cleanup, cellmask_from_polygon_set real(kind=dp), allocatable :: xpmin_cellmask(:), ypmin_cellmask(:) !< Polygon bounding box min coordinates real(kind=dp), allocatable :: xpmax_cellmask(:), ypmax_cellmask(:) !< Polygon bounding box max coordinates @@ -51,7 +50,7 @@ module m_dbpinpol_cellmask !> Initialize module-level cellmask polygon data structures, such as the bounding boxes and iistart/iiend ! this keeps the actual calculation routines elemental. - subroutine dbpinpol_cellmask_init(NPL, xpl, ypl, zpl) + subroutine cellmask_from_polygon_set_init(NPL, xpl, ypl, zpl) use m_alloc use geometry_module, only: get_startend @@ -62,7 +61,7 @@ subroutine dbpinpol_cellmask_init(NPL, xpl, ypl, zpl) integer :: ipoint, istart, iend, ipoly if (cellmask_initialized) then - call dbpinpol_cellmask_cleanup() + call cellmask_from_polygon_set_cleanup() end if if (NPL == 0) then @@ -119,10 +118,10 @@ subroutine dbpinpol_cellmask_init(NPL, xpl, ypl, zpl) enclosures_present = any(zpl_cellmask(1:Npoly_cellmask) < 0.0_dp .and. zpl_cellmask(1:Npoly_cellmask) /= dmiss) cellmask_initialized = .true. - end subroutine dbpinpol_cellmask_init + end subroutine cellmask_from_polygon_set_init !> Check if a point should be masked, either inside a dry-area polygon or outside an enclosure polygon. - elemental function dbpinpol_cellmask(xp, yp) result(mask) + elemental function cellmask_from_polygon_set(xp, yp) result(mask) integer :: mask real(kind=dp), intent(in) :: xp, yp !< Point coordinates @@ -178,10 +177,10 @@ elemental function dbpinpol_cellmask(xp, yp) result(mask) mask = 1 end if - end function dbpinpol_cellmask + end function cellmask_from_polygon_set !> Clean up module-level cellmask polygon data structures. - subroutine dbpinpol_cellmask_cleanup() + subroutine cellmask_from_polygon_set_cleanup() if (allocated(xpmin_cellmask)) deallocate (xpmin_cellmask) if (allocated(xpmax_cellmask)) deallocate (xpmax_cellmask) @@ -194,7 +193,7 @@ subroutine dbpinpol_cellmask_cleanup() Npoly_cellmask = 0 cellmask_initialized = .false. - end subroutine dbpinpol_cellmask_cleanup + end subroutine cellmask_from_polygon_set_cleanup !> Optimized elemental point-in-polygon test using ray casting algorithm. !! Accesses polygon data via module arrays. @@ -258,4 +257,4 @@ elemental function pinpok_elemental(xl, yl, ipoly) result(inside) end function pinpok_elemental -end module m_dbpinpol_cellmask +end module m_cellmask_from_polygon_set diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.F90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.F90 index ef09524ffc..6359085b36 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.F90 +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.F90 @@ -34,7 +34,7 @@ module m_pol_to_cellmask use precision, only: dp - use m_dbpinpol_cellmask, only: dbpinpol_cellmask_init, dbpinpol_cellmask_cleanup, dbpinpol_cellmask + use m_cellmask_from_polygon_set, only: cellmask_from_polygon_set_init, cellmask_from_polygon_set_cleanup, cellmask_from_polygon_set implicit none private @@ -58,7 +58,7 @@ subroutine pol_to_cellmask() if (NPL == 0) return ! Initialize once - call dbpinpol_cellmask_init(NPL, xpl, ypl, zpl) + call cellmask_from_polygon_set_init(NPL, xpl, ypl, zpl) #ifdef _OPENMP temp_threads = omp_get_max_threads() !> Save old number of threads if (jampi == 0) then @@ -67,7 +67,7 @@ subroutine pol_to_cellmask() #endif !$OMP PARALLEL DO SCHEDULE(DYNAMIC, 100) do k = 1, nump - cellmask(k) = dbpinpol_cellmask(xzw(k), yzw(k)) + cellmask(k) = cellmask_from_polygon_set(xzw(k), yzw(k)) end do !$OMP END PARALLEL DO #ifdef _OPENMP @@ -75,7 +75,7 @@ subroutine pol_to_cellmask() #endif ! Cleanup - call dbpinpol_cellmask_cleanup() + call cellmask_from_polygon_set_cleanup() end subroutine pol_to_cellmask diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/test_pol_to_cellmask.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/test_pol_to_cellmask.f90 index 73f4c6898f..a7eda9ae7e 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/test_pol_to_cellmask.f90 +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/test_pol_to_cellmask.f90 @@ -3,7 +3,7 @@ module test_pol_to_cellmask use precision, only: dp use m_missing, only: dmiss use network_data, only: cellmask, nump1d2d, npl, nump, xzw, yzw, xpl, ypl, zpl - use m_dbpinpol_cellmask, only: dbpinpol_cellmask_init, dbpinpol_cellmask, dbpinpol_cellmask_cleanup + use m_cellmask_from_polygon_set, only: cellmask_from_polygon_set_init, cellmask_from_polygon_set, cellmask_from_polygon_set_cleanup implicit none contains @@ -62,14 +62,14 @@ subroutine test_mixed_polygon() bind(C) xpl(15) = dmiss; ypl(15) = dmiss; zpl(15) = dmiss ! Initialize polygon data structures - call dbpinpol_cellmask_init(NPL, xpl, ypl, zpl) + call cellmask_from_polygon_set_init(NPL, xpl, ypl, zpl) ! Process all cells cellmask = 0 - cellmask = dbpinpol_cellmask(xzw, yzw) + cellmask = cellmask_from_polygon_set(xzw, yzw) ! Cleanup - call dbpinpol_cellmask_cleanup() + call cellmask_from_polygon_set_cleanup() ! Check results: ! Cell at (5,5) - inside enclosure, outside dry point -> mask=0 @@ -143,14 +143,14 @@ subroutine test_nested_drypoint_polygons() bind(C) xpl(13) = dmiss; ypl(13) = dmiss; zpl(13) = dmiss ! Initialize polygon data structures - call dbpinpol_cellmask_init(NPL, xpl, ypl, zpl) + call cellmask_from_polygon_set_init(NPL, xpl, ypl, zpl) ! Process all cells cellmask = 0 - cellmask = dbpinpol_cellmask(xzw, yzw) + cellmask = cellmask_from_polygon_set(xzw, yzw) ! Cleanup - call dbpinpol_cellmask_cleanup() + call cellmask_from_polygon_set_cleanup() ! Check results (odd-even rule): ! Cell at (5,5) - inside 1 polygon (outer) -> mask=1 From 083259331c22d0174ce27e3b8a518133a4428a88 Mon Sep 17 00:00:00 2001 From: FlorisBuwaldaDeltares Date: Mon, 8 Dec 2025 14:37:30 +0100 Subject: [PATCH 21/32] is_inside is now logical, small style fixes --- .../prepost/dbpinpol_cellmask.f90 | 265 ++++++++++++++++++ 1 file changed, 265 insertions(+) create mode 100644 src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 new file mode 100644 index 0000000000..f9b011c913 --- /dev/null +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 @@ -0,0 +1,265 @@ +!----- AGPL -------------------------------------------------------------------- +! +! Copyright (C) Stichting Deltares, 2017-2025. +! +! This file is part of Delft3D (D-Flow Flexible Mesh component). +! +! Delft3D is free software: you can redistribute it and/or modify +! it under the terms of the GNU Affero General Public License as +! published by the Free Software Foundation version 3. +! +! Delft3D is distributed in the hope that it will be useful, +! but WITHOUT ANY WARRANTY; without even the implied warranty of +! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +! GNU Affero General Public License for more details. +! +! You should have received a copy of the GNU Affero General Public License +! along with Delft3D. If not, see . +! +! contact: delft3d.support@deltares.nl +! Stichting Deltares +! P.O. Box 177 +! 2600 MH Delft, The Netherlands +! +! All indications and logos of, and references to, "Delft3D", +! "D-Flow Flexible Mesh" and "Deltares" are registered trademarks of Stichting +! Deltares, and remain the property of Stichting Deltares. All rights reserved. +! +!------------------------------------------------------------------------------- + +module m_cellmask_from_polygon_set + use m_missing, only: jins, dmiss + use precision, only: dp + use m_polygon, only: xpl, ypl, npl + + implicit none + + private + + public :: cellmask_from_polygon_set_init, cellmask_from_polygon_set_cleanup, cellmask_from_polygon_set + + real(kind=dp), allocatable :: xpmin_cellmask(:), ypmin_cellmask(:) !< Polygon bounding box min coordinates + real(kind=dp), allocatable :: xpmax_cellmask(:), ypmax_cellmask(:) !< Polygon bounding box max coordinates + real(kind=dp), allocatable :: zpl_cellmask(:) !< Polygon coordinate arrays + integer, allocatable :: iistart_cellmask(:), iiend_cellmask(:) !< Polygon start and end indices in coordinate arrays (dim = number of polygons) + integer :: Npoly_cellmask = 0 !< Number of polygons stored in module arrays + logical :: cellmask_initialized = .false. !< Flag indicating if cellmask data structures have been initialized for safety + logical :: enclosures_present = .false. !< Flag indicating if any enclosures are present in the polygon dataset + +contains + + !> Initialize module-level cellmask polygon data structures, such as the bounding boxes and iistart/iiend + ! this keeps the actual calculation routines elemental. + subroutine cellmask_from_polygon_set_init(NPL, xpl, ypl, zpl) + use m_alloc + use geometry_module, only: get_startend + + integer, intent(in) :: NPL !< Number of polygon points + real(kind=dp), intent(in) :: xpl(NPL), ypl(NPL), zpl(NPL) !< Polygon coordinate arrays + + integer :: polygon_buffer_size !> polygon arrays buffer size, increases 10% every time to avoid realloc at every polygon + integer :: ipoint, istart, iend, ipoly + + if (cellmask_initialized) then + call cellmask_from_polygon_set_cleanup() + end if + + if (NPL == 0) then + cellmask_initialized = .true. + return + end if + + polygon_buffer_size = 1000 + call realloc(xpmin_cellmask, polygon_buffer_size, keepExisting=.false.) + call realloc(xpmax_cellmask, polygon_buffer_size, keepExisting=.false.) + call realloc(ypmin_cellmask, polygon_buffer_size, keepExisting=.false.) + call realloc(ypmax_cellmask, polygon_buffer_size, keepExisting=.false.) + call realloc(iistart_cellmask, polygon_buffer_size, keepExisting=.false.) + call realloc(iiend_cellmask, polygon_buffer_size, keepExisting=.false.) + call realloc(zpl_cellmask, polygon_buffer_size, keepExisting=.false.) + + ipoint = 1 + ipoly = 0 + + do while (ipoint < NPL) + ipoly = ipoly + 1 + if (ipoly > polygon_buffer_size) then + polygon_buffer_size = ceiling(polygon_buffer_size * 1.1_dp) + call realloc(xpmin_cellmask, polygon_buffer_size, keepExisting=.true.) + call realloc(xpmax_cellmask, polygon_buffer_size, keepExisting=.true.) + call realloc(ypmin_cellmask, polygon_buffer_size, keepExisting=.true.) + call realloc(ypmax_cellmask, polygon_buffer_size, keepExisting=.true.) + call realloc(iistart_cellmask, polygon_buffer_size, keepExisting=.true.) + call realloc(iiend_cellmask, polygon_buffer_size, keepExisting=.true.) + call realloc(zpl_cellmask, polygon_buffer_size, keepExisting=.true.) + end if + + call get_startend(NPL - ipoint + 1, xpl(ipoint:NPL), ypl(ipoint:NPL), istart, iend, dmiss) + istart = istart + ipoint - 1 + iend = iend + ipoint - 1 + + if (istart >= iend .or. iend > NPL) exit + + xpmin_cellmask(ipoly) = minval(xpl(istart:iend)) + xpmax_cellmask(ipoly) = maxval(xpl(istart:iend)) + ypmin_cellmask(ipoly) = minval(ypl(istart:iend)) + ypmax_cellmask(ipoly) = maxval(ypl(istart:iend)) + + iistart_cellmask(ipoly) = istart + iiend_cellmask(ipoly) = iend + zpl_cellmask(ipoly) = zpl(istart) + + ipoint = iend + 2 + end do + + Npoly_cellmask = ipoly + + ! check if there are any enclosure polygons + enclosures_present = any(zpl_cellmask(1:Npoly_cellmask) < 0.0_dp .and. zpl_cellmask(1:Npoly_cellmask) /= dmiss) + cellmask_initialized = .true. + + end subroutine cellmask_from_polygon_set_init + + !> Check if a point should be masked, either is_inside a dry-area polygon or outside an enclosure polygon. + elemental function cellmask_from_polygon_set(xp, yp) result(mask) + + integer :: mask + real(kind=dp), intent(in) :: xp, yp !< Point coordinates + + integer :: count_drypoint, ipoly, num_enclosures + logical :: found_inside_enclosure, is_inside + real(kind=dp) :: zpl_val + + mask = 0 + if (.not. cellmask_initialized) then + return + end if + + num_enclosures = 0 + count_drypoint = 0 + found_inside_enclosure = .false. + is_inside = .false. + + ! Single loop over all polygons + do ipoly = 1, Npoly_cellmask + zpl_val = zpl_cellmask(ipoly) + + ! Bounding box check + if (xp < xpmin_cellmask(ipoly) .or. xp > xpmax_cellmask(ipoly) .or. & + yp < ypmin_cellmask(ipoly) .or. yp > ypmax_cellmask(ipoly)) then + cycle + end if + + ! Point-in-polygon test + is_inside = pinpok_elemental(xp, yp, ipoly) + + if (zpl_val == dmiss .or. zpl_val > 0.0_dp) then + ! Dry point polygon + if (is_inside) then + count_drypoint = count_drypoint + 1 + end if + else if (zpl_val < 0.0_dp .and. is_inside) then + found_inside_enclosure = .true. + end if + end do + + ! Apply odd-even rule only if counting was needed + if (JINS == 1) then + if (mod(count_drypoint, 2) == 1) then + mask = 1 + end if + else + if (mod(count_drypoint, 2) == 0) then + mask = 1 + end if + end if + + ! if an enclosure is present, the point must lie is_inside at least one + ! NOTE: this means we do not handle nested enclosure polygons. + if (enclosures_present .and. .not. found_inside_enclosure) then + mask = 1 + end if + + end function cellmask_from_polygon_set + + !> Clean up module-level cellmask polygon data structures. + subroutine cellmask_from_polygon_set_cleanup() + + if (allocated(xpmin_cellmask)) deallocate (xpmin_cellmask) + if (allocated(xpmax_cellmask)) deallocate (xpmax_cellmask) + if (allocated(ypmin_cellmask)) deallocate (ypmin_cellmask) + if (allocated(ypmax_cellmask)) deallocate (ypmax_cellmask) + if (allocated(zpl_cellmask)) deallocate (zpl_cellmask) + if (allocated(iistart_cellmask)) deallocate (iistart_cellmask) + if (allocated(iiend_cellmask)) deallocate (iiend_cellmask) + + Npoly_cellmask = 0 + cellmask_initialized = .false. + + end subroutine cellmask_from_polygon_set_cleanup + + !> Optimized elemental point-in-polygon test using ray casting algorithm. + !! Accesses polygon data via module arrays. + elemental function pinpok_elemental(xl, yl, ipoly) result(is_inside) + + real(kind=dp), intent(in) :: xl, yl !< Point coordinates (scalar) + integer, intent(in) :: ipoly !< Polygon index + logical :: is_inside !< Result: .true.=is_inside, .false.=outside + + integer :: i, j, istart, iend, crossings + real(kind=dp) :: x1, x2, y1, y2, xinters + + is_inside = .false. + + ! Get polygon bounds from module variables + istart = iistart_cellmask(ipoly) + iend = iiend_cellmask(ipoly) + + if (iend - istart + 1 <= 2) then + is_inside = .true. + goto 999 + end if + + ! Ray-casting algorithm + crossings = 0 + j = iend + + do i = istart, iend + if (xpl(i) == dmiss) then + exit + end if + + x1 = xpl(j) + y1 = ypl(j) + x2 = xpl(i) + y2 = ypl(i) + + ! Check if point is on vertex + if (xl == x1 .and. yl == y1) then + is_inside = .true. + goto 999 + end if + + ! Check if ray crosses edge + if ((y1 > yl) .neqv. (y2 > yl)) then + xinters = x1 + (yl - y1) * (x2 - x1) / (y2 - y1) + + if (xl < xinters) then + crossings = crossings + 1 + else if (xl == xinters) then + is_inside = .true. + goto 999 + end if + end if + j = i + end do + + is_inside = mod(crossings, 2) == 1 +999 continue + if (jins == 0) then + is_inside = .not. is_inside + end if + + end function pinpok_elemental + +end module m_cellmask_from_polygon_set From d0ff67d61421671969b1d66b4da673accbd8dc9d Mon Sep 17 00:00:00 2001 From: FlorisBuwaldaDeltares Date: Mon, 8 Dec 2025 14:37:53 +0100 Subject: [PATCH 22/32] accidental space --- .../src/dflowfm_kernel/prepost/pol_to_cellmask.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.F90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.F90 index 6359085b36..b58e993d8a 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.F90 +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.F90 @@ -1,4 +1,4 @@ - !----- AGPL -------------------------------------------------------------------- +!----- AGPL -------------------------------------------------------------------- ! ! Copyright (C) Stichting Deltares, 2017-2025. ! From 449d699beef802ad5bba430409a61921f088d10d Mon Sep 17 00:00:00 2001 From: FlorisBuwaldaDeltares Date: Mon, 8 Dec 2025 14:51:29 +0100 Subject: [PATCH 23/32] style guide --- .../prepost/dbpinpol_cellmask.f90 | 112 +++++++++--------- 1 file changed, 56 insertions(+), 56 deletions(-) diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 index f9b011c913..43c8d3c179 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 @@ -41,8 +41,8 @@ module m_cellmask_from_polygon_set real(kind=dp), allocatable :: xpmin_cellmask(:), ypmin_cellmask(:) !< Polygon bounding box min coordinates real(kind=dp), allocatable :: xpmax_cellmask(:), ypmax_cellmask(:) !< Polygon bounding box max coordinates real(kind=dp), allocatable :: zpl_cellmask(:) !< Polygon coordinate arrays - integer, allocatable :: iistart_cellmask(:), iiend_cellmask(:) !< Polygon start and end indices in coordinate arrays (dim = number of polygons) - integer :: Npoly_cellmask = 0 !< Number of polygons stored in module arrays + integer, allocatable :: i_start_cellmask(:), i_end_cellmask(:) !< Polygon start and end indices in coordinate arrays (dim = number of polygons) + integer :: polygons = 0 !< Number of polygons stored in module arrays logical :: cellmask_initialized = .false. !< Flag indicating if cellmask data structures have been initialized for safety logical :: enclosures_present = .false. !< Flag indicating if any enclosures are present in the polygon dataset @@ -58,7 +58,7 @@ subroutine cellmask_from_polygon_set_init(NPL, xpl, ypl, zpl) real(kind=dp), intent(in) :: xpl(NPL), ypl(NPL), zpl(NPL) !< Polygon coordinate arrays integer :: polygon_buffer_size !> polygon arrays buffer size, increases 10% every time to avoid realloc at every polygon - integer :: ipoint, istart, iend, ipoly + integer :: i_point, i_start, i_end, i_poly if (cellmask_initialized) then call cellmask_from_polygon_set_cleanup() @@ -74,59 +74,59 @@ subroutine cellmask_from_polygon_set_init(NPL, xpl, ypl, zpl) call realloc(xpmax_cellmask, polygon_buffer_size, keepExisting=.false.) call realloc(ypmin_cellmask, polygon_buffer_size, keepExisting=.false.) call realloc(ypmax_cellmask, polygon_buffer_size, keepExisting=.false.) - call realloc(iistart_cellmask, polygon_buffer_size, keepExisting=.false.) - call realloc(iiend_cellmask, polygon_buffer_size, keepExisting=.false.) + call realloc(i_start_cellmask, polygon_buffer_size, keepExisting=.false.) + call realloc(i_end_cellmask, polygon_buffer_size, keepExisting=.false.) call realloc(zpl_cellmask, polygon_buffer_size, keepExisting=.false.) - ipoint = 1 - ipoly = 0 + i_point = 1 + i_poly = 0 - do while (ipoint < NPL) - ipoly = ipoly + 1 - if (ipoly > polygon_buffer_size) then + do while (i_point < NPL) + i_poly = i_poly + 1 + if (i_poly > polygon_buffer_size) then polygon_buffer_size = ceiling(polygon_buffer_size * 1.1_dp) call realloc(xpmin_cellmask, polygon_buffer_size, keepExisting=.true.) call realloc(xpmax_cellmask, polygon_buffer_size, keepExisting=.true.) call realloc(ypmin_cellmask, polygon_buffer_size, keepExisting=.true.) call realloc(ypmax_cellmask, polygon_buffer_size, keepExisting=.true.) - call realloc(iistart_cellmask, polygon_buffer_size, keepExisting=.true.) - call realloc(iiend_cellmask, polygon_buffer_size, keepExisting=.true.) + call realloc(i_start_cellmask, polygon_buffer_size, keepExisting=.true.) + call realloc(i_end_cellmask, polygon_buffer_size, keepExisting=.true.) call realloc(zpl_cellmask, polygon_buffer_size, keepExisting=.true.) end if - call get_startend(NPL - ipoint + 1, xpl(ipoint:NPL), ypl(ipoint:NPL), istart, iend, dmiss) - istart = istart + ipoint - 1 - iend = iend + ipoint - 1 + call get_startend(NPL - i_point + 1, xpl(i_point:NPL), ypl(i_point:NPL), i_start, i_end, dmiss) + i_start = i_start + i_point - 1 + i_end = i_end + i_point - 1 - if (istart >= iend .or. iend > NPL) exit + if (i_start >= i_end .or. i_end > NPL) exit - xpmin_cellmask(ipoly) = minval(xpl(istart:iend)) - xpmax_cellmask(ipoly) = maxval(xpl(istart:iend)) - ypmin_cellmask(ipoly) = minval(ypl(istart:iend)) - ypmax_cellmask(ipoly) = maxval(ypl(istart:iend)) + xpmin_cellmask(i_poly) = minval(xpl(i_start:i_end)) + xpmax_cellmask(i_poly) = maxval(xpl(i_start:i_end)) + ypmin_cellmask(i_poly) = minval(ypl(i_start:i_end)) + ypmax_cellmask(i_poly) = maxval(ypl(i_start:i_end)) - iistart_cellmask(ipoly) = istart - iiend_cellmask(ipoly) = iend - zpl_cellmask(ipoly) = zpl(istart) + i_start_cellmask(i_poly) = i_start + i_end_cellmask(i_poly) = i_end + zpl_cellmask(i_poly) = zpl(i_start) - ipoint = iend + 2 + i_point = i_end + 2 end do - Npoly_cellmask = ipoly + polygons = i_poly ! check if there are any enclosure polygons - enclosures_present = any(zpl_cellmask(1:Npoly_cellmask) < 0.0_dp .and. zpl_cellmask(1:Npoly_cellmask) /= dmiss) + enclosures_present = any(zpl_cellmask(1:polygons) < 0.0_dp .and. zpl_cellmask(1:polygons) /= dmiss) cellmask_initialized = .true. end subroutine cellmask_from_polygon_set_init !> Check if a point should be masked, either is_inside a dry-area polygon or outside an enclosure polygon. - elemental function cellmask_from_polygon_set(xp, yp) result(mask) + elemental function cellmask_from_polygon_set(x, y) result(mask) integer :: mask - real(kind=dp), intent(in) :: xp, yp !< Point coordinates + real(kind=dp), intent(in) :: x, y !< Point coordinates - integer :: count_drypoint, ipoly, num_enclosures + integer :: count_drypoint, i_poly, num_enclosures logical :: found_inside_enclosure, is_inside real(kind=dp) :: zpl_val @@ -141,17 +141,17 @@ elemental function cellmask_from_polygon_set(xp, yp) result(mask) is_inside = .false. ! Single loop over all polygons - do ipoly = 1, Npoly_cellmask - zpl_val = zpl_cellmask(ipoly) + do i_poly = 1, polygons + zpl_val = zpl_cellmask(i_poly) ! Bounding box check - if (xp < xpmin_cellmask(ipoly) .or. xp > xpmax_cellmask(ipoly) .or. & - yp < ypmin_cellmask(ipoly) .or. yp > ypmax_cellmask(ipoly)) then + if (x < xpmin_cellmask(i_poly) .or. x > xpmax_cellmask(i_poly) .or. & + y < ypmin_cellmask(i_poly) .or. y > ypmax_cellmask(i_poly)) then cycle end if ! Point-in-polygon test - is_inside = pinpok_elemental(xp, yp, ipoly) + is_inside = pinpok_elemental(x, y, i_poly) if (zpl_val == dmiss .or. zpl_val > 0.0_dp) then ! Dry point polygon @@ -190,63 +190,63 @@ subroutine cellmask_from_polygon_set_cleanup() if (allocated(ypmin_cellmask)) deallocate (ypmin_cellmask) if (allocated(ypmax_cellmask)) deallocate (ypmax_cellmask) if (allocated(zpl_cellmask)) deallocate (zpl_cellmask) - if (allocated(iistart_cellmask)) deallocate (iistart_cellmask) - if (allocated(iiend_cellmask)) deallocate (iiend_cellmask) + if (allocated(i_start_cellmask)) deallocate (i_start_cellmask) + if (allocated(i_end_cellmask)) deallocate (i_end_cellmask) - Npoly_cellmask = 0 + polygons = 0 cellmask_initialized = .false. end subroutine cellmask_from_polygon_set_cleanup !> Optimized elemental point-in-polygon test using ray casting algorithm. !! Accesses polygon data via module arrays. - elemental function pinpok_elemental(xl, yl, ipoly) result(is_inside) + elemental function pinpok_elemental(x, y, i_poly) result(is_inside) - real(kind=dp), intent(in) :: xl, yl !< Point coordinates (scalar) - integer, intent(in) :: ipoly !< Polygon index + real(kind=dp), intent(in) :: x, y !< Point coordinates (scalar) + integer, intent(in) :: i_poly !< Polygon index logical :: is_inside !< Result: .true.=is_inside, .false.=outside - integer :: i, j, istart, iend, crossings - real(kind=dp) :: x1, x2, y1, y2, xinters + integer :: i, j, i_start, i_end, crossings + real(kind=dp) :: x_j, x_i, y_j, y_i, x_intersect is_inside = .false. ! Get polygon bounds from module variables - istart = iistart_cellmask(ipoly) - iend = iiend_cellmask(ipoly) + i_start = i_start_cellmask(i_poly) + i_end = i_end_cellmask(i_poly) - if (iend - istart + 1 <= 2) then + if (i_end - i_start + 1 <= 2) then is_inside = .true. goto 999 end if ! Ray-casting algorithm crossings = 0 - j = iend + j = i_end - do i = istart, iend + do i = i_start, i_end if (xpl(i) == dmiss) then exit end if - x1 = xpl(j) - y1 = ypl(j) - x2 = xpl(i) - y2 = ypl(i) + x_j = xpl(j) + y_j = ypl(j) + x_i = xpl(i) + y_i = ypl(i) ! Check if point is on vertex - if (xl == x1 .and. yl == y1) then + if (x == x_j .and. y == y_j) then is_inside = .true. goto 999 end if ! Check if ray crosses edge - if ((y1 > yl) .neqv. (y2 > yl)) then - xinters = x1 + (yl - y1) * (x2 - x1) / (y2 - y1) + if ((y_j > y) .neqv. (y_i > y)) then + x_intersect = x_j + (y - y_j) * (x_i - x_j) / (y_i - y_j) - if (xl < xinters) then + if (x < x_intersect) then crossings = crossings + 1 - else if (xl == xinters) then + else if (x == x_intersect) then is_inside = .true. goto 999 end if From 021dc0400b5cb567a48d783daf0e1f9f0da5e133 Mon Sep 17 00:00:00 2001 From: FlorisBuwaldaDeltares Date: Mon, 8 Dec 2025 14:57:57 +0100 Subject: [PATCH 24/32] style guide renaming --- .../prepost/dbpinpol_cellmask.f90 | 98 +++++++++---------- 1 file changed, 49 insertions(+), 49 deletions(-) diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 index 43c8d3c179..a2a89f7cb5 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 @@ -30,7 +30,6 @@ module m_cellmask_from_polygon_set use m_missing, only: jins, dmiss use precision, only: dp - use m_polygon, only: xpl, ypl, npl implicit none @@ -38,10 +37,10 @@ module m_cellmask_from_polygon_set public :: cellmask_from_polygon_set_init, cellmask_from_polygon_set_cleanup, cellmask_from_polygon_set - real(kind=dp), allocatable :: xpmin_cellmask(:), ypmin_cellmask(:) !< Polygon bounding box min coordinates - real(kind=dp), allocatable :: xpmax_cellmask(:), ypmax_cellmask(:) !< Polygon bounding box max coordinates - real(kind=dp), allocatable :: zpl_cellmask(:) !< Polygon coordinate arrays - integer, allocatable :: i_start_cellmask(:), i_end_cellmask(:) !< Polygon start and end indices in coordinate arrays (dim = number of polygons) + real(kind=dp), allocatable :: x_poly_min(:), y_poly_min(:) !< Polygon bounding box min coordinates + real(kind=dp), allocatable :: x_poly_max(:), y_poly_max(:) !< Polygon bounding box max coordinates + real(kind=dp), allocatable :: z_poly_cellmask(:) !< Polygon coordinate arrays + integer, allocatable :: i_poly_start(:), i_poly_end(:) !< Polygon start and end indices in coordinate arrays (dim = number of polygons) integer :: polygons = 0 !< Number of polygons stored in module arrays logical :: cellmask_initialized = .false. !< Flag indicating if cellmask data structures have been initialized for safety logical :: enclosures_present = .false. !< Flag indicating if any enclosures are present in the polygon dataset @@ -50,12 +49,12 @@ module m_cellmask_from_polygon_set !> Initialize module-level cellmask polygon data structures, such as the bounding boxes and iistart/iiend ! this keeps the actual calculation routines elemental. - subroutine cellmask_from_polygon_set_init(NPL, xpl, ypl, zpl) + subroutine cellmask_from_polygon_set_init(polygon_points, x_poly, y_poly, z_poly) use m_alloc use geometry_module, only: get_startend - integer, intent(in) :: NPL !< Number of polygon points - real(kind=dp), intent(in) :: xpl(NPL), ypl(NPL), zpl(NPL) !< Polygon coordinate arrays + integer, intent(in) :: polygon_points !< Number of polygon points + real(kind=dp), intent(in) :: x_poly(polygon_points), y_poly(polygon_points), z_poly(polygon_points) !< Polygon coordinate arrays integer :: polygon_buffer_size !> polygon arrays buffer size, increases 10% every time to avoid realloc at every polygon integer :: i_point, i_start, i_end, i_poly @@ -64,50 +63,50 @@ subroutine cellmask_from_polygon_set_init(NPL, xpl, ypl, zpl) call cellmask_from_polygon_set_cleanup() end if - if (NPL == 0) then + if (polygon_points == 0) then cellmask_initialized = .true. return end if polygon_buffer_size = 1000 - call realloc(xpmin_cellmask, polygon_buffer_size, keepExisting=.false.) - call realloc(xpmax_cellmask, polygon_buffer_size, keepExisting=.false.) - call realloc(ypmin_cellmask, polygon_buffer_size, keepExisting=.false.) - call realloc(ypmax_cellmask, polygon_buffer_size, keepExisting=.false.) - call realloc(i_start_cellmask, polygon_buffer_size, keepExisting=.false.) - call realloc(i_end_cellmask, polygon_buffer_size, keepExisting=.false.) - call realloc(zpl_cellmask, polygon_buffer_size, keepExisting=.false.) + call realloc(x_poly_min, polygon_buffer_size, keepExisting=.false.) + call realloc(x_poly_max, polygon_buffer_size, keepExisting=.false.) + call realloc(y_poly_min, polygon_buffer_size, keepExisting=.false.) + call realloc(y_poly_max, polygon_buffer_size, keepExisting=.false.) + call realloc(i_poly_start, polygon_buffer_size, keepExisting=.false.) + call realloc(i_poly_end, polygon_buffer_size, keepExisting=.false.) + call realloc(z_poly_cellmask, polygon_buffer_size, keepExisting=.false.) i_point = 1 i_poly = 0 - do while (i_point < NPL) + do while (i_point < polygon_points) i_poly = i_poly + 1 if (i_poly > polygon_buffer_size) then polygon_buffer_size = ceiling(polygon_buffer_size * 1.1_dp) - call realloc(xpmin_cellmask, polygon_buffer_size, keepExisting=.true.) - call realloc(xpmax_cellmask, polygon_buffer_size, keepExisting=.true.) - call realloc(ypmin_cellmask, polygon_buffer_size, keepExisting=.true.) - call realloc(ypmax_cellmask, polygon_buffer_size, keepExisting=.true.) - call realloc(i_start_cellmask, polygon_buffer_size, keepExisting=.true.) - call realloc(i_end_cellmask, polygon_buffer_size, keepExisting=.true.) - call realloc(zpl_cellmask, polygon_buffer_size, keepExisting=.true.) + call realloc(x_poly_min, polygon_buffer_size, keepExisting=.true.) + call realloc(x_poly_max, polygon_buffer_size, keepExisting=.true.) + call realloc(y_poly_min, polygon_buffer_size, keepExisting=.true.) + call realloc(y_poly_max, polygon_buffer_size, keepExisting=.true.) + call realloc(i_poly_start, polygon_buffer_size, keepExisting=.true.) + call realloc(i_poly_end, polygon_buffer_size, keepExisting=.true.) + call realloc(z_poly_cellmask, polygon_buffer_size, keepExisting=.true.) end if - call get_startend(NPL - i_point + 1, xpl(i_point:NPL), ypl(i_point:NPL), i_start, i_end, dmiss) + call get_startend(polygon_points - i_point + 1, x_poly(i_point:polygon_points), y_poly(i_point:polygon_points), i_start, i_end, dmiss) i_start = i_start + i_point - 1 i_end = i_end + i_point - 1 - if (i_start >= i_end .or. i_end > NPL) exit + if (i_start >= i_end .or. i_end > polygon_points) exit - xpmin_cellmask(i_poly) = minval(xpl(i_start:i_end)) - xpmax_cellmask(i_poly) = maxval(xpl(i_start:i_end)) - ypmin_cellmask(i_poly) = minval(ypl(i_start:i_end)) - ypmax_cellmask(i_poly) = maxval(ypl(i_start:i_end)) + x_poly_min(i_poly) = minval(x_poly(i_start:i_end)) + x_poly_max(i_poly) = maxval(x_poly(i_start:i_end)) + y_poly_min(i_poly) = minval(y_poly(i_start:i_end)) + y_poly_max(i_poly) = maxval(y_poly(i_start:i_end)) - i_start_cellmask(i_poly) = i_start - i_end_cellmask(i_poly) = i_end - zpl_cellmask(i_poly) = zpl(i_start) + i_poly_start(i_poly) = i_start + i_poly_end(i_poly) = i_end + z_poly_cellmask(i_poly) = z_poly(i_start) i_point = i_end + 2 end do @@ -115,7 +114,7 @@ subroutine cellmask_from_polygon_set_init(NPL, xpl, ypl, zpl) polygons = i_poly ! check if there are any enclosure polygons - enclosures_present = any(zpl_cellmask(1:polygons) < 0.0_dp .and. zpl_cellmask(1:polygons) /= dmiss) + enclosures_present = any(z_poly_cellmask(1:polygons) < 0.0_dp .and. z_poly_cellmask(1:polygons) /= dmiss) cellmask_initialized = .true. end subroutine cellmask_from_polygon_set_init @@ -128,7 +127,7 @@ elemental function cellmask_from_polygon_set(x, y) result(mask) integer :: count_drypoint, i_poly, num_enclosures logical :: found_inside_enclosure, is_inside - real(kind=dp) :: zpl_val + real(kind=dp) :: z_poly_val mask = 0 if (.not. cellmask_initialized) then @@ -142,23 +141,23 @@ elemental function cellmask_from_polygon_set(x, y) result(mask) ! Single loop over all polygons do i_poly = 1, polygons - zpl_val = zpl_cellmask(i_poly) + z_poly_val = z_poly_cellmask(i_poly) ! Bounding box check - if (x < xpmin_cellmask(i_poly) .or. x > xpmax_cellmask(i_poly) .or. & - y < ypmin_cellmask(i_poly) .or. y > ypmax_cellmask(i_poly)) then + if (x < x_poly_min(i_poly) .or. x > x_poly_max(i_poly) .or. & + y < y_poly_min(i_poly) .or. y > y_poly_max(i_poly)) then cycle end if ! Point-in-polygon test is_inside = pinpok_elemental(x, y, i_poly) - if (zpl_val == dmiss .or. zpl_val > 0.0_dp) then + if (z_poly_val == dmiss .or. z_poly_val > 0.0_dp) then ! Dry point polygon if (is_inside) then count_drypoint = count_drypoint + 1 end if - else if (zpl_val < 0.0_dp .and. is_inside) then + else if (z_poly_val < 0.0_dp .and. is_inside) then found_inside_enclosure = .true. end if end do @@ -185,13 +184,13 @@ end function cellmask_from_polygon_set !> Clean up module-level cellmask polygon data structures. subroutine cellmask_from_polygon_set_cleanup() - if (allocated(xpmin_cellmask)) deallocate (xpmin_cellmask) - if (allocated(xpmax_cellmask)) deallocate (xpmax_cellmask) - if (allocated(ypmin_cellmask)) deallocate (ypmin_cellmask) - if (allocated(ypmax_cellmask)) deallocate (ypmax_cellmask) - if (allocated(zpl_cellmask)) deallocate (zpl_cellmask) - if (allocated(i_start_cellmask)) deallocate (i_start_cellmask) - if (allocated(i_end_cellmask)) deallocate (i_end_cellmask) + if (allocated(x_poly_min)) deallocate (x_poly_min) + if (allocated(x_poly_max)) deallocate (x_poly_max) + if (allocated(y_poly_min)) deallocate (y_poly_min) + if (allocated(y_poly_max)) deallocate (y_poly_max) + if (allocated(z_poly_cellmask)) deallocate (z_poly_cellmask) + if (allocated(i_poly_start)) deallocate (i_poly_start) + if (allocated(i_poly_end)) deallocate (i_poly_end) polygons = 0 cellmask_initialized = .false. @@ -201,6 +200,7 @@ end subroutine cellmask_from_polygon_set_cleanup !> Optimized elemental point-in-polygon test using ray casting algorithm. !! Accesses polygon data via module arrays. elemental function pinpok_elemental(x, y, i_poly) result(is_inside) + use m_polygon, only: xpl, ypl real(kind=dp), intent(in) :: x, y !< Point coordinates (scalar) integer, intent(in) :: i_poly !< Polygon index @@ -212,8 +212,8 @@ elemental function pinpok_elemental(x, y, i_poly) result(is_inside) is_inside = .false. ! Get polygon bounds from module variables - i_start = i_start_cellmask(i_poly) - i_end = i_end_cellmask(i_poly) + i_start = i_poly_start(i_poly) + i_end = i_poly_end(i_poly) if (i_end - i_start + 1 <= 2) then is_inside = .true. From dc90fdfdf4bfe471fa88ecb5265cc202919adb84 Mon Sep 17 00:00:00 2001 From: FlorisBuwaldaDeltares Date: Mon, 8 Dec 2025 14:59:02 +0100 Subject: [PATCH 25/32] remove fancy array statement to prevent stack array allocation --- .../src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 index a2a89f7cb5..9b19612b93 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 @@ -114,7 +114,12 @@ subroutine cellmask_from_polygon_set_init(polygon_points, x_poly, y_poly, z_poly polygons = i_poly ! check if there are any enclosure polygons - enclosures_present = any(z_poly_cellmask(1:polygons) < 0.0_dp .and. z_poly_cellmask(1:polygons) /= dmiss) + do i_poly = 1, polygons + if (z_poly_cellmask(i_poly) < 0.0_dp .and. z_poly_cellmask(i_poly) /= dmiss) then + enclosures_present = .true. + exit + end if + end do cellmask_initialized = .true. end subroutine cellmask_from_polygon_set_init From 71bdea579c92eea048f95b03448f3123f37dd265 Mon Sep 17 00:00:00 2001 From: FlorisBuwaldaDeltares Date: Mon, 8 Dec 2025 15:01:12 +0100 Subject: [PATCH 26/32] remove single line if --- .../src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 index 9b19612b93..4810d968ff 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 @@ -97,7 +97,9 @@ subroutine cellmask_from_polygon_set_init(polygon_points, x_poly, y_poly, z_poly i_start = i_start + i_point - 1 i_end = i_end + i_point - 1 - if (i_start >= i_end .or. i_end > polygon_points) exit + if (i_start >= i_end .or. i_end > polygon_points) then + exit + end if x_poly_min(i_poly) = minval(x_poly(i_start:i_end)) x_poly_max(i_poly) = maxval(x_poly(i_start:i_end)) From 8fe1b77ca734d7190c2cf570f9eb20e29dc049ab Mon Sep 17 00:00:00 2001 From: FlorisBuwaldaDeltares Date: Mon, 8 Dec 2025 15:03:42 +0100 Subject: [PATCH 27/32] rename zpl_cellmask --- .../dflowfm_kernel/prepost/dbpinpol_cellmask.f90 | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 index 4810d968ff..d269d1a534 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 @@ -39,7 +39,7 @@ module m_cellmask_from_polygon_set real(kind=dp), allocatable :: x_poly_min(:), y_poly_min(:) !< Polygon bounding box min coordinates real(kind=dp), allocatable :: x_poly_max(:), y_poly_max(:) !< Polygon bounding box max coordinates - real(kind=dp), allocatable :: z_poly_cellmask(:) !< Polygon coordinate arrays + real(kind=dp), allocatable :: polygon_type(:) !< Polygon type, positive or dmiss = drypoint , negative = enclosure, integer, allocatable :: i_poly_start(:), i_poly_end(:) !< Polygon start and end indices in coordinate arrays (dim = number of polygons) integer :: polygons = 0 !< Number of polygons stored in module arrays logical :: cellmask_initialized = .false. !< Flag indicating if cellmask data structures have been initialized for safety @@ -75,7 +75,7 @@ subroutine cellmask_from_polygon_set_init(polygon_points, x_poly, y_poly, z_poly call realloc(y_poly_max, polygon_buffer_size, keepExisting=.false.) call realloc(i_poly_start, polygon_buffer_size, keepExisting=.false.) call realloc(i_poly_end, polygon_buffer_size, keepExisting=.false.) - call realloc(z_poly_cellmask, polygon_buffer_size, keepExisting=.false.) + call realloc(polygon_type, polygon_buffer_size, keepExisting=.false.) i_point = 1 i_poly = 0 @@ -90,7 +90,7 @@ subroutine cellmask_from_polygon_set_init(polygon_points, x_poly, y_poly, z_poly call realloc(y_poly_max, polygon_buffer_size, keepExisting=.true.) call realloc(i_poly_start, polygon_buffer_size, keepExisting=.true.) call realloc(i_poly_end, polygon_buffer_size, keepExisting=.true.) - call realloc(z_poly_cellmask, polygon_buffer_size, keepExisting=.true.) + call realloc(polygon_type, polygon_buffer_size, keepExisting=.true.) end if call get_startend(polygon_points - i_point + 1, x_poly(i_point:polygon_points), y_poly(i_point:polygon_points), i_start, i_end, dmiss) @@ -108,7 +108,7 @@ subroutine cellmask_from_polygon_set_init(polygon_points, x_poly, y_poly, z_poly i_poly_start(i_poly) = i_start i_poly_end(i_poly) = i_end - z_poly_cellmask(i_poly) = z_poly(i_start) + polygon_type(i_poly) = z_poly(i_start) i_point = i_end + 2 end do @@ -117,7 +117,7 @@ subroutine cellmask_from_polygon_set_init(polygon_points, x_poly, y_poly, z_poly ! check if there are any enclosure polygons do i_poly = 1, polygons - if (z_poly_cellmask(i_poly) < 0.0_dp .and. z_poly_cellmask(i_poly) /= dmiss) then + if (polygon_type(i_poly) < 0.0_dp .and. polygon_type(i_poly) /= dmiss) then enclosures_present = .true. exit end if @@ -148,7 +148,7 @@ elemental function cellmask_from_polygon_set(x, y) result(mask) ! Single loop over all polygons do i_poly = 1, polygons - z_poly_val = z_poly_cellmask(i_poly) + z_poly_val = polygon_type(i_poly) ! Bounding box check if (x < x_poly_min(i_poly) .or. x > x_poly_max(i_poly) .or. & @@ -195,7 +195,7 @@ subroutine cellmask_from_polygon_set_cleanup() if (allocated(x_poly_max)) deallocate (x_poly_max) if (allocated(y_poly_min)) deallocate (y_poly_min) if (allocated(y_poly_max)) deallocate (y_poly_max) - if (allocated(z_poly_cellmask)) deallocate (z_poly_cellmask) + if (allocated(polygon_type)) deallocate (polygon_type) if (allocated(i_poly_start)) deallocate (i_poly_start) if (allocated(i_poly_end)) deallocate (i_poly_end) From 6b9a3da610bdec4ea807a4a4e125dc2c1c51cc18 Mon Sep 17 00:00:00 2001 From: FlorisBuwaldaDeltares Date: Mon, 8 Dec 2025 15:07:28 +0100 Subject: [PATCH 28/32] lowercase jins, doxygen comments --- .../src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 index d269d1a534..abdbefbf64 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 @@ -37,11 +37,11 @@ module m_cellmask_from_polygon_set public :: cellmask_from_polygon_set_init, cellmask_from_polygon_set_cleanup, cellmask_from_polygon_set - real(kind=dp), allocatable :: x_poly_min(:), y_poly_min(:) !< Polygon bounding box min coordinates - real(kind=dp), allocatable :: x_poly_max(:), y_poly_max(:) !< Polygon bounding box max coordinates - real(kind=dp), allocatable :: polygon_type(:) !< Polygon type, positive or dmiss = drypoint , negative = enclosure, - integer, allocatable :: i_poly_start(:), i_poly_end(:) !< Polygon start and end indices in coordinate arrays (dim = number of polygons) - integer :: polygons = 0 !< Number of polygons stored in module arrays + integer :: polygons = 0 !< Number of polygons stored in module arrays xpl, ypl, zpl + real(kind=dp), allocatable :: x_poly_min(:), y_poly_min(:) !< Polygon bounding box min coordinates, (dim = polygons) + real(kind=dp), allocatable :: x_poly_max(:), y_poly_max(:) !< Polygon bounding box max coordinates, (dim = polygons) + real(kind=dp), allocatable :: polygon_type(:) !< Polygon type, positive or dmiss = drypoint , negative = enclosure (dim = polygons) + integer, allocatable :: i_poly_start(:), i_poly_end(:) !< Polygon start and end indices in coordinate arrays (dim = polygons) logical :: cellmask_initialized = .false. !< Flag indicating if cellmask data structures have been initialized for safety logical :: enclosures_present = .false. !< Flag indicating if any enclosures are present in the polygon dataset @@ -170,7 +170,7 @@ elemental function cellmask_from_polygon_set(x, y) result(mask) end do ! Apply odd-even rule only if counting was needed - if (JINS == 1) then + if (jins == 1) then if (mod(count_drypoint, 2) == 1) then mask = 1 end if From 2e53c904479920615e6eba71c005bf7a31825f8e Mon Sep 17 00:00:00 2001 From: FlorisBuwaldaDeltares Date: Mon, 8 Dec 2025 15:13:56 +0100 Subject: [PATCH 29/32] style guide changes in pol_to_cellmask --- .../prepost/pol_to_cellmask.F90 | 21 ++++++++----------- 1 file changed, 9 insertions(+), 12 deletions(-) diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.F90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.F90 index b58e993d8a..ccd474a9e0 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.F90 +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.F90 @@ -1,4 +1,4 @@ -!----- AGPL -------------------------------------------------------------------- + !----- AGPL -------------------------------------------------------------------- ! ! Copyright (C) Stichting Deltares, 2017-2025. ! @@ -27,14 +27,11 @@ ! !------------------------------------------------------------------------------- -!> update cellmask from samples -!> a cell is dry when it is: -!> 1) inside ANY "1"-polygon (drypnts), OR -!> 2) outside ALL "-1"-polygons (enclosures) +!> Wrapper around cellmask_from_polygon_set that uses OpenMP to parallelize the loop over all points if not in MPI mode module m_pol_to_cellmask - use precision, only: dp use m_cellmask_from_polygon_set, only: cellmask_from_polygon_set_init, cellmask_from_polygon_set_cleanup, cellmask_from_polygon_set + implicit none private @@ -46,18 +43,19 @@ module m_pol_to_cellmask subroutine pol_to_cellmask() use network_data, only: cellmask, nump1d2d, npl, nump, xzw, yzw, xpl, ypl, zpl use m_partitioninfo, only: jampi + use m_alloc, only: realloc #ifdef _OPENMP use omp_lib integer :: temp_threads #endif integer :: k - if (allocated(cellmask)) deallocate (cellmask) - allocate (cellmask(nump1d2d)) - cellmask = 0 - if (NPL == 0) return + if (NPL == 0) then + return + end if + + call realloc(cellmask, nump1d2d, keepexisting=.false., fill=0) - ! Initialize once call cellmask_from_polygon_set_init(NPL, xpl, ypl, zpl) #ifdef _OPENMP temp_threads = omp_get_max_threads() !> Save old number of threads @@ -74,7 +72,6 @@ subroutine pol_to_cellmask() call omp_set_num_threads(temp_threads) #endif - ! Cleanup call cellmask_from_polygon_set_cleanup() end subroutine pol_to_cellmask From da9f59e33e85da748c06482109d0dcfd9befce9b Mon Sep 17 00:00:00 2001 From: FlorisBuwaldaDeltares Date: Mon, 8 Dec 2025 15:40:49 +0100 Subject: [PATCH 30/32] retry rename of dbpinpol_cellmask --- .../prepost/dbpinpol_cellmask.f90 | 272 ------------------ .../prepost/m_cellmask_from_polygon_set.f90 | 222 +++++++------- 2 files changed, 117 insertions(+), 377 deletions(-) delete mode 100644 src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 deleted file mode 100644 index abdbefbf64..0000000000 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90 +++ /dev/null @@ -1,272 +0,0 @@ -!----- AGPL -------------------------------------------------------------------- -! -! Copyright (C) Stichting Deltares, 2017-2025. -! -! This file is part of Delft3D (D-Flow Flexible Mesh component). -! -! Delft3D is free software: you can redistribute it and/or modify -! it under the terms of the GNU Affero General Public License as -! published by the Free Software Foundation version 3. -! -! Delft3D is distributed in the hope that it will be useful, -! but WITHOUT ANY WARRANTY; without even the implied warranty of -! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the -! GNU Affero General Public License for more details. -! -! You should have received a copy of the GNU Affero General Public License -! along with Delft3D. If not, see . -! -! contact: delft3d.support@deltares.nl -! Stichting Deltares -! P.O. Box 177 -! 2600 MH Delft, The Netherlands -! -! All indications and logos of, and references to, "Delft3D", -! "D-Flow Flexible Mesh" and "Deltares" are registered trademarks of Stichting -! Deltares, and remain the property of Stichting Deltares. All rights reserved. -! -!------------------------------------------------------------------------------- - -module m_cellmask_from_polygon_set - use m_missing, only: jins, dmiss - use precision, only: dp - - implicit none - - private - - public :: cellmask_from_polygon_set_init, cellmask_from_polygon_set_cleanup, cellmask_from_polygon_set - - integer :: polygons = 0 !< Number of polygons stored in module arrays xpl, ypl, zpl - real(kind=dp), allocatable :: x_poly_min(:), y_poly_min(:) !< Polygon bounding box min coordinates, (dim = polygons) - real(kind=dp), allocatable :: x_poly_max(:), y_poly_max(:) !< Polygon bounding box max coordinates, (dim = polygons) - real(kind=dp), allocatable :: polygon_type(:) !< Polygon type, positive or dmiss = drypoint , negative = enclosure (dim = polygons) - integer, allocatable :: i_poly_start(:), i_poly_end(:) !< Polygon start and end indices in coordinate arrays (dim = polygons) - logical :: cellmask_initialized = .false. !< Flag indicating if cellmask data structures have been initialized for safety - logical :: enclosures_present = .false. !< Flag indicating if any enclosures are present in the polygon dataset - -contains - - !> Initialize module-level cellmask polygon data structures, such as the bounding boxes and iistart/iiend - ! this keeps the actual calculation routines elemental. - subroutine cellmask_from_polygon_set_init(polygon_points, x_poly, y_poly, z_poly) - use m_alloc - use geometry_module, only: get_startend - - integer, intent(in) :: polygon_points !< Number of polygon points - real(kind=dp), intent(in) :: x_poly(polygon_points), y_poly(polygon_points), z_poly(polygon_points) !< Polygon coordinate arrays - - integer :: polygon_buffer_size !> polygon arrays buffer size, increases 10% every time to avoid realloc at every polygon - integer :: i_point, i_start, i_end, i_poly - - if (cellmask_initialized) then - call cellmask_from_polygon_set_cleanup() - end if - - if (polygon_points == 0) then - cellmask_initialized = .true. - return - end if - - polygon_buffer_size = 1000 - call realloc(x_poly_min, polygon_buffer_size, keepExisting=.false.) - call realloc(x_poly_max, polygon_buffer_size, keepExisting=.false.) - call realloc(y_poly_min, polygon_buffer_size, keepExisting=.false.) - call realloc(y_poly_max, polygon_buffer_size, keepExisting=.false.) - call realloc(i_poly_start, polygon_buffer_size, keepExisting=.false.) - call realloc(i_poly_end, polygon_buffer_size, keepExisting=.false.) - call realloc(polygon_type, polygon_buffer_size, keepExisting=.false.) - - i_point = 1 - i_poly = 0 - - do while (i_point < polygon_points) - i_poly = i_poly + 1 - if (i_poly > polygon_buffer_size) then - polygon_buffer_size = ceiling(polygon_buffer_size * 1.1_dp) - call realloc(x_poly_min, polygon_buffer_size, keepExisting=.true.) - call realloc(x_poly_max, polygon_buffer_size, keepExisting=.true.) - call realloc(y_poly_min, polygon_buffer_size, keepExisting=.true.) - call realloc(y_poly_max, polygon_buffer_size, keepExisting=.true.) - call realloc(i_poly_start, polygon_buffer_size, keepExisting=.true.) - call realloc(i_poly_end, polygon_buffer_size, keepExisting=.true.) - call realloc(polygon_type, polygon_buffer_size, keepExisting=.true.) - end if - - call get_startend(polygon_points - i_point + 1, x_poly(i_point:polygon_points), y_poly(i_point:polygon_points), i_start, i_end, dmiss) - i_start = i_start + i_point - 1 - i_end = i_end + i_point - 1 - - if (i_start >= i_end .or. i_end > polygon_points) then - exit - end if - - x_poly_min(i_poly) = minval(x_poly(i_start:i_end)) - x_poly_max(i_poly) = maxval(x_poly(i_start:i_end)) - y_poly_min(i_poly) = minval(y_poly(i_start:i_end)) - y_poly_max(i_poly) = maxval(y_poly(i_start:i_end)) - - i_poly_start(i_poly) = i_start - i_poly_end(i_poly) = i_end - polygon_type(i_poly) = z_poly(i_start) - - i_point = i_end + 2 - end do - - polygons = i_poly - - ! check if there are any enclosure polygons - do i_poly = 1, polygons - if (polygon_type(i_poly) < 0.0_dp .and. polygon_type(i_poly) /= dmiss) then - enclosures_present = .true. - exit - end if - end do - cellmask_initialized = .true. - - end subroutine cellmask_from_polygon_set_init - - !> Check if a point should be masked, either is_inside a dry-area polygon or outside an enclosure polygon. - elemental function cellmask_from_polygon_set(x, y) result(mask) - - integer :: mask - real(kind=dp), intent(in) :: x, y !< Point coordinates - - integer :: count_drypoint, i_poly, num_enclosures - logical :: found_inside_enclosure, is_inside - real(kind=dp) :: z_poly_val - - mask = 0 - if (.not. cellmask_initialized) then - return - end if - - num_enclosures = 0 - count_drypoint = 0 - found_inside_enclosure = .false. - is_inside = .false. - - ! Single loop over all polygons - do i_poly = 1, polygons - z_poly_val = polygon_type(i_poly) - - ! Bounding box check - if (x < x_poly_min(i_poly) .or. x > x_poly_max(i_poly) .or. & - y < y_poly_min(i_poly) .or. y > y_poly_max(i_poly)) then - cycle - end if - - ! Point-in-polygon test - is_inside = pinpok_elemental(x, y, i_poly) - - if (z_poly_val == dmiss .or. z_poly_val > 0.0_dp) then - ! Dry point polygon - if (is_inside) then - count_drypoint = count_drypoint + 1 - end if - else if (z_poly_val < 0.0_dp .and. is_inside) then - found_inside_enclosure = .true. - end if - end do - - ! Apply odd-even rule only if counting was needed - if (jins == 1) then - if (mod(count_drypoint, 2) == 1) then - mask = 1 - end if - else - if (mod(count_drypoint, 2) == 0) then - mask = 1 - end if - end if - - ! if an enclosure is present, the point must lie is_inside at least one - ! NOTE: this means we do not handle nested enclosure polygons. - if (enclosures_present .and. .not. found_inside_enclosure) then - mask = 1 - end if - - end function cellmask_from_polygon_set - - !> Clean up module-level cellmask polygon data structures. - subroutine cellmask_from_polygon_set_cleanup() - - if (allocated(x_poly_min)) deallocate (x_poly_min) - if (allocated(x_poly_max)) deallocate (x_poly_max) - if (allocated(y_poly_min)) deallocate (y_poly_min) - if (allocated(y_poly_max)) deallocate (y_poly_max) - if (allocated(polygon_type)) deallocate (polygon_type) - if (allocated(i_poly_start)) deallocate (i_poly_start) - if (allocated(i_poly_end)) deallocate (i_poly_end) - - polygons = 0 - cellmask_initialized = .false. - - end subroutine cellmask_from_polygon_set_cleanup - - !> Optimized elemental point-in-polygon test using ray casting algorithm. - !! Accesses polygon data via module arrays. - elemental function pinpok_elemental(x, y, i_poly) result(is_inside) - use m_polygon, only: xpl, ypl - - real(kind=dp), intent(in) :: x, y !< Point coordinates (scalar) - integer, intent(in) :: i_poly !< Polygon index - logical :: is_inside !< Result: .true.=is_inside, .false.=outside - - integer :: i, j, i_start, i_end, crossings - real(kind=dp) :: x_j, x_i, y_j, y_i, x_intersect - - is_inside = .false. - - ! Get polygon bounds from module variables - i_start = i_poly_start(i_poly) - i_end = i_poly_end(i_poly) - - if (i_end - i_start + 1 <= 2) then - is_inside = .true. - goto 999 - end if - - ! Ray-casting algorithm - crossings = 0 - j = i_end - - do i = i_start, i_end - if (xpl(i) == dmiss) then - exit - end if - - x_j = xpl(j) - y_j = ypl(j) - x_i = xpl(i) - y_i = ypl(i) - - ! Check if point is on vertex - if (x == x_j .and. y == y_j) then - is_inside = .true. - goto 999 - end if - - ! Check if ray crosses edge - if ((y_j > y) .neqv. (y_i > y)) then - x_intersect = x_j + (y - y_j) * (x_i - x_j) / (y_i - y_j) - - if (x < x_intersect) then - crossings = crossings + 1 - else if (x == x_intersect) then - is_inside = .true. - goto 999 - end if - end if - j = i - end do - - is_inside = mod(crossings, 2) == 1 -999 continue - if (jins == 0) then - is_inside = .not. is_inside - end if - - end function pinpok_elemental - -end module m_cellmask_from_polygon_set diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/m_cellmask_from_polygon_set.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/m_cellmask_from_polygon_set.f90 index 2ad317616a..abdbefbf64 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/m_cellmask_from_polygon_set.f90 +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/m_cellmask_from_polygon_set.f90 @@ -30,7 +30,6 @@ module m_cellmask_from_polygon_set use m_missing, only: jins, dmiss use precision, only: dp - use m_polygon, only: xpl, ypl, npl implicit none @@ -38,11 +37,11 @@ module m_cellmask_from_polygon_set public :: cellmask_from_polygon_set_init, cellmask_from_polygon_set_cleanup, cellmask_from_polygon_set - real(kind=dp), allocatable :: xpmin_cellmask(:), ypmin_cellmask(:) !< Polygon bounding box min coordinates - real(kind=dp), allocatable :: xpmax_cellmask(:), ypmax_cellmask(:) !< Polygon bounding box max coordinates - real(kind=dp), allocatable :: zpl_cellmask(:) !< Polygon coordinate arrays - integer, allocatable :: iistart_cellmask(:), iiend_cellmask(:) !< Polygon start and end indices in coordinate arrays (dim = number of polygons) - integer :: Npoly_cellmask = 0 !< Number of polygons stored in module arrays + integer :: polygons = 0 !< Number of polygons stored in module arrays xpl, ypl, zpl + real(kind=dp), allocatable :: x_poly_min(:), y_poly_min(:) !< Polygon bounding box min coordinates, (dim = polygons) + real(kind=dp), allocatable :: x_poly_max(:), y_poly_max(:) !< Polygon bounding box max coordinates, (dim = polygons) + real(kind=dp), allocatable :: polygon_type(:) !< Polygon type, positive or dmiss = drypoint , negative = enclosure (dim = polygons) + integer, allocatable :: i_poly_start(:), i_poly_end(:) !< Polygon start and end indices in coordinate arrays (dim = polygons) logical :: cellmask_initialized = .false. !< Flag indicating if cellmask data structures have been initialized for safety logical :: enclosures_present = .false. !< Flag indicating if any enclosures are present in the polygon dataset @@ -50,118 +49,128 @@ module m_cellmask_from_polygon_set !> Initialize module-level cellmask polygon data structures, such as the bounding boxes and iistart/iiend ! this keeps the actual calculation routines elemental. - subroutine cellmask_from_polygon_set_init(NPL, xpl, ypl, zpl) + subroutine cellmask_from_polygon_set_init(polygon_points, x_poly, y_poly, z_poly) use m_alloc use geometry_module, only: get_startend - integer, intent(in) :: NPL !< Number of polygon points - real(kind=dp), intent(in) :: xpl(NPL), ypl(NPL), zpl(NPL) !< Polygon coordinate arrays + integer, intent(in) :: polygon_points !< Number of polygon points + real(kind=dp), intent(in) :: x_poly(polygon_points), y_poly(polygon_points), z_poly(polygon_points) !< Polygon coordinate arrays - integer :: MAXPOLY - integer :: ipoint, istart, iend, ipoly + integer :: polygon_buffer_size !> polygon arrays buffer size, increases 10% every time to avoid realloc at every polygon + integer :: i_point, i_start, i_end, i_poly if (cellmask_initialized) then call cellmask_from_polygon_set_cleanup() end if - if (NPL == 0) then + if (polygon_points == 0) then cellmask_initialized = .true. return end if - MAXPOLY = 1000 - call realloc(xpmin_cellmask, maxpoly, keepExisting=.false.) - call realloc(xpmax_cellmask, maxpoly, keepExisting=.false.) - call realloc(ypmin_cellmask, maxpoly, keepExisting=.false.) - call realloc(ypmax_cellmask, maxpoly, keepExisting=.false.) - call realloc(iistart_cellmask, maxpoly, keepExisting=.false.) - call realloc(iiend_cellmask, maxpoly, keepExisting=.false.) - call realloc(zpl_cellmask, maxpoly, keepExisting=.false.) - - ipoint = 1 - ipoly = 0 - - do while (ipoint < NPL) - ipoly = ipoly + 1 - if (ipoly > maxpoly) then - maxpoly = ceiling(maxpoly * 1.1) - call realloc(xpmin_cellmask, maxpoly, keepExisting=.true.) - call realloc(xpmax_cellmask, maxpoly, keepExisting=.true.) - call realloc(ypmin_cellmask, maxpoly, keepExisting=.true.) - call realloc(ypmax_cellmask, maxpoly, keepExisting=.true.) - call realloc(iistart_cellmask, maxpoly, keepExisting=.true.) - call realloc(iiend_cellmask, maxpoly, keepExisting=.true.) - call realloc(zpl_cellmask, maxpoly, keepExisting=.true.) + polygon_buffer_size = 1000 + call realloc(x_poly_min, polygon_buffer_size, keepExisting=.false.) + call realloc(x_poly_max, polygon_buffer_size, keepExisting=.false.) + call realloc(y_poly_min, polygon_buffer_size, keepExisting=.false.) + call realloc(y_poly_max, polygon_buffer_size, keepExisting=.false.) + call realloc(i_poly_start, polygon_buffer_size, keepExisting=.false.) + call realloc(i_poly_end, polygon_buffer_size, keepExisting=.false.) + call realloc(polygon_type, polygon_buffer_size, keepExisting=.false.) + + i_point = 1 + i_poly = 0 + + do while (i_point < polygon_points) + i_poly = i_poly + 1 + if (i_poly > polygon_buffer_size) then + polygon_buffer_size = ceiling(polygon_buffer_size * 1.1_dp) + call realloc(x_poly_min, polygon_buffer_size, keepExisting=.true.) + call realloc(x_poly_max, polygon_buffer_size, keepExisting=.true.) + call realloc(y_poly_min, polygon_buffer_size, keepExisting=.true.) + call realloc(y_poly_max, polygon_buffer_size, keepExisting=.true.) + call realloc(i_poly_start, polygon_buffer_size, keepExisting=.true.) + call realloc(i_poly_end, polygon_buffer_size, keepExisting=.true.) + call realloc(polygon_type, polygon_buffer_size, keepExisting=.true.) end if - call get_startend(NPL - ipoint + 1, xpl(ipoint:NPL), ypl(ipoint:NPL), istart, iend, dmiss) - istart = istart + ipoint - 1 - iend = iend + ipoint - 1 + call get_startend(polygon_points - i_point + 1, x_poly(i_point:polygon_points), y_poly(i_point:polygon_points), i_start, i_end, dmiss) + i_start = i_start + i_point - 1 + i_end = i_end + i_point - 1 - if (istart >= iend .or. iend > NPL) exit + if (i_start >= i_end .or. i_end > polygon_points) then + exit + end if - xpmin_cellmask(ipoly) = minval(xpl(istart:iend)) - xpmax_cellmask(ipoly) = maxval(xpl(istart:iend)) - ypmin_cellmask(ipoly) = minval(ypl(istart:iend)) - ypmax_cellmask(ipoly) = maxval(ypl(istart:iend)) + x_poly_min(i_poly) = minval(x_poly(i_start:i_end)) + x_poly_max(i_poly) = maxval(x_poly(i_start:i_end)) + y_poly_min(i_poly) = minval(y_poly(i_start:i_end)) + y_poly_max(i_poly) = maxval(y_poly(i_start:i_end)) - iistart_cellmask(ipoly) = istart - iiend_cellmask(ipoly) = iend - zpl_cellmask(ipoly) = zpl(istart) + i_poly_start(i_poly) = i_start + i_poly_end(i_poly) = i_end + polygon_type(i_poly) = z_poly(i_start) - ipoint = iend + 2 + i_point = i_end + 2 end do - Npoly_cellmask = ipoly + polygons = i_poly ! check if there are any enclosure polygons - enclosures_present = any(zpl_cellmask(1:Npoly_cellmask) < 0.0_dp .and. zpl_cellmask(1:Npoly_cellmask) /= dmiss) + do i_poly = 1, polygons + if (polygon_type(i_poly) < 0.0_dp .and. polygon_type(i_poly) /= dmiss) then + enclosures_present = .true. + exit + end if + end do cellmask_initialized = .true. end subroutine cellmask_from_polygon_set_init - !> Check if a point should be masked, either inside a dry-area polygon or outside an enclosure polygon. - elemental function cellmask_from_polygon_set(xp, yp) result(mask) + !> Check if a point should be masked, either is_inside a dry-area polygon or outside an enclosure polygon. + elemental function cellmask_from_polygon_set(x, y) result(mask) integer :: mask - real(kind=dp), intent(in) :: xp, yp !< Point coordinates + real(kind=dp), intent(in) :: x, y !< Point coordinates - integer :: ipoly, in_test - integer :: count_drypoint - logical :: found_inside_enclosure - integer :: num_enclosures - real(kind=dp) :: zpl_val + integer :: count_drypoint, i_poly, num_enclosures + logical :: found_inside_enclosure, is_inside + real(kind=dp) :: z_poly_val mask = 0 - if (.not. cellmask_initialized) return + if (.not. cellmask_initialized) then + return + end if num_enclosures = 0 count_drypoint = 0 found_inside_enclosure = .false. + is_inside = .false. ! Single loop over all polygons - do ipoly = 1, Npoly_cellmask - zpl_val = zpl_cellmask(ipoly) + do i_poly = 1, polygons + z_poly_val = polygon_type(i_poly) ! Bounding box check - if (xp < xpmin_cellmask(ipoly) .or. xp > xpmax_cellmask(ipoly) .or. & - yp < ypmin_cellmask(ipoly) .or. yp > ypmax_cellmask(ipoly)) cycle + if (x < x_poly_min(i_poly) .or. x > x_poly_max(i_poly) .or. & + y < y_poly_min(i_poly) .or. y > y_poly_max(i_poly)) then + cycle + end if ! Point-in-polygon test - in_test = pinpok_elemental(xp, yp, ipoly) + is_inside = pinpok_elemental(x, y, i_poly) - if (zpl_val == dmiss .or. zpl_val > 0.0_dp) then + if (z_poly_val == dmiss .or. z_poly_val > 0.0_dp) then ! Dry point polygon - if (in_test == 1) then + if (is_inside) then count_drypoint = count_drypoint + 1 end if - else if (zpl_val < 0.0_dp .and. in_test == 1) then + else if (z_poly_val < 0.0_dp .and. is_inside) then found_inside_enclosure = .true. end if end do ! Apply odd-even rule only if counting was needed - if (JINS == 1) then + if (jins == 1) then if (mod(count_drypoint, 2) == 1) then mask = 1 end if @@ -171,7 +180,7 @@ elemental function cellmask_from_polygon_set(xp, yp) result(mask) end if end if - ! if an enclosure is present, the point must lie inside at least one + ! if an enclosure is present, the point must lie is_inside at least one ! NOTE: this means we do not handle nested enclosure polygons. if (enclosures_present .and. .not. found_inside_enclosure) then mask = 1 @@ -182,78 +191,81 @@ end function cellmask_from_polygon_set !> Clean up module-level cellmask polygon data structures. subroutine cellmask_from_polygon_set_cleanup() - if (allocated(xpmin_cellmask)) deallocate (xpmin_cellmask) - if (allocated(xpmax_cellmask)) deallocate (xpmax_cellmask) - if (allocated(ypmin_cellmask)) deallocate (ypmin_cellmask) - if (allocated(ypmax_cellmask)) deallocate (ypmax_cellmask) - if (allocated(zpl_cellmask)) deallocate (zpl_cellmask) - if (allocated(iistart_cellmask)) deallocate (iistart_cellmask) - if (allocated(iiend_cellmask)) deallocate (iiend_cellmask) + if (allocated(x_poly_min)) deallocate (x_poly_min) + if (allocated(x_poly_max)) deallocate (x_poly_max) + if (allocated(y_poly_min)) deallocate (y_poly_min) + if (allocated(y_poly_max)) deallocate (y_poly_max) + if (allocated(polygon_type)) deallocate (polygon_type) + if (allocated(i_poly_start)) deallocate (i_poly_start) + if (allocated(i_poly_end)) deallocate (i_poly_end) - Npoly_cellmask = 0 + polygons = 0 cellmask_initialized = .false. end subroutine cellmask_from_polygon_set_cleanup !> Optimized elemental point-in-polygon test using ray casting algorithm. !! Accesses polygon data via module arrays. - elemental function pinpok_elemental(xl, yl, ipoly) result(inside) - - implicit none + elemental function pinpok_elemental(x, y, i_poly) result(is_inside) + use m_polygon, only: xpl, ypl - real(kind=dp), intent(in) :: xl, yl !< Point coordinates (scalar) - integer, intent(in) :: ipoly !< Polygon index - integer :: inside !< Result: 1=inside, 0=outside + real(kind=dp), intent(in) :: x, y !< Point coordinates (scalar) + integer, intent(in) :: i_poly !< Polygon index + logical :: is_inside !< Result: .true.=is_inside, .false.=outside - integer :: i, j, istart, iend, crossings - real(kind=dp) :: x1, x2, y1, y2, xinters + integer :: i, j, i_start, i_end, crossings + real(kind=dp) :: x_j, x_i, y_j, y_i, x_intersect - inside = 0 + is_inside = .false. ! Get polygon bounds from module variables - istart = iistart_cellmask(ipoly) - iend = iiend_cellmask(ipoly) + i_start = i_poly_start(i_poly) + i_end = i_poly_end(i_poly) - if (iend - istart + 1 <= 2) then - inside = 1 + if (i_end - i_start + 1 <= 2) then + is_inside = .true. goto 999 end if ! Ray-casting algorithm crossings = 0 - j = iend + j = i_end - do i = istart, iend - if (xpl(i) == dmiss) exit + do i = i_start, i_end + if (xpl(i) == dmiss) then + exit + end if - x1 = xpl(j) - y1 = ypl(j) - x2 = xpl(i) - y2 = ypl(i) + x_j = xpl(j) + y_j = ypl(j) + x_i = xpl(i) + y_i = ypl(i) ! Check if point is on vertex - if (xl == x1 .and. yl == y1) then - inside = 1 + if (x == x_j .and. y == y_j) then + is_inside = .true. goto 999 end if ! Check if ray crosses edge - if ((y1 > yl) .neqv. (y2 > yl)) then - xinters = x1 + (yl - y1) * (x2 - x1) / (y2 - y1) + if ((y_j > y) .neqv. (y_i > y)) then + x_intersect = x_j + (y - y_j) * (x_i - x_j) / (y_i - y_j) - if (xl < xinters) then + if (x < x_intersect) then crossings = crossings + 1 - else if (xl == xinters) then - inside = 1 + else if (x == x_intersect) then + is_inside = .true. goto 999 end if end if j = i end do - inside = mod(crossings, 2) + is_inside = mod(crossings, 2) == 1 999 continue - if (jins == 0) inside = 1 - inside + if (jins == 0) then + is_inside = .not. is_inside + end if end function pinpok_elemental From 02bfdb9eea1dfd41a3ffcf623581ade163714edd Mon Sep 17 00:00:00 2001 From: FlorisBuwaldaDeltares Date: Mon, 8 Dec 2025 15:58:13 +0100 Subject: [PATCH 31/32] no separators --- .../test/test_pol_to_cellmask.f90 | 112 +++++++++++++----- 1 file changed, 84 insertions(+), 28 deletions(-) diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/test_pol_to_cellmask.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/test_pol_to_cellmask.f90 index a7eda9ae7e..eb735c0405 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/test_pol_to_cellmask.f90 +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/test_pol_to_cellmask.f90 @@ -39,27 +39,57 @@ subroutine test_mixed_polygon() bind(C) allocate(xpl(npl), ypl(npl), zpl(npl)) ! Enclosure polygon (rectangle from 0,0 to 30,30) with zpl=-1 - xpl(1) = 0.0_dp; ypl(1) = 0.0_dp; zpl(1) = -1.0_dp - xpl(2) = 30.0_dp; ypl(2) = 0.0_dp; zpl(2) = -1.0_dp - xpl(3) = 30.0_dp; ypl(3) = 30.0_dp; zpl(3) = -1.0_dp - xpl(4) = 0.0_dp; ypl(4) = 30.0_dp; zpl(4) = -1.0_dp - xpl(5) = 0.0_dp; ypl(5) = 0.0_dp; zpl(5) = -1.0_dp + xpl(1) = 0.0_dp + ypl(1) = 0.0_dp + zpl(1) = -1.0_dp + xpl(2) = 30.0_dp + ypl(2) = 0.0_dp + zpl(2) = -1.0_dp + xpl(3) = 30.0_dp + ypl(3) = 30.0_dp + zpl(3) = -1.0_dp + xpl(4) = 0.0_dp + ypl(4) = 30.0_dp + zpl(4) = -1.0_dp + xpl(5) = 0.0_dp + ypl(5) = 0.0_dp + zpl(5) = -1.0_dp ! Separator - xpl(6) = dmiss; ypl(6) = dmiss; zpl(6) = dmiss + xpl(6) = dmiss + ypl(6) = dmiss + zpl(6) = dmiss ! Dry point polygon (rectangle from 10,10 to 20,20) with zpl=1 - xpl(7) = 10.0_dp; ypl(7) = 10.0_dp; zpl(7) = 1.0_dp - xpl(8) = 20.0_dp; ypl(8) = 10.0_dp; zpl(8) = 1.0_dp - xpl(9) = 20.0_dp; ypl(9) = 20.0_dp; zpl(9) = 1.0_dp - xpl(10) = 10.0_dp; ypl(10) = 20.0_dp; zpl(10) = 1.0_dp - xpl(11) = 10.0_dp; ypl(11) = 10.0_dp; zpl(11) = 1.0_dp + xpl(7) = 10.0_dp + ypl(7) = 10.0_dp + zpl(7) = 1.0_dp + xpl(8) = 20.0_dp + ypl(8) = 10.0_dp + zpl(8) = 1.0_dp + xpl(9) = 20.0_dp + ypl(9) = 20.0_dp + zpl(9) = 1.0_dp + xpl(10) = 10.0_dp + ypl(10) = 20.0_dp + zpl(10) = 1.0_dp + xpl(11) = 10.0_dp + ypl(11) = 10.0_dp + zpl(11) = 1.0_dp ! Separators - xpl(12) = dmiss; ypl(12) = dmiss; zpl(12) = dmiss - xpl(13) = dmiss; ypl(13) = dmiss; zpl(13) = dmiss - xpl(14) = dmiss; ypl(14) = dmiss; zpl(14) = dmiss - xpl(15) = dmiss; ypl(15) = dmiss; zpl(15) = dmiss + xpl(12) = dmiss + ypl(12) = dmiss + zpl(12) = dmiss + xpl(13) = dmiss + ypl(13) = dmiss + zpl(13) = dmiss + xpl(14) = dmiss + ypl(14) = dmiss + zpl(14) = dmiss + xpl(15) = dmiss + ypl(15) = dmiss + zpl(15) = dmiss ! Initialize polygon data structures call cellmask_from_polygon_set_init(NPL, xpl, ypl, zpl) @@ -122,25 +152,51 @@ subroutine test_nested_drypoint_polygons() bind(C) allocate(xpl(npl), ypl(npl), zpl(npl)) ! Outer dry point polygon (rectangle from 0,0 to 40,40) - xpl(1) = 0.0_dp; ypl(1) = 0.0_dp; zpl(1) = 1.0_dp - xpl(2) = 40.0_dp; ypl(2) = 0.0_dp; zpl(2) = 1.0_dp - xpl(3) = 40.0_dp; ypl(3) = 40.0_dp; zpl(3) = 1.0_dp - xpl(4) = 0.0_dp; ypl(4) = 40.0_dp; zpl(4) = 1.0_dp - xpl(5) = 0.0_dp; ypl(5) = 0.0_dp; zpl(5) = 1.0_dp + xpl(1) = 0.0_dp + ypl(1) = 0.0_dp + zpl(1) = 1.0_dp + xpl(2) = 40.0_dp + ypl(2) = 0.0_dp + zpl(2) = 1.0_dp + xpl(3) = 40.0_dp + ypl(3) = 40.0_dp + zpl(3) = 1.0_dp + xpl(4) = 0.0_dp + ypl(4) = 40.0_dp + zpl(4) = 1.0_dp + xpl(5) = 0.0_dp + ypl(5) = 0.0_dp + zpl(5) = 1.0_dp ! Separator - xpl(6) = dmiss; ypl(6) = dmiss; zpl(6) = dmiss + xpl(6) = dmiss + ypl(6) = dmiss + zpl(6) = dmiss ! Inner dry point polygon (rectangle from 10,10 to 30,30) - xpl(7) = 10.0_dp; ypl(7) = 10.0_dp; zpl(7) = 1.0_dp - xpl(8) = 30.0_dp; ypl(8) = 10.0_dp; zpl(8) = 1.0_dp - xpl(9) = 30.0_dp; ypl(9) = 30.0_dp; zpl(9) = 1.0_dp - xpl(10) = 10.0_dp; ypl(10) = 30.0_dp; zpl(10) = 1.0_dp - xpl(11) = 10.0_dp; ypl(11) = 10.0_dp; zpl(11) = 1.0_dp + xpl(7) = 10.0_dp + ypl(7) = 10.0_dp + zpl(7) = 1.0_dp + xpl(8) = 30.0_dp + ypl(8) = 10.0_dp + zpl(8) = 1.0_dp + xpl(9) = 30.0_dp + ypl(9) = 30.0_dp + zpl(9) = 1.0_dp + xpl(10) = 10.0_dp + ypl(10) = 30.0_dp + zpl(10) = 1.0_dp + xpl(11) = 10.0_dp + ypl(11) = 10.0_dp + zpl(11) = 1.0_dp ! Separators - xpl(12) = dmiss; ypl(12) = dmiss; zpl(12) = dmiss - xpl(13) = dmiss; ypl(13) = dmiss; zpl(13) = dmiss + xpl(12) = dmiss + ypl(12) = dmiss + zpl(12) = dmiss + xpl(13) = dmiss + ypl(13) = dmiss + zpl(13) = dmiss ! Initialize polygon data structures call cellmask_from_polygon_set_init(NPL, xpl, ypl, zpl) From c78b37ae005f22d1075d2ac9eba5db72b889ff83 Mon Sep 17 00:00:00 2001 From: FlorisBuwaldaDeltares Date: Mon, 8 Dec 2025 16:10:34 +0100 Subject: [PATCH 32/32] forgot to reset "enclosures_present" boolean --- .../src/dflowfm_kernel/prepost/m_cellmask_from_polygon_set.f90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/m_cellmask_from_polygon_set.f90 b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/m_cellmask_from_polygon_set.f90 index abdbefbf64..5b60cd5979 100644 --- a/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/m_cellmask_from_polygon_set.f90 +++ b/src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/m_cellmask_from_polygon_set.f90 @@ -201,6 +201,7 @@ subroutine cellmask_from_polygon_set_cleanup() polygons = 0 cellmask_initialized = .false. + enclosures_present = .false. end subroutine cellmask_from_polygon_set_cleanup