refactor(disv): improve cell area calculation

* improve cell area calculation * close #1346
MODFLOW-USGS · Sep 9, 2023 · 0ca18e7 · 0ca18e7
1 parent efafa5c
commit 0ca18e7
Show file tree

Hide file tree

Showing 2 changed files with 16 additions and 4 deletions.
diff --git a/src/Model/GroundWaterFlow/gwf3disv8.f90 b/src/Model/GroundWaterFlow/gwf3disv8.f90
@@ -1325,21 +1325,27 @@ function get_cell2d_area(this, icell2d) result(area)
     integer(I4B) :: ivert
     integer(I4B) :: nvert
     integer(I4B) :: icount
+    integer(I4B) :: iv1
     real(DP) :: x
     real(DP) :: y
+    real(DP) :: x1
+    real(DP) :: y1
 ! ------------------------------------------------------------------------------
     !
     area = DZERO
     nvert = this%iavert(icell2d + 1) - this%iavert(icell2d)
     icount = 1
+    iv1 = this%javert(this%iavert(icell2d))
+    x1 = this%vertices(1, iv1)
+    y1 = this%vertices(2, iv1)
     do ivert = this%iavert(icell2d), this%iavert(icell2d + 1) - 1
       x = this%vertices(1, this%javert(ivert))
       if (icount < nvert) then
         y = this%vertices(2, this%javert(ivert + 1))
       else
         y = this%vertices(2, this%javert(this%iavert(icell2d)))
       end if
-      area = area + x * y
+      area = area + (x - x1) * (y - y1)
       icount = icount + 1
     end do
     !
@@ -1351,7 +1357,7 @@ function get_cell2d_area(this, icell2d) result(area)
       else
         x = this%vertices(1, this%javert(this%iavert(icell2d)))
       end if
-      area = area - x * y
+      area = area - (x - x1) * (y - y1)
       icount = icount + 1
     end do
     !

diff --git a/src/Model/ModelUtilities/DisvGeom.f90 b/src/Model/ModelUtilities/DisvGeom.f90
@@ -335,21 +335,27 @@ function get_area(this) result(area)
     integer(I4B) :: ivert
     integer(I4B) :: nvert
     integer(I4B) :: icount
+    integer(I4B) :: iv1
     real(DP) :: x
     real(DP) :: y
+    real(DP) :: x1
+    real(DP) :: y1
 ! ------------------------------------------------------------------------------
     !
     area = DZERO
     nvert = this%iavert(this%j + 1) - this%iavert(this%j)
     icount = 1
+    iv1 = this%javert(this%iavert(this%j))
+    x1 = this%vertex_grid(1, iv1)
+    y1 = this%vertex_grid(2, iv1)
     do ivert = this%iavert(this%j), this%iavert(this%j + 1) - 1
       x = this%vertex_grid(1, this%javert(ivert))
       if (icount < nvert) then
         y = this%vertex_grid(2, this%javert(ivert + 1))
       else
         y = this%vertex_grid(2, this%javert(this%iavert(this%j)))
       end if
-      area = area + x * y
+      area = area + (x - x1) * (y - y1)
       icount = icount + 1
     end do
     !
@@ -361,7 +367,7 @@ function get_area(this) result(area)
       else
         x = this%vertex_grid(1, this%javert(this%iavert(this%j)))
       end if
-      area = area - x * y
+      area = area - (x - x1) * (y - y1)
       icount = icount + 1
     end do
     !