All/task/unst 9570 fast dbpinpol

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

@@ -49,7 +49,8 @@ module m_find1dcells
 !>    it is assumed that kc has been allocated
 !>    it is assumed that findcells has already been called (for 2d cells)
    subroutine find1dcells()
-      use m_cellmask_from_polygon_set, only: init_cell_geom_as_polylines, point_find_netcell, cleanup_cell_geom_polylines
+      use m_cellmask_from_polygon_set, only: point_find_netcell, cleanup_cell_geom_polylines
+      use m_pol_to_cellmask, only: init_cell_geom_as_polylines
 
       implicit none
 

src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/m_pol_to_cellmask.F90

@@ -0,0 +1,159 @@
+ !----- AGPL --------------------------------------------------------------------
+!
+!  Copyright (C)  Stichting Deltares, 2017-2026.
+!
+!  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.
+!
+!-------------------------------------------------------------------------------
+
+!> 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
+
+   public :: pol_to_cellmask, init_cell_geom_as_polylines, find_cells_crossed_by_polyline
+
+
+contains
+
+   subroutine pol_to_cellmask()
+      use network_data, only: cellmask, nump1d2d, npl, nump, xzw, yzw, xpl, ypl, zpl
+      use m_alloc, only: realloc
+
+      integer :: k
+
+      if (NPL == 0) then
+         return
+      end if
+
+      call realloc(cellmask, nump1d2d, keepexisting=.false., fill=0)
+
+      call cellmask_from_polygon_set_init(NPL, xpl, ypl, zpl)
+
+      !> Dynamic scheduling in case of unequal work, chunksize guided
+      !$OMP PARALLEL DO SCHEDULE(GUIDED)
+      do k = 1, nump
+         cellmask(k) = cellmask_from_polygon_set(xzw(k), yzw(k))
+      end do
+      !$OMP END PARALLEL DO
+
+      call cellmask_from_polygon_set_cleanup()
+
+   end subroutine pol_to_cellmask
+
+!> Initialize xpl, ypl, zpl arrays with all netcell geometries (called once)
+   subroutine init_cell_geom_as_polylines()
+      use network_data
+      use m_alloc
+      USE m_cellmask_from_polygon_set, ONLY: cellmask_from_polygon_set_init, cellmask_initialized
+      use m_missing, only: dmiss
+
+      integer :: k, n, k1, total_points, ipoint
+
+      if (cellmask_initialized) then !> reuse cellmask cache boolean
+         return
+      end if
+
+      call savepol()
+
+      ! calculate total points needed: sum(netcell(k)%n + 1) for all cells
+      ! +1 for dmiss separator after each polygon
+      total_points = 0
+      do k = 1, nump
+         total_points = total_points + netcell(k)%n + 1 ! +1 for dmiss
+      end do
+
+      ! allocate or reallocate xpl, ypl, zpl
+      call realloc(xpl, total_points, keepexisting=.false.)
+      call realloc(ypl, total_points, keepexisting=.false.)
+      call realloc(zpl, total_points, keepexisting=.false.)
+
+      ! fill arrays with netcell geometry
+      ipoint = 0
+      do k = 1, nump
+         do n = 1, netcell(k)%n
+            ipoint = ipoint + 1
+            k1 = netcell(k)%nod(n)
+            xpl(ipoint) = xk(k1)
+            ypl(ipoint) = yk(k1)
+            zpl(ipoint) = real(k, dp) ! store cell index as z-value
+         end do
+
+         ! add separator
+         ipoint = ipoint + 1
+         xpl(ipoint) = dmiss
+         ypl(ipoint) = dmiss
+         zpl(ipoint) = dmiss
+      end do
+
+      npl = ipoint
+
+      ! initialize the cellmask module with these polygons
+      ! this builds bounding boxes and polygon indices
+      call cellmask_from_polygon_set_init(npl, xpl, ypl, zpl)
+
+   end subroutine init_cell_geom_as_polylines
+
+!> Find all cells crossed by polyline using brute force on cached geometry. The routine is inclusive of edge cases (touching edges or vertices).
+   subroutine find_cells_crossed_by_polyline(xpoly, ypoly, crossed_cells, error)
+      use m_alloc, only: realloc
+      use network_data, only: cellmask, nump
+      use m_missing, only: dmiss
+      use m_cellmask_from_polygon_set, only: cellmask_initialized, find_cells_for_segment
+
+      real(kind=dp), intent(in) :: xpoly(:) !< Polyline x-coordinates
+      real(kind=dp), intent(in) :: ypoly(:) !< Polyline y-coordinates
+      integer, allocatable, intent(out) :: crossed_cells(:) !> Indices of crossed cells in network_data::netcells
+      character, dimension(:), allocatable, intent(out) :: error !> Error message, empty if no error, to be handled at call site
+
+      integer :: npoly, i
+
+      error = ''
+
+      npoly = size(xpoly)
+      if (any(xpoly == dmiss) .or. any(ypoly == dmiss) .or. npoly < 2) then
+         error = 'Invalid polyline input'
+         return
+      end if
+
+      if (.not. cellmask_initialized) then
+         call init_cell_geom_as_polylines()
+      end if
+
+      call realloc(cellmask, nump, keepexisting=.false., fill=0)
+
+      ! Process each segment and put the result in cellmask
+      do i = 1, npoly - 1
+         call find_cells_for_segment(xpoly(i), ypoly(i), xpoly(i + 1), ypoly(i + 1), cellmask)
+      end do
+
+      crossed_cells = pack([(i, i=1, nump)], mask=(cellmask == 1))
+
+   end subroutine find_cells_crossed_by_polyline
+
+end module m_pol_to_cellmask

