Fm/task/unst 9476 pol to cellmask

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

src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/dbpinpol_cellmask.f90

@@ -0,0 +1,261 @@
+!----- 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 <http://www.gnu.org/licenses/>.
+!
+!  contact: [email protected]
+!  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
+   use m_polygon, only: xpl, ypl, npl
+
+   implicit none
+
+   private
+
+   !> dbpinpol routines are public to avoid PetSC dependency in unit tests
+   public :: dbpinpol_cellmask_init, dbpinpol_cellmask_cleanup, 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 :: 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
+
+      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
+
+      if (cellmask_initialized) then
+         call dbpinpol_cellmask_cleanup()
+      end if
+
+      if (NPL == 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.)
+         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 dbpinpol_cellmask_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)
+
+      integer :: mask
+      real(kind=dp), intent(in) :: xp, yp !< Point coordinates
+
+      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
+      found_inside_enclosure = .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)) 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 .and. in_test == 1) 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 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 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(i) == dmiss) exit
+
+         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
+            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

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))

src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/pol_to_cellmask.f90 → 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,16 +27,14 @@
 !
 !-------------------------------------------------------------------------------
 
-!
-!
-
 !> 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
@@ -47,47 +45,38 @@ module m_pol_to_cellmask
 
    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
-
-      integer :: in, k, KMOD
-      integer :: num
-
-      if (allocated(cellmask)) then
-         deallocate (cellmask)
-      end if
+      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 ! no polygon
-
-      call READYY('Applying polygon cellmask', 0.0_dp)
-      KMOD = max(1, NUMP / 100)
-
-      !  generate cell mask
-      in = -1
+      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
-!     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)
-
-         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
-         end if
-
-         if (in > 0) then
-!        mask cell
-            cellmask(k) = 1
-         end if
-         if (mod(k, KMOD) == 0) then
-            call READYY(' ', min(1.0_dp, real(k, kind=dp) / nump))
-         end if
+         cellmask(k) = dbpinpol_cellmask(xzw(k), yzw(k))
       end do
-      call READYY(' ', -1.0_dp)
+      !$OMP END PARALLEL DO
+#ifdef _OPENMP
+      call omp_set_num_threads(temp_threads)
+#endif
+
+      ! Cleanup
+      call dbpinpol_cellmask_cleanup()
 
-      return
    end subroutine pol_to_cellmask
 
 end module m_pol_to_cellmask

src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/CMakeLists.txt

@@ -15,3 +15,16 @@ 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
+)
+
+# 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)

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"