Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

feature(discretization): add get_polyverts routine to DIS/DISV #1444

Merged
merged 2 commits into from
Nov 15, 2023
+386 −663
Merged

feature(discretization): add get_polyverts routine to DIS/DISV #1444

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
5595394
feature(discretization): add get_polyverts routine to DIS/DISV
wpbonelli Nov 14, 2023
ab00c9f
revisions:
wpbonelli Nov 15, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
198 changes: 100 additions & 98 deletions autotest/TestGeomUtil.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -23,6 +23,8 @@ subroutine collect_geomutil(testsuite)
]
end subroutine collect_geomutil

! 2D arrays for polygons and check points use column-major indexing

subroutine test_point_in_polygon(error, shape, &
poly, in_pts, out_pts, vert_pts, face_pts)
type(error_type), allocatable, intent(inout) :: error
Expand All @@ -36,39 +38,39 @@ subroutine test_point_in_polygon(error, shape, &
real(DP) :: x, y

! test inside points
do i = 1, size(in_pts, 1)
x = in_pts(i, 1)
y = in_pts(i, 2)
do i = 1, size(in_pts, 2)
x = in_pts(1, i)
y = in_pts(2, i)
call check(error, point_in_polygon(x, y, poly), &
"point inside "//shape//" failed: " &
//to_string(x)//", "//to_string(y))
if (allocated(error)) return
end do

! test outside points
do i = 1, size(out_pts, 1)
x = out_pts(i, 1)
y = out_pts(i, 2)
do i = 1, size(out_pts, 2)
x = out_pts(1, i)
y = out_pts(2, i)
call check(error, (.not. point_in_polygon(x, y, poly)), &
"point outside "//shape//" failed: " &
//to_string(x)//", "//to_string(y))
if (allocated(error)) return
end do

! test vertex points
do i = 1, size(vert_pts, 1)
x = vert_pts(i, 1)
y = vert_pts(i, 2)
do i = 1, size(vert_pts, 2)
x = vert_pts(1, i)
y = vert_pts(2, i)
call check(error, point_in_polygon(x, y, poly), &
"point on "//shape//" vertex failed: " &
//to_string(x)//", "//to_string(y))
if (allocated(error)) return
end do

! test face points
do i = 1, size(face_pts, 1)
x = face_pts(i, 1)
y = face_pts(i, 2)
do i = 1, size(face_pts, 2)
x = face_pts(1, i)
y = face_pts(2, i)
call check(error, point_in_polygon(x, y, poly), &
"point on "//shape//" face failed: " &
//to_string(x)//", "//to_string(y))
Expand All @@ -85,41 +87,41 @@ subroutine test_point_in_polygon_sq(error)
real(DP), allocatable :: vert_pts(:, :)
real(DP), allocatable :: face_pts(:, :)

allocate (poly(4, 2))
allocate (poly(2, 4))

allocate (in_pts(3, 2))
in_pts(1, :) = (/0.99_DP, 0.01_DP/)
in_pts(2, :) = (/0.5_DP, 0.5_DP/)
in_pts(3, :) = (/0.0001_DP, 0.9999_DP/)
allocate (in_pts(2, 3))
in_pts(:, 1) = (/0.99_DP, 0.01_DP/)
in_pts(:, 2) = (/0.5_DP, 0.5_DP/)
in_pts(:, 3) = (/0.0001_DP, 0.9999_DP/)

allocate (out_pts(2, 2))
out_pts(1, :) = (/0.5_DP, 1.00001_DP/)
out_pts(2, :) = (/-0.5_DP, 34.0_DP/)

allocate (vert_pts(4, 2))
vert_pts(1, :) = (/0.0_DP, 0.0_DP/)
vert_pts(2, :) = (/1.0_DP, 0.0_DP/)
vert_pts(3, :) = (/0.0_DP, 1.0_DP/)
vert_pts(4, :) = (/1.0_DP, 1.0_DP/)

allocate (face_pts(4, 2))
face_pts(1, :) = (/0.0_DP, 0.5_DP/)
face_pts(2, :) = (/0.5_DP, 0.0_DP/)
face_pts(3, :) = (/1.0_DP, 0.5_DP/)
face_pts(4, :) = (/0.5_DP, 1.0_DP/)

poly(1, :) = (/0.0_DP, 0.0_DP/)
poly(2, :) = (/0.0_DP, 1.0_DP/)
poly(3, :) = (/1.0_DP, 1.0_DP/)
poly(4, :) = (/1.0_DP, 0.0_DP/)
out_pts(:, 1) = (/0.5_DP, 1.00001_DP/)
out_pts(:, 2) = (/-0.5_DP, 34.0_DP/)

allocate (vert_pts(2, 4))
vert_pts(:, 1) = (/0.0_DP, 0.0_DP/)
vert_pts(:, 2) = (/1.0_DP, 0.0_DP/)
vert_pts(:, 3) = (/0.0_DP, 1.0_DP/)
vert_pts(:, 4) = (/1.0_DP, 1.0_DP/)

allocate (face_pts(2, 4))
face_pts(:, 1) = (/0.0_DP, 0.5_DP/)
face_pts(:, 2) = (/0.5_DP, 0.0_DP/)
face_pts(:, 3) = (/1.0_DP, 0.5_DP/)
face_pts(:, 4) = (/0.5_DP, 1.0_DP/)

poly(:, 1) = (/0.0_DP, 0.0_DP/)
poly(:, 2) = (/0.0_DP, 1.0_DP/)
poly(:, 3) = (/1.0_DP, 1.0_DP/)
poly(:, 4) = (/1.0_DP, 0.0_DP/)
call test_point_in_polygon(error, "clockwise square", &
poly, in_pts, out_pts, vert_pts, face_pts)
if (allocated(error)) return

poly(1, :) = (/0.0_DP, 0.0_DP/)
poly(2, :) = (/1.0_DP, 0.0_DP/)
poly(3, :) = (/1.0_DP, 1.0_DP/)
poly(4, :) = (/0.0_DP, 1.0_DP/)
poly(:, 1) = (/0.0_DP, 0.0_DP/)
poly(:, 2) = (/1.0_DP, 0.0_DP/)
poly(:, 3) = (/1.0_DP, 1.0_DP/)
poly(:, 4) = (/0.0_DP, 1.0_DP/)
call test_point_in_polygon(error, "counter-clockwise square", &
poly, in_pts, out_pts, vert_pts, face_pts)
if (allocated(error)) return
Expand All @@ -140,37 +142,37 @@ subroutine test_point_in_polygon_tri(error)
real(DP), allocatable :: vert_pts(:, :)
real(DP), allocatable :: face_pts(:, :)

allocate (poly(3, 2))
allocate (poly(2, 3))

allocate (in_pts(3, 2))
in_pts(1, :) = (/0.8_DP, 0.0001_DP/)
in_pts(2, :) = (/0.5_DP, 0.49999_DP/)
in_pts(3, :) = (/0.0001_DP, 0.8_DP/)
allocate (in_pts(2, 3))
in_pts(:, 1) = (/0.8_DP, 0.0001_DP/)
in_pts(:, 2) = (/0.5_DP, 0.49999_DP/)
in_pts(:, 3) = (/0.0001_DP, 0.8_DP/)

allocate (out_pts(2, 2))
out_pts(1, :) = (/0.5_DP, 0.50001_DP/)
out_pts(2, :) = (/-0.5_DP, 34.0_DP/)

allocate (vert_pts(3, 2))
vert_pts(1, :) = (/0.0_DP, 0.0_DP/)
vert_pts(2, :) = (/1.0_DP, 0.0_DP/)
vert_pts(3, :) = (/0.0_DP, 1.0_DP/)

allocate (face_pts(3, 2))
face_pts(1, :) = (/0.0_DP, 0.5_DP/)
face_pts(2, :) = (/0.5_DP, 0.0_DP/)
face_pts(3, :) = (/0.5_DP, 0.5_DP/)

poly(1, :) = (/0.0_DP, 0.0_DP/)
poly(2, :) = (/0.0_DP, 1.0_DP/)
poly(3, :) = (/1.0_DP, 0.0_DP/)
out_pts(:, 1) = (/0.5_DP, 0.50001_DP/)
out_pts(:, 2) = (/-0.5_DP, 34.0_DP/)

allocate (vert_pts(2, 3))
vert_pts(:, 1) = (/0.0_DP, 0.0_DP/)
vert_pts(:, 2) = (/1.0_DP, 0.0_DP/)
vert_pts(:, 3) = (/0.0_DP, 1.0_DP/)

allocate (face_pts(2, 3))
face_pts(:, 1) = (/0.0_DP, 0.5_DP/)
face_pts(:, 2) = (/0.5_DP, 0.0_DP/)
face_pts(:, 3) = (/0.5_DP, 0.5_DP/)

poly(:, 1) = (/0.0_DP, 0.0_DP/)
poly(:, 2) = (/0.0_DP, 1.0_DP/)
poly(:, 3) = (/1.0_DP, 0.0_DP/)
call test_point_in_polygon(error, "clockwise triangle", &
poly, in_pts, out_pts, vert_pts, face_pts)
if (allocated(error)) return

poly(1, :) = (/0.0_DP, 0.0_DP/)
poly(2, :) = (/1.0_DP, 0.0_DP/)
poly(3, :) = (/0.0_DP, 1.0_DP/)
poly(:, 1) = (/0.0_DP, 0.0_DP/)
poly(:, 2) = (/1.0_DP, 0.0_DP/)
poly(:, 3) = (/0.0_DP, 1.0_DP/)
call test_point_in_polygon(error, "counter-clockwise triangle", &
poly, in_pts, out_pts, vert_pts, face_pts)
if (allocated(error)) return
Expand All @@ -191,45 +193,45 @@ subroutine test_point_in_polygon_irr(error)
real(DP), allocatable :: vert_pts(:, :)
real(DP), allocatable :: face_pts(:, :)

allocate (poly(5, 2))

allocate (in_pts(3, 2))
in_pts(1, :) = (/0.5_DP, 0.1_DP/)
in_pts(2, :) = (/0.5_DP, 0.49_DP/)
in_pts(3, :) = (/1.999_DP, 1.999_DP/)

allocate (out_pts(3, 2))
out_pts(1, :) = (/0.5_DP, -0.1_DP/)
out_pts(2, :) = (/0.5_DP, 0.51_DP/)
out_pts(3, :) = (/-0.5_DP, 34.0_DP/)

allocate (vert_pts(5, 2))
vert_pts(1, :) = (/0.0_DP, 0.0_DP/)
vert_pts(2, :) = (/1.0_DP, 1.0_DP/)
vert_pts(3, :) = (/1.0_DP, 2.0_DP/)
vert_pts(4, :) = (/2.0_DP, 2.0_DP/)
vert_pts(5, :) = (/2.0_DP, 0.0_DP/)

allocate (face_pts(3, 2))
face_pts(1, :) = (/0.5_DP, 0.5_DP/)
face_pts(2, :) = (/2.0_DP, 1.0_DP/)
face_pts(3, :) = (/1.5_DP, 2.0_DP/)

poly(1, :) = (/0.0_DP, 0.0_DP/)
poly(2, :) = (/1.0_DP, 1.0_DP/)
poly(3, :) = (/1.0_DP, 2.0_DP/)
poly(4, :) = (/2.0_DP, 2.0_DP/)
poly(5, :) = (/2.0_DP, 0.0_DP/)
allocate (poly(2, 5))

allocate (in_pts(2, 3))
in_pts(:, 1) = (/0.5_DP, 0.1_DP/)
in_pts(:, 2) = (/0.5_DP, 0.49_DP/)
in_pts(:, 3) = (/1.999_DP, 1.999_DP/)

allocate (out_pts(2, 3))
out_pts(:, 1) = (/0.5_DP, -0.1_DP/)
out_pts(:, 2) = (/0.5_DP, 0.51_DP/)
out_pts(:, 3) = (/-0.5_DP, 34.0_DP/)

allocate (vert_pts(2, 5))
vert_pts(:, 1) = (/0.0_DP, 0.0_DP/)
vert_pts(:, 2) = (/1.0_DP, 1.0_DP/)
vert_pts(:, 3) = (/1.0_DP, 2.0_DP/)
vert_pts(:, 4) = (/2.0_DP, 2.0_DP/)
vert_pts(:, 5) = (/2.0_DP, 0.0_DP/)

allocate (face_pts(2, 3))
face_pts(:, 1) = (/0.5_DP, 0.5_DP/)
face_pts(:, 2) = (/2.0_DP, 1.0_DP/)
face_pts(:, 3) = (/1.5_DP, 2.0_DP/)

poly(:, 1) = (/0.0_DP, 0.0_DP/)
poly(:, 2) = (/1.0_DP, 1.0_DP/)
poly(:, 3) = (/1.0_DP, 2.0_DP/)
poly(:, 4) = (/2.0_DP, 2.0_DP/)
poly(:, 5) = (/2.0_DP, 0.0_DP/)
call test_point_in_polygon(error, &
"clockwise irregular polygon", &
poly, in_pts, out_pts, vert_pts, face_pts)
if (allocated(error)) return

poly(1, :) = (/0.0_DP, 0.0_DP/)
poly(2, :) = (/2.0_DP, 0.0_DP/)
poly(3, :) = (/2.0_DP, 2.0_DP/)
poly(4, :) = (/1.0_DP, 2.0_DP/)
poly(5, :) = (/1.0_DP, 1.0_DP/)
poly(:, 1) = (/0.0_DP, 0.0_DP/)
poly(:, 2) = (/2.0_DP, 0.0_DP/)
poly(:, 3) = (/2.0_DP, 2.0_DP/)
poly(:, 4) = (/1.0_DP, 2.0_DP/)
poly(:, 5) = (/1.0_DP, 1.0_DP/)
call test_point_in_polygon(error, &
"counter-clockwise irregular polygon", &
poly, in_pts, out_pts, vert_pts, face_pts)
Expand Down
55 changes: 52 additions & 3 deletions src/Model/GroundWaterFlow/gwf3dis8.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
module GwfDisModule

use ArrayReadersModule, only: ReadArray
use KindModule, only: DP, I4B
use KindModule, only: DP, I4B, LGP
use ConstantsModule, only: LINELENGTH, DHALF, DZERO, LENMEMPATH, LENVARNAME
use BaseDisModule, only: DisBaseType
use InputOutputModule, only: get_node, URWORD, ulasav, ulaprufw, ubdsv1, &
Expand Down Expand Up @@ -45,6 +45,7 @@ module GwfDisModule
procedure :: nodeu_from_cellid
procedure :: supports_layers
procedure :: get_ncpl
procedure :: get_polyverts
procedure :: connection_vector
procedure :: connection_normal
! -- private
Expand Down Expand Up @@ -1308,15 +1309,63 @@ subroutine connection_vector(this, noden, nodem, nozee, satn, satm, ihc, &
return
end subroutine

! return discretization type
!> @brief return discretization type
subroutine get_dis_type(this, dis_type)
class(GwfDisType), intent(in) :: this
character(len=*), intent(out) :: dis_type

dis_type = "DIS"

end subroutine get_dis_type

!> @brief Get a 2D array of polygon vertices, listed in
!! clockwise order beginning with the lower left corner
subroutine get_polyverts(this, ic, polyverts, closed)
! -- modules
use InputOutputModule, only: get_ijk
! -- dummy
class(GwfDisType), intent(inout) :: this
integer(I4B), intent(in) :: ic !< cell number (reduced)
real(DP), allocatable, intent(out) :: polyverts(:, :) !< polygon vertices (column-major indexing)
logical(LGP), intent(in), optional :: closed !< whether to close the polygon, duplicating a vertex
! -- local
integer(I4B) :: icu, nverts, irow, jcol, klay
real(DP) :: cellx, celly, dxhalf, dyhalf
logical(LGP) :: lclosed

nverts = 4

! check closed option
if (.not. (present(closed))) then
lclosed = .false.
else
lclosed = closed
end if

! allocate vertices array
if (lclosed) then
allocate (polyverts(2, nverts + 1))
else
allocate (polyverts(2, nverts))
end if

! set vertices
icu = this%get_nodeuser(ic)
call get_ijk(icu, this%nrow, this%ncol, this%nlay, irow, jcol, klay)
cellx = this%cellx(jcol)
celly = this%celly(irow)
dxhalf = DHALF * this%delr(jcol)
dyhalf = DHALF * this%delc(irow)
polyverts(:, 1) = (/cellx - dxhalf, celly - dyhalf/) ! SW
polyverts(:, 2) = (/cellx - dxhalf, celly + dyhalf/) ! NW
polyverts(:, 3) = (/cellx + dxhalf, celly + dyhalf/) ! NE
polyverts(:, 4) = (/cellx + dxhalf, celly - dyhalf/) ! SE

! close if enabled
if (lclosed) &
polyverts(:, nverts + 1) = polyverts(:, 1)

end subroutine

subroutine read_int_array(this, line, lloc, istart, istop, iout, in, &
iarray, aname)
! ******************************************************************************
Expand Down
Loading
Loading