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
Merged

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

Show file tree
Hide file tree
Changes from 1 commit
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
51 changes: 48 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,59 @@ 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
wpbonelli marked this conversation as resolved.
Show resolved Hide resolved
subroutine get_polyverts(this, ic, polyverts, wrap)
! -- 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
logical(LGP), intent(in), optional :: wrap !< whether to wrap around (duplicating a vertex)
wpbonelli marked this conversation as resolved.
Show resolved Hide resolved
! -- local
integer(I4B) :: icu, nverts, irow, jcol, klay
real(DP) :: cellx, celly, dxhalf, dyhalf
logical(LGP) :: lwrap

! check wraparound option
if (.not. (present(wrap))) then
lwrap = .false.
nverts = 4
else
lwrap = wrap
if (lwrap) &
nverts = 5
end if

! allocate vertices array
allocate (polyverts(nverts, 2))
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should this be allocated and accessed as polyverts(2, nverts) to be consistent with Fortran column-major ordering? Probably doesn't matter a whole lot, but it would make sense to be consistent with the rest of the code.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

switched to column-major indexing for get_polyverts as well as for point_in_polygon and its tests


! 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

! wraparound if enabled
if (lwrap) &
polyverts(5, :) = polyverts(1, :)

end subroutine

subroutine read_int_array(this, line, lloc, istart, istop, iout, in, &
iarray, aname)
! ******************************************************************************
Expand Down
46 changes: 45 additions & 1 deletion src/Model/GroundWaterFlow/gwf3disv8.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
module GwfDisvModule

use ArrayReadersModule, only: ReadArray
use KindModule, only: DP, I4B
use KindModule, only: DP, I4B, LGP
use ConstantsModule, only: LINELENGTH, LENMEMPATH, LENVARNAME, DZERO, DONE, &
DHALF
use BaseDisModule, only: DisBaseType
Expand Down Expand Up @@ -49,6 +49,7 @@ module GwfDisvModule
procedure :: connection_vector
procedure :: supports_layers
procedure :: get_ncpl
procedure :: get_polyverts
! -- private
procedure :: source_options
procedure :: source_dimensions
Expand Down Expand Up @@ -1575,6 +1576,49 @@ function get_ncpl(this)
return
end function get_ncpl

!> @brief Get a 2D array of polygon vertices
subroutine get_polyverts(this, ic, polyverts, wrap)
! -- dummy
class(GwfDisvType), intent(inout) :: this
integer(I4B), intent(in) :: ic !< cell number (reduced)
real(DP), allocatable, intent(out) :: polyverts(:, :) !< polygon vertices
logical(LGP), intent(in), optional :: wrap !< whether to wrap around (duplicating a vertex)
! -- local
integer(I4B) :: icu, icu2d, iavert, ncpl, nverts, m, j
logical(LGP) :: lwrap

! count vertices
ncpl = this%get_ncpl()
icu = this%get_nodeuser(ic)
icu2d = icu - ((icu - 1) / ncpl) * ncpl
nverts = this%iavert(icu2d + 1) - this%iavert(icu2d) - 1
if (nverts .le. 0) nverts = nverts + size(this%javert)

! check wraparound option
if (.not. (present(wrap))) then
lwrap = .false.
else
lwrap = wrap
if (lwrap) &
nverts = nverts + 1
end if

! allocate vertices array
allocate (polyverts(nverts, 2))

! set vertices
iavert = this%iavert(icu2d)
do m = 1, nverts - 1
j = this%javert(iavert - 1 + m)
polyverts(m, :) = (/this%vertices(1, j), this%vertices(2, j)/)
end do

! wraparound if enabled
if (lwrap) &
polyverts(nverts, :) = polyverts(1, :)

end subroutine

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