feat(voronoi): add VoronoiGrid class (#1034)

Add a VoronoiGrid class for building MODFLOW models with voronoi grids.  The class requires a set of points as input.  These points can be generated using Triangle, for example.

autotest/t075_ugridtests.py

@@ -5,8 +5,9 @@
 
 import os
 import numpy as np
-from flopy.discretization import UnstructuredGrid
+from flopy.discretization import UnstructuredGrid, VertexGrid
 from flopy.utils.triangle import Triangle
+from flopy.utils.voronoi import VoronoiGrid
 
 tpth = os.path.join("temp", "t075")
 if not os.path.isdir(tpth):
@@ -235,6 +236,32 @@ def test_triangle_unstructured_grid():
     return
 
 
+def test_voronoi_vertex_grid():
+    xmin = 0.0
+    xmax = 2.0
+    ymin = 0.0
+    ymax = 1.0
+    area_max = 0.05
+    tri = Triangle(maximum_area=area_max, angle=30, model_ws=tpth)
+    poly = np.array(((xmin, ymin), (xmax, ymin), (xmax, ymax), (xmin, ymax)))
+    tri.add_polygon(poly)
+    tri.build(verbose=False)
+
+    # create vor object and VertexGrid
+    vor = VoronoiGrid(tri.verts)
+    gridprops = vor.get_gridprops_vertexgrid()
+    vgrid = VertexGrid(**gridprops, nlay=1)
+    assert vgrid.is_valid
+
+    # arguments for creating a mf6 disv package
+    gridprops = vor.get_disv_gridprops()
+    assert gridprops["ncpl"] == 43
+    assert gridprops["nvert"] == 83
+    assert len(gridprops["vertices"]) == 83
+    assert len(gridprops["cell2d"]) == 43
+    return
+
+
 if __name__ == "__main__":
     test_unstructured_grid_shell()
     test_unstructured_grid_dimensions()
@@ -243,3 +270,4 @@ def test_triangle_unstructured_grid():
     test_loading_argus_meshes()
     test_create_unstructured_grid_from_verts()
     test_triangle_unstructured_grid()
+    test_voronoi_vertex_grid()

flopy/discretization/unstructuredgrid.py

@@ -58,15 +58,17 @@ class UnstructuredGrid(Grid):
     This class handles spatial representation of unstructured grids.  It is
     based on the concept of being able to support multiple model layers that
     may have a different number of cells in each layer.  The array ncpl is of
-    size nlay and and its sum must equal nodes.  If the grid_varies_by_layer
-    flag is set to true, then the iverts array must be of size ncpl[0] and all
-    values in the ncpl array must be equal to nodes / nlay.  The xcenters and
-    ycenters arrays must also be of size ncpl[0].  This makes it
-    possible to efficiently store spatial grid information for multiple layers.
+    size nlay and and its sum must equal nodes.  If the length of iverts is
+    equal to ncpl[0] and the number of cells per layer is the same for each
+    layer, then it is assumed that the grid does not vary by layer.  In this
+    case, the xcenters and ycenters arrays must also be of size ncpl[0].
+    This makes it possible to efficiently store spatial grid information
+    for multiple layers.
 
     If the spatial grid is different for each model layer, then the
-    grid_varies_by_layer flag should be false, and iverts must be of size
-    nodes. The arrays for xcenters and ycenters must also be of size nodes.
+    grid_varies_by_layer flag will automatically be set to false, and iverts
+    must be of size nodes. The arrays for xcenters and ycenters must also
+    be of size nodes.
 
     """
 

flopy/utils/cvfdutil.py

@@ -345,20 +345,18 @@ def gridlist_to_verts(gridlist):
     return verts, iverts
 
 
-def get_disv_gridprops(verts, iverts):
+def get_disv_gridprops(verts, iverts, xcyc=None):
     """
 
-    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.
+    Calculate disv grid properties from verts and iverts
 
     Parameters
     ----------
-    gridlist : list
-        List of flopy.discretization.modelgrid.  Must be of type structured
-        grids
+    verts : ndarray
+        2d array of x, y vertices
+
+    iverts : list
+        list of size ncpl, with a list of vertex numbers for each cell
 
     Returns
     -------
@@ -369,12 +367,15 @@ def get_disv_gridprops(verts, iverts):
     """
     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)))
+    if xcyc is None:
+        xcyc = np.empty((ncpl, 2), 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)
+    else:
+        assert xcyc.shape == (ncpl, 2)
     vertices = []
     for i in range(nvert):
         vertices.append((i, verts[i, 0], verts[i, 1]))

flopy/utils/gridgen.py

@@ -1450,7 +1450,7 @@ def get_gridprops_disu5(self):
 
     def get_gridprops_disu6(self, repair_asymmetry=True):
         """
-        Get a dictionary of information needed to create a MODFLOW-USG DISU
+        Get a dictionary of information needed to create a MODFLOW 6 DISU
         Package.  The returned dictionary can be unpacked directly into the
         ModflowGwfdisu constructor.
 

flopy/utils/triangle.py

@@ -26,6 +26,10 @@ class Triangle(object):
     angle : float
         Triangle will continue to add vertices until no angle is less than
         this specified value.  (default is 20 degrees)
+    nodes : ndarray
+        Two dimensional array of shape (npoints, 2) with x and y positions
+        of fixed node locations to include in the resulting triangular mesh.
+        (default is None)
     additional_args : list
         list of additional command line switches to pass to triangle
 
@@ -41,6 +45,7 @@ def __init__(
         exe_name="triangle",
         maximum_area=None,
         angle=20.0,
+        nodes=None,
         additional_args=None,
     ):
         self.model_ws = model_ws
@@ -50,6 +55,7 @@ def __init__(
         self.exe_name = os.path.abspath(exe_name)
         self.angle = angle
         self.maximum_area = maximum_area
+        self._nodes = nodes
         self.additional_args = additional_args
         self._initialize_vars()
         return
@@ -741,6 +747,8 @@ def _write_nodefile(self, fname):
         nvert = 0
         for p in self._polygons:
             nvert += len(p)
+        if self._nodes is not None:
+            nvert += self._nodes.shape[0]
         s = "{} {} {} {}\n".format(nvert, 2, 0, 0)
         f.write(s)
         ip = 0
@@ -749,6 +757,13 @@ def _write_nodefile(self, fname):
                 s = "{} {} {}\n".format(ip, vertex[0], vertex[1])
                 f.write(s)
                 ip += 1
+        if self._nodes is not None:
+            for i in range(self._nodes.shape[0]):
+                s = "{} {} {}\n".format(
+                    ip, self._nodes[i, 0], self._nodes[i, 1]
+                )
+                f.write(s)
+                ip += 1
         f.close()
 
     def _write_polyfile(self, fname):