fix(get_disu_kwargs): incorrect indexing of delr and delc (#2011)

* fix(get_disu_kwargs): incorrect indexing of delr and delc

* add vertices to the return of get_disu_kwargs

* ensure nvert has proper value

* handle delr/delc and ensure size is correct

* isort

* Modify to allow mf6 autotest to pass

autotest/test_gridutil.py

@@ -82,19 +82,28 @@ def test_get_lni_infers_layer_count_when_int_ncpl(ncpl, nodes, expected):
             np.array([-10]),
             np.array([-30.0, -50.0]),
         ),
+        (
+            1,  # nlay
+            3,  # nrow
+            4,  # ncol
+            np.array(4 * [4.]),  # delr
+            np.array(3 * [3.]),  # delc
+            np.array([-10]),  # top
+            np.array([-30.0]),  # botm
+        ),
     ],
 )
 def test_get_disu_kwargs(nlay, nrow, ncol, delr, delc, tp, botm):
     kwargs = get_disu_kwargs(
-        nlay=nlay, nrow=nrow, ncol=ncol, delr=delr, delc=delc, tp=tp, botm=botm
+        nlay=nlay, nrow=nrow, ncol=ncol, delr=delr, delc=delc, tp=tp, botm=botm, return_vertices=True
     )
 
     from pprint import pprint
 
     pprint(kwargs["area"])
 
     assert kwargs["nodes"] == nlay * nrow * ncol
-    assert kwargs["nvert"] == None
+    assert kwargs["nvert"] == (nrow + 1) * (ncol + 1)
 
     area = np.array([dr * dc for (dr, dc) in product(delr, delc)], dtype=float)
     area = np.array(nlay * [area]).flatten()

flopy/utils/gridutil.py

@@ -6,7 +6,7 @@
 
 import numpy as np
 
-from .cvfdutil import get_disv_gridprops
+from .cvfdutil import centroid_of_polygon, get_disv_gridprops
 
 
 def get_lni(ncpl, nodes) -> List[Tuple[int, int]]:
@@ -68,6 +68,7 @@ def get_disu_kwargs(
     delc,
     tp,
     botm,
+    return_vertices=False,
 ):
     """
     Create args needed to construct a DISU package for a regular
@@ -89,11 +90,20 @@ def get_disu_kwargs(
         Top elevation(s) of cells in the model's top layer
     botm : numpy.ndarray
         Bottom elevation(s) for each layer
+    return_vertices: bool
+        If true, then include vertices and cell2d in kwargs
     """
 
     def get_nn(k, i, j):
         return k * nrow * ncol + i * ncol + j
 
+    if not isinstance(delr, np.ndarray):
+        delr = np.array(delr)
+    if not isinstance(delc, np.ndarray):
+        delc = np.array(delc)
+    assert delr.shape == (ncol,)
+    assert delc.shape == (nrow,)
+
     nodes = nlay * nrow * ncol
     iac = np.zeros((nodes), dtype=int)
     ja = []
@@ -110,8 +120,8 @@ def get_nn(k, i, j):
                 n = get_nn(k, i, j)
                 ja.append(n)
                 iac[n] += 1
-                area[n] = delr[i] * delc[j]
-                ihc.append(n + 1)
+                area[n] = delr[j] * delc[i]
+                ihc.append(k + 1)  # put layer in diagonal for flopy plotting
                 cl12.append(n + 1)
                 hwva.append(n + 1)
                 if k == 0:
@@ -126,7 +136,7 @@ def get_nn(k, i, j):
                     ihc.append(0)
                     dz = botm[k - 1] - botm[k]
                     cl12.append(0.5 * dz)
-                    hwva.append(delr[i] * delc[j])
+                    hwva.append(delr[j] * delc[i])
                 # back
                 if i > 0:
                     ja.append(get_nn(k, i - 1, j))
@@ -165,14 +175,45 @@ def get_nn(k, i, j):
                     else:
                         dz = botm[k - 1] - botm[k]
                     cl12.append(0.5 * dz)
-                    hwva.append(delr[i] * delc[j])
+                    hwva.append(delr[j] * delc[i])
     ja = np.array(ja, dtype=int)
     nja = ja.shape[0]
     hwva = np.array(hwva, dtype=float)
+
+    # build vertices
+    nvert = None
+    if return_vertices:
+        xv = np.cumsum(delr)
+        xv = np.array([0] + list(xv))
+        ymax = delc.sum()
+        yv = np.cumsum(delc)
+        yv = ymax - np.array([0] + list(yv))
+        xmg, ymg = np.meshgrid(xv, yv)
+        nvert = xv.shape[0] * yv.shape[0]
+        verts = np.array(list(zip(xmg.flatten(), ymg.flatten())))
+        vertices = []
+        for i in range(nvert):
+            vertices.append((i, verts[i, 0], verts[i, 1]))
+
+        cell2d = []
+        icell = 0
+        for k in range(nlay):
+            for i in range(nrow):
+                for j in range(ncol):
+                    iv0 = j + i * (ncol + 1)  # upper left vertex
+                    iv1 = iv0 + 1  # upper right vertex
+                    iv3 = iv0 + ncol + 1  # lower left vertex
+                    iv2 = iv3 + 1  # lower right vertex
+                    iverts = [iv0, iv1, iv2, iv3]
+                    vlist = [(verts[iv, 0], verts[iv, 1]) for iv in iverts]
+                    xc, yc = centroid_of_polygon(vlist)
+                    cell2d.append([icell, xc, yc, len(iverts)] + iverts)
+                    icell += 1
+
     kw = {}
     kw["nodes"] = nodes
     kw["nja"] = nja
-    kw["nvert"] = None
+    kw["nvert"] = nvert
     kw["top"] = top
     kw["bot"] = bot
     kw["area"] = area
@@ -181,6 +222,9 @@ def get_nn(k, i, j):
     kw["ihc"] = ihc
     kw["cl12"] = cl12
     kw["hwva"] = hwva
+    if return_vertices:
+        kw["vertices"] = vertices
+        kw["cell2d"] = cell2d
     return kw