feat(cvfdutil): capability to create disv nested grids (#997)

The modflow6 disv constructor in flopy needs a list of vertices and cell2d information. The flopy.util.cvfdutil package was updated to include a new method that takes a list of flopy.discretization.structuredgrid.StructuredGrid objects and stitches them together to create a dictionary of arguments that can be passed info flopy.modflow6.ModflowGwfdisv. For nested grids, idomain must be set to zero in any parent grids in order for the child grid to be patched in properly. This new tool can be used to construct vertices and cell2d information needed to make disv nested grids.

autotest/t073_test_cvfd.py

@@ -0,0 +1,78 @@
+import numpy as np
+import flopy
+from flopy.utils.cvfdutil import to_cvfd, gridlist_to_disv_gridprops
+
+
+def test_tocvfd1():
+    vertdict = {}
+    vertdict[0] = [(0, 0), (100, 0), (100, 100), (0, 100), (0, 0)]
+    vertdict[1] = [(100, 0), (120, 0), (120, 20), (100, 20), (100, 0)]
+    verts, iverts = to_cvfd(vertdict)
+    assert 6 in iverts[0]
+
+
+def test_tocvfd2():
+    vertdict = {}
+    vertdict[0] = [(0, 0), (1, 0), (1, 1), (0, 1), (0, 0)]
+    vertdict[1] = [(1, 0), (3, 0), (3, 2), (1, 2), (1, 0)]
+    verts, iverts = to_cvfd(vertdict)
+    assert [1, 4, 5, 6, 2, 1] in iverts
+
+
+def test_tocvfd3():
+    # create the nested grid described in the modflow-usg documentation
+
+    # outer grid
+    nlay = 1
+    nrow = ncol = 7
+    delr = 100.0 * np.ones(ncol)
+    delc = 100.0 * np.ones(nrow)
+    tp = np.zeros((nrow, ncol))
+    bt = -100.0 * np.ones((nlay, nrow, ncol))
+    idomain = np.ones((nlay, nrow, ncol))
+    idomain[:, 2:5, 2:5] = 0
+    sg1 = flopy.discretization.StructuredGrid(
+        delr=delr, delc=delc, top=tp, botm=bt, idomain=idomain
+    )
+    # inner grid
+    nlay = 1
+    nrow = ncol = 9
+    delr = 100.0 / 3.0 * np.ones(ncol)
+    delc = 100.0 / 3.0 * np.ones(nrow)
+    tp = np.zeros((nrow, ncol))
+    bt = -100 * np.ones((nlay, nrow, ncol))
+    idomain = np.ones((nlay, nrow, ncol))
+    sg2 = flopy.discretization.StructuredGrid(
+        delr=delr,
+        delc=delc,
+        top=tp,
+        botm=bt,
+        xoff=200.0,
+        yoff=200,
+        idomain=idomain,
+    )
+    gridprops = gridlist_to_disv_gridprops([sg1, sg2])
+    assert "ncpl" in gridprops
+    assert "nvert" in gridprops
+    assert "vertices" in gridprops
+    assert "cell2d" in gridprops
+
+    ncpl = gridprops["ncpl"]
+    nvert = gridprops["nvert"]
+    vertices = gridprops["vertices"]
+    cell2d = gridprops["cell2d"]
+    assert ncpl == 121
+    assert nvert == 148
+    assert len(vertices) == nvert
+    assert len(cell2d) == 121
+
+    # spot check information for cell 28 (zero based)
+    answer = [28, 250.0, 150.0, 7, 38, 142, 143, 45, 46, 44, 38]
+    for i, j in zip(cell2d[28], answer):
+        assert i == j, "{} not equal {}".format(i, j)
+
+
+if __name__ == "__main__":
+    test_tocvfd1()
+    test_tocvfd2()
+    test_tocvfd3()

flopy/utils/cvfdutil.py

@@ -70,32 +70,61 @@ def shared_face(ivlist1, ivlist2):
 
 
 def segment_face(ivert, ivlist1, ivlist2, vertices):
-    ic1pos = ivlist1.index(ivert)
-    if ic1pos == 0:  # if ivert is first, then must also be last
-        ic1pos = len(ivlist1) - 1
-    ic1v1 = ivlist1[ic1pos - 1]
-    ic1v2 = ivlist1[ic1pos]
-
-    x, y = vertices[ic1v1]
-    a = Point(x, y)
+    """
+    Check the vertex lists for cell 1 and cell 2.  Add a new vertex to cell 1
+    if necessary.
 
-    x, y = vertices[ic1v2]
-    b = Point(x, y)
+    Parameters
+    ----------
+    ivert : int
+        vertex number to check
+    ivlist1 : list
+        list of vertices for cell 1.  Add a new vertex to this cell if needed.
+    ivlist2 : list
+        list of vertices for cell2.
+    vertices : ndarray
+        array of x, y vertices
 
-    ic2pos = ivlist2.index(ivert)
-    ic2v2 = ivlist2[ic2pos + 1]
-    x, y = vertices[ic2v2]
-    c = Point(x, y)
+    Returns
+    -------
+    segmented : bool
+        Return True if a face in cell 1 was split up by adding a new vertex
 
-    if ic1v1 == ic2v2 or ic1v2 == ic2v2:
-        return
+    """
 
-    # print('Checking segment {} {} with point {}'.format(ic1v1, ic1v2, ic2v2))
+    # go through ivlist1 and find faces that have ivert
+    faces_to_check = []
+    for ipos in range(len(ivlist1) - 1):
+        face = (ivlist1[ipos], ivlist1[ipos + 1])
+        if ivert in face:
+            faces_to_check.append(face)
+
+    # go through ivlist2 and find points to check
+    points_to_check = []
+    for ipos in range(len(ivlist2) - 1):
+        if ivlist2[ipos] == ivert:
+            points_to_check.append(ivlist2[ipos + 1])
+        elif ivlist2[ipos + 1] == ivert:
+            points_to_check.append(ivlist2[ipos])
+
+    for face in faces_to_check:
+        iva, ivb = face
+        x, y = vertices[iva]
+        a = Point(x, y)
+        x, y = vertices[ivb]
+        b = Point(x, y)
+        for ivc in points_to_check:
+            if ivc not in face:
+                x, y = vertices[ivc]
+                c = Point(x, y)
+                if isBetween(a, b, c):
+                    ipos = ivlist1.index(ivb)
+                    if ipos == 0:
+                        ipos = len(ivlist1) - 1
+                    ivlist1.insert(ipos, ivc)
+                    return True
 
-    if isBetween(a, b, c):
-        # print('between: ', a, b, c)
-        ivlist1.insert(ic1pos, ic2v2)
-    return
+    return False
 
 
 def to_cvfd(
@@ -106,7 +135,7 @@ def to_cvfd(
     verbose=False,
 ):
     """
-    Convert a vertex dictionary
+    Convert a vertex dictionary into verts and iverts
 
     Parameters
     ----------
@@ -193,35 +222,45 @@ def to_cvfd(
         ivertlist = vertexlist[icell]
         for ivert in ivertlist:
             if ivert in vertex_cell_dict:
-                vertex_cell_dict[ivert].append(icell)
+                if icell not in vertex_cell_dict[ivert]:
+                    vertex_cell_dict[ivert].append(icell)
             else:
                 vertex_cell_dict[ivert] = [icell]
+    if verbose:
+        print("Done creating dict of vertices with their associated cells")
 
     # Now, go through each vertex and look at the cells that use the vertex.
     # For quadtree-like grids, there may be a need to add a new hanging node
     # vertex to the larger cell.
-    if verbose:
-        print("Done creating dict of vertices with their associated cells")
-        print("Checking for hanging nodes.")
-    vertexdict_keys = list(vertexdict.keys())
-    for ivert, cell_list in vertex_cell_dict.items():
-        for icell1 in cell_list:
-            for icell2 in cell_list:
-
-                # skip if same cell
-                if icell1 == icell2:
-                    continue
-
-                # skip if share face already
-                ivertlist1 = vertexlist[icell1]
-                ivertlist2 = vertexlist[icell2]
-                if shared_face(ivertlist1, ivertlist2):
-                    continue
-
-                # don't share a face, so need to segment if necessary
-                segment_face(ivert, ivertlist1, ivertlist2, vertexdict_keys)
-    if verbose:
-        print("Done checking for hanging nodes.")
+    if not skip_hanging_node_check:
+        if verbose:
+            print("Checking for hanging nodes.")
+        vertexdict_keys = list(vertexdict.keys())
+        finished = False
+        while not finished:
+            finished = True
+            for ivert, cell_list in vertex_cell_dict.items():
+                for icell1 in cell_list:
+                    for icell2 in cell_list:
+
+                        # skip if same cell
+                        if icell1 == icell2:
+                            continue
+
+                        # skip if share face already
+                        ivertlist1 = vertexlist[icell1]
+                        ivertlist2 = vertexlist[icell2]
+                        if shared_face(ivertlist1, ivertlist2):
+                            continue
+
+                        # don't share a face, so need to segment if necessary
+                        segmented = segment_face(
+                            ivert, ivertlist1, ivertlist2, vertexdict_keys
+                        )
+                        if segmented:
+                            finished = False
+        if verbose:
+            print("Done checking for hanging nodes.")
 
     verts = np.array(vertexdict_keys)
     iverts = vertexlist
@@ -272,3 +311,109 @@ def shapefile_to_xcyc(shp):
         xcyc[icell, 0] = xc
         xcyc[icell, 1] = yc
     return xcyc
+
+
+def gridlist_to_verts(gridlist):
+    """
+
+    Take a list of flopy structured model grids and convert them into vertices.
+    The idomain can be set to remove cells in a parent grid.  Cells from a
+    child grid will patched in to make a single set of vertices.  Cells will
+    be numbered according to consecutive numbering of active cells in the
+    grid list.
+
+    Parameters
+    ----------
+    gridlist : list
+        List of flopy.discretization.modelgrid.  Must be of type structured
+        grids
+
+    Returns
+    -------
+    verts, iverts : np.ndarray, list
+        vertices and list of cells and which vertices comprise the cells
+    """
+    vertdict = {}
+    icell = 0
+    for sg in gridlist:
+        ilays, irows, icols = np.where(sg.idomain > 0)
+        for _, i, j in zip(ilays, irows, icols):
+            v = sg.get_cell_vertices(i, j)
+            vertdict[icell] = v + [v[0]]
+            icell += 1
+    verts, iverts = to_cvfd(vertdict, verbose=False)
+    return verts, iverts
+
+
+def get_disv_gridprops(verts, iverts):
+    """
+
+    Take a list of flopy structured model grids and convert them into vertices.
+    The idomain can be set to remove cells in a parent grid.  Cells from a
+    child grid will patched in to make a single set of vertices.  Cells will
+    be numbered according to consecutive numbering of active cells in the
+    grid list.
+
+    Parameters
+    ----------
+    gridlist : list
+        List of flopy.discretization.modelgrid.  Must be of type structured
+        grids
+
+    Returns
+    -------
+    gridprops : dict
+        Dictionary containing entries that can be passed directly into the
+        modflow6 disv package.
+
+    """
+    nvert = verts.shape[0]
+    ncpl = len(iverts)
+    xcyc = np.empty((ncpl, 2), dtype=np.float)
+    area = np.empty(ncpl, dtype=np.float)
+    for icell in range(ncpl):
+        vlist = [(verts[ivert, 0], verts[ivert, 1]) for ivert in iverts[icell]]
+        xcyc[icell, 0], xcyc[icell, 1] = centroid_of_polygon(vlist)
+        area[icell] = abs(area_of_polygon(*zip(*vlist)))
+    vertices = []
+    for i in range(nvert):
+        vertices.append((i, verts[i, 0], verts[i, 1]))
+    cell2d = []
+    for i in range(ncpl):
+        cell2d.append(
+            [i, xcyc[i, 0], xcyc[i, 1], len(iverts[i])]
+            + [iv for iv in iverts[i]]
+        )
+    gridprops = {}
+    gridprops["ncpl"] = ncpl
+    gridprops["nvert"] = nvert
+    gridprops["vertices"] = vertices
+    gridprops["cell2d"] = cell2d
+    return gridprops
+
+
+def gridlist_to_disv_gridprops(gridlist):
+    """
+
+    Take a list of flopy structured model grids and convert them into a
+    dictionary that can be passed into the modflow6 disv package.  Cells from a
+    child grid will patched in to make a single set of vertices.  Cells will
+    be numbered according to consecutive numbering of active cells in the
+    grid list.
+
+    Parameters
+    ----------
+    gridlist : list
+        List of flopy.discretization.modelgrid.  Must be of type structured
+        grids
+
+    Returns
+    -------
+    gridprops : dict
+        Dictionary containing entries that can be passed directly into the
+        modflow6 disv package.
+
+    """
+    verts, iverts = gridlist_to_verts(gridlist)
+    gridprops = get_disv_gridprops(verts, iverts)
+    return gridprops