src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_gui/m_samples_to_cellmask.f90 → src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/prepost/m_samples_to_cellmask.f90

@@ -39,7 +39,8 @@ subroutine samples_to_cellmask()
 
       use network_data, only: cellmask, nump1d2d
       use m_samples, only: ns, xs, ys
-      use m_cellmask_from_polygon_set, only: init_cell_geom_as_polylines, point_find_netcell, cleanup_cell_geom_polylines
+      use m_cellmask_from_polygon_set, only: point_find_netcell, cleanup_cell_geom_polylines
+      use m_pol_to_cellmask, only: init_cell_geom_as_polylines
       use m_alloc, only: realloc
 
       integer :: i, k

src/engines_gpl/dflowfm/packages/dflowfm_kernel/src/dflowfm_kernel/timespace/fm_external_forcings.f90

@@ -2509,7 +2509,8 @@ subroutine finalize()
       use m_filez, only: doclose
       use m_physcoef, only: constant_dicoww, dicoww
       use m_array_or_scalar, only: realloc
-      use m_cellmask_from_polygon_set, only: init_cell_geom_as_polylines, point_find_netcell, cleanup_cell_geom_polylines
+      use m_cellmask_from_polygon_set, only: point_find_netcell, cleanup_cell_geom_polylines
+      use m_pol_to_cellmask, only: init_cell_geom_as_polylines
 
       integer :: j, k, ierr, l, n, itp, kk, k1, k2, kb, kt, nstor, i, ja
       integer :: imba, needextramba, needextrambar

src/engines_gpl/dflowfm/packages/dflowfm_kernel/test/test_pol_to_cellmask.f90

@@ -3,7 +3,8 @@ module test_pol_to_cellmask
    use precision, only: dp
    use m_missing, only: dmiss
    use network_data, only: cellmask, npl, nump, xzw, yzw, xpl, ypl, zpl
-   use m_cellmask_from_polygon_set, only: cellmask_from_polygon_set_init, cellmask_from_polygon_set, cellmask_from_polygon_set_cleanup
+   use m_cellmask_from_polygon_set
+   use m_pol_to_cellmask
    use geometry_module, only: pinpok_legacy, pinpok_raycast
 
    implicit none(external)
@@ -406,7 +407,8 @@ end subroutine test_pinpok_raycast_complex
    subroutine test_incells_basic_functionality() bind(C)
       ! Test basic incells functionality: point inside/outside netcells
       use gridoperations, only: incells
-      use m_cellmask_from_polygon_set, only: init_cell_geom_as_polylines, point_find_netcell, cleanup_cell_geom_polylines
+      use m_cellmask_from_polygon_set, only: point_find_netcell, cleanup_cell_geom_polylines
+      use m_pol_to_cellmask, only: init_cell_geom_as_polylines
       use network_data, only: netcell, nump, xk, yk
       use m_alloc, only: realloc
 
@@ -471,7 +473,8 @@ end subroutine test_incells_basic_functionality
    subroutine test_incells_complex_geometry() bind(C)
       ! Test incells with more complex netcell geometries (triangles, pentagons, hexagons)
       use gridoperations, only: incells
-      use m_cellmask_from_polygon_set, only: init_cell_geom_as_polylines, point_find_netcell, cleanup_cell_geom_polylines
+      use m_cellmask_from_polygon_set, only: point_find_netcell, cleanup_cell_geom_polylines
+      use m_pol_to_cellmask, only: init_cell_geom_as_polylines
       use network_data, only: netcell, nump, xk, yk
 
       integer :: kin_old, kin_new
@@ -524,7 +527,8 @@ end subroutine test_incells_complex_geometry
    subroutine test_incells_large_grid() bind(C)
       ! Test incells with a larger grid (performance sanity check)
       use gridoperations, only: incells
-      use m_cellmask_from_polygon_set, only: init_cell_geom_as_polylines, point_find_netcell, cleanup_cell_geom_polylines
+      use m_cellmask_from_polygon_set, only: point_find_netcell, cleanup_cell_geom_polylines
+      use m_pol_to_cellmask, only: init_cell_geom_as_polylines
       use network_data, only: netcell, nump
 
       integer :: kin_old, kin_new, i, mismatches
@@ -566,7 +570,8 @@ end subroutine test_incells_large_grid
    subroutine test_incells_cache_consistency() bind(C)
       ! Test that cache initialization produces consistent results
       use gridoperations, only: incells
-      use m_cellmask_from_polygon_set, only: init_cell_geom_as_polylines, point_find_netcell, cleanup_cell_geom_polylines
+      use m_cellmask_from_polygon_set, only: point_find_netcell, cleanup_cell_geom_polylines
+      use m_pol_to_cellmask, only: init_cell_geom_as_polylines
       use network_data, only: nump
 
       integer :: kin1, kin2, kin_old
@@ -614,7 +619,8 @@ end subroutine test_incells_cache_consistency
    subroutine test_incells_edge_cases() bind(C)
       ! Test edge cases: empty grid, single cell, point at infinity
       use gridoperations, only: incells
-      use m_cellmask_from_polygon_set, only: init_cell_geom_as_polylines, point_find_netcell, cleanup_cell_geom_polylines
+      use m_cellmask_from_polygon_set, only: point_find_netcell, cleanup_cell_geom_polylines
+      use m_pol_to_cellmask, only: init_cell_geom_as_polylines
       use network_data, only: netcell, nump
 
       integer :: kin_old, kin_new
@@ -819,7 +825,6 @@ end subroutine test_samples_to_cellmask_edge_cases
 !$f90tw TESTCODE(TEST, test_pol_to_cellmask, test_find_cells_crossed_by_polyline_simple, test_find_cells_crossed_by_polyline_simple,
    subroutine test_find_cells_crossed_by_polyline_simple() bind(C)
       ! Test find_cells_crossed_by_polyline: polyline crosses some cells, misses others
-      use m_cellmask_from_polygon_set, only: init_cell_geom_as_polylines, find_cells_crossed_by_polyline, cleanup_cell_geom_polylines
       use network_data, only: nump
       use m_alloc, only: realloc
 
@@ -888,7 +893,6 @@ end subroutine test_find_cells_crossed_by_polyline_simple
 !$f90tw TESTCODE(TEST, test_pol_to_cellmask, test_find_cells_crossed_by_polyline_edge_cases, test_find_cells_crossed_by_polyline_edge_cases,
    subroutine test_find_cells_crossed_by_polyline_edge_cases() bind(C)
       ! Test find_cells_crossed_by_polyline: polyline crosses some cells, misses others
-      use m_cellmask_from_polygon_set, only: init_cell_geom_as_polylines, find_cells_crossed_by_polyline, cleanup_cell_geom_polylines
       use network_data, only: nump
       use m_alloc, only: realloc
 
@@ -952,6 +956,50 @@ end subroutine test_find_cells_crossed_by_polyline_edge_cases
 ! Helper subroutines for setting up test geometries
 ! ============================================================================
 
+!$f90tw TESTCODE(TEST, test_pol_to_cellmask, test_dbpinpol_basic_functionality, test_dbpinpol_basic_functionality,
+   subroutine test_dbpinpol_basic_functionality() bind(C)
+      ! Test dbpinpol with known inputs/outputs (no legacy comparison)
+      use geometry_module, only: dbpinpol
+      use m_polygon, only: xpl, ypl, zpl, npl
+      use m_alloc, only: realloc
+
+      real(kind=dp) :: xp, yp
+      integer :: in_result
+
+      ! Setup simple rectangular polygon
+      npl = 5
+      call realloc(xpl, npl, keepexisting=.false.)
+      call realloc(ypl, npl, keepexisting=.false.)
+      call realloc(zpl, npl, keepexisting=.false.)
+
+      xpl = [0.0_dp, 10.0_dp, 10.0_dp, 0.0_dp, 0.0_dp]
+      ypl = [0.0_dp, 0.0_dp, 10.0_dp, 10.0_dp, 0.0_dp]
+      zpl = [1.0_dp, 1.0_dp, 1.0_dp, 1.0_dp, 1.0_dp]
+
+      ! Initialize AND test first point
+      in_result = -1
+      xp = 0.0_dp
+      yp = 0.0_dp
+      call dbpinpol(xp, yp, in_result, dmiss, 1, npl, xpl, ypl, zpl)
+      call f90_expect_eq(in_result, 1, "Point (0,0) should be inside polygon")
+
+      ! Test second point (reusing cached data)
+      xp = 5.0_dp
+      yp = 5.0_dp
+      call dbpinpol(xp, yp, in_result, dmiss, 1, npl, xpl, ypl, zpl)
+      call f90_expect_eq(in_result, 1, "Point (5,5) should be inside polygon")
+
+      ! Test outside point
+      xp = 15.0_dp
+      yp = 15.0_dp
+      call dbpinpol(xp, yp, in_result, dmiss, 1, npl, xpl, ypl, zpl)
+      call f90_expect_eq(in_result, 0, "Point (15,15) should be outside polygon")
+
+      deallocate (xpl, ypl, zpl)
+
+   end subroutine
+!$f90tw)
+
    ! Helper subroutine: setup 5 cells in a row (0-10, 10-20, 20-30, 30-40, 40-50)
    subroutine setup_row_netcells()
       use network_data, only: netcell, nump, xk, yk, numk, nump1d2d