feat(cvfdutil): capability to create 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