feat(unstructured): improve unstructured grid support for MODFLOW 6 a…

…nd MODFLOW-USG (#1021)

* Revise and extend the VertexGrid and UnstructuredGrid objects for plotting and export to shapefiles; these are the objects that are accessed through gwf.modelgrid, for example
* Support layered plotting for DISU models when the grid is different for each layer
* For DISU grids, if the plottable layer number is assigned to the diagonal position of the IHC array, then these numbers will be used as layers for layered plots, such as those shown by PlotMapView()
* Revise the gridgen wrapper to pass the plotabble layer number through the IHC array so that layered plotting makes sense
* model.plot(), package.plot(), and array.plot() functionality now works for DIS, DISV, and DISU grids
* Update the flopy3_grigen.ipynb notebook to demonstrate use of Gridgen to construct mf6 DISV and DISU simulations and mfusg disu simulation of a simple flow problem (the same one shown on the main flopy repo page
* Add additional tests for the new capabilities
* More work is needed for vtk and netcdf output; better vtk support for unstructured grids is being worked on separately
* Cross sections are not yet supported for unstructured grids
* Add support for MODFLOW-USG plotting, especially through gridgen
* May help with some MODFLOW-USG issues(#896, #1010, #1024) but more work is needed
* Closes #1018

autotest/t061_test_gridgen.py

@@ -122,8 +122,8 @@ def test_gridgen():
         # test the different gridprops dictionaries, which contain all the
         # information needed to make the different discretization packages
         gridprops = g.get_gridprops_disv()
-        gridprops = g.get_gridprops()
-        #gridprops = g.get_gridprops_disu6()
+        gridprops = g.get_gridprops_disu5()
+        gridprops = g.get_gridprops_disu6()
 
         # test the gridgen point intersection
         points = [(4750., 5250.)]
@@ -147,18 +147,6 @@ def test_gridgen():
         mu = flopy.modflow.Modflow(version='mfusg', structured=False)
         disu = g.get_disu(mu)
 
-        # test writing a modflow 6 disu package
-        fname = os.path.join(cpth, 'mymf6model.disu')
-        g6.to_disu6(fname)
-        assert os.path.isfile(fname), \
-            'MF6 disu file not created: {}'.format(fname)
-
-        # test writing a modflow 6 disv package
-        fname = os.path.join(cpth, 'mymf6model.disv')
-        g6.to_disv6(fname)
-        assert os.path.isfile(fname), \
-            'MF6 disv file not created: {}'.format(fname)
-
         # test mfusg with vertical pass-through (True above at instantiation)
         gu.build()
         disu_vp = gu.get_disu(ms_u)

autotest/t075_ugridtests.py

@@ -0,0 +1,245 @@
+"""
+Unstructured grid tests
+
+"""
+
+import os
+import numpy as np
+from flopy.discretization import UnstructuredGrid
+from flopy.utils.triangle import Triangle
+
+tpth = os.path.join("temp", "t075")
+if not os.path.isdir(tpth):
+    os.mkdir(tpth)
+
+
+def test_unstructured_grid_shell():
+    # constructor with no arguments.  incomplete shell should exist
+    g = UnstructuredGrid()
+    assert g.nlay is None
+    assert g.nnodes is None
+    assert g.ncpl is None
+    assert not g.grid_varies_by_layer
+    assert not g.is_valid
+    assert not g.is_complete
+    return
+
+
+def test_unstructured_grid_dimensions():
+    # constructor with just dimensions
+    ncpl = [1, 10, 1]
+    g = UnstructuredGrid(ncpl=ncpl)
+    assert np.allclose(g.ncpl, ncpl)
+    assert g.nlay == 3
+    assert g.nnodes == 12
+    assert not g.is_valid
+    assert not g.is_complete
+    assert not g.grid_varies_by_layer
+    return
+
+
+def test_unstructured_minimal_grid():
+
+    # pass in simple 2 cell minimal grid to make grid valid
+    vertices = [
+        [0, 0.0, 1.0],
+        [1, 1.0, 1.0],
+        [2, 2.0, 1.0],
+        [3, 0.0, 0.0],
+        [4, 1.0, 0.0],
+        [5, 2.0, 0.0],
+    ]
+    iverts = [[0, 1, 4, 3], [1, 2, 5, 4]]
+    xcenters = [0.5, 1.5]
+    ycenters = [0.5, 0.5]
+    g = UnstructuredGrid(
+        vertices=vertices, iverts=iverts, xcenters=xcenters, ycenters=ycenters
+    )
+    assert np.allclose(g.ncpl, np.array([2], dtype=np.int))
+    assert g.nlay == 1
+    assert g.nnodes == 2
+    assert g.is_valid
+    assert not g.is_complete
+    assert not g.grid_varies_by_layer
+    assert g._vertices == vertices
+    assert g._iverts == iverts
+    assert g._xc == xcenters
+    assert g._yc == ycenters
+    grid_lines = [
+        [(0.0, 0), (0.0, 1.0)],
+        [(0.0, 1), (1.0, 1.0)],
+        [(1.0, 1), (1.0, 0.0)],
+        [(1.0, 0), (0.0, 0.0)],
+        [(1.0, 0), (1.0, 1.0)],
+        [(1.0, 1), (2.0, 1.0)],
+        [(2.0, 1), (2.0, 0.0)],
+        [(2.0, 0), (1.0, 0.0)],
+    ]
+    assert g.grid_lines == grid_lines, "\n{} \n /=   \n{}".format(
+        g.grid_lines, grid_lines
+    )
+    assert g.extent == (0, 2, 0, 1)
+    xv, yv, zv = g.xyzvertices
+    assert xv == [[0, 1, 1, 0], [1, 2, 2, 1]]
+    assert yv == [[1, 1, 0, 0], [1, 1, 0, 0]]
+    assert zv is None
+
+    return
+
+
+def test_unstructured_complete_grid():
+
+    # pass in simple 2 cell complete grid to make grid valid, and put each
+    # cell in a different layer
+    vertices = [
+        [0, 0.0, 1.0],
+        [1, 1.0, 1.0],
+        [2, 2.0, 1.0],
+        [3, 0.0, 0.0],
+        [4, 1.0, 0.0],
+        [5, 2.0, 0.0],
+    ]
+    iverts = [[0, 1, 4, 3], [1, 2, 5, 4]]
+    xcenters = [0.5, 1.5]
+    ycenters = [0.5, 0.5]
+    ncpl = [1, 1]
+    top = [1, 0]
+    top = np.array(top)
+    botm = [0, -1]
+    botm = np.array(botm)
+    g = UnstructuredGrid(
+        vertices=vertices,
+        iverts=iverts,
+        xcenters=xcenters,
+        ycenters=ycenters,
+        ncpl=ncpl,
+        top=top,
+        botm=botm,
+    )
+    assert np.allclose(g.ncpl, np.array([1, 1], dtype=np.int))
+    assert g.nlay == 2
+    assert g.nnodes == 2
+    assert g.is_valid
+    assert not g.is_complete
+    assert g.grid_varies_by_layer
+    assert g._vertices == vertices
+    assert g._iverts == iverts
+    assert g._xc == xcenters
+    assert g._yc == ycenters
+    grid_lines = {
+        0: [
+            [(0.0, 0.0), (0.0, 1.0)],
+            [(0.0, 1.0), (1.0, 1.0)],
+            [(1.0, 1.0), (1.0, 0.0)],
+            [(1.0, 0.0), (0.0, 0.0)],
+        ],
+        1: [
+            [(1.0, 0.0), (1.0, 1.0)],
+            [(1.0, 1.0), (2.0, 1.0)],
+            [(2.0, 1.0), (2.0, 0.0)],
+            [(2.0, 0.0), (1.0, 0.0)],
+        ],
+    }
+    assert isinstance(g.grid_lines, dict)
+    assert g.grid_lines == grid_lines, "\n{} \n /=   \n{}".format(
+        g.grid_lines, grid_lines
+    )
+    assert g.extent == (0, 2, 0, 1)
+    xv, yv, zv = g.xyzvertices
+    assert xv == [[0, 1, 1, 0], [1, 2, 2, 1]]
+    assert yv == [[1, 1, 0, 0], [1, 1, 0, 0]]
+    assert np.allclose(zv, np.array([[1, 0], [0, -1]]))
+
+    return
+
+
+def test_loading_argus_meshes():
+    datapth = os.path.join("..", "examples", "data", "unstructured")
+    fnames = [fname for fname in os.listdir(datapth) if fname.endswith(".exp")]
+    for fname in fnames:
+        fname = os.path.join(datapth, fname)
+        print("Loading Argus mesh ({}) into UnstructuredGrid".format(fname))
+        g = UnstructuredGrid.from_argus_export(fname)
+        print("  Number of nodes: {}".format(g.nnodes))
+
+
+def test_create_unstructured_grid_from_verts():
+
+    datapth = os.path.join("..", "examples", "data", "unstructured")
+
+    # simple functions to load vertices and incidence lists
+    def load_verts(fname):
+        print("Loading vertices from: {}".format(fname))
+        verts = np.genfromtxt(
+            fname, dtype=[int, float, float], names=["iv", "x", "y"]
+        )
+        verts["iv"] -= 1  # zero based
+        return verts
+
+    def load_iverts(fname):
+        print("Loading iverts from: {}".format(fname))
+        f = open(fname, "r")
+        iverts = []
+        xc = []
+        yc = []
+        for line in f:
+            ll = line.strip().split()
+            iverts.append([int(i) - 1 for i in ll[4:]])
+            xc.append(float(ll[1]))
+            yc.append(float(ll[2]))
+        return iverts, np.array(xc), np.array(yc)
+
+    # load vertices
+    fname = os.path.join(datapth, "ugrid_verts.dat")
+    verts = load_verts(fname)
+
+    # load the incidence list into iverts
+    fname = os.path.join(datapth, "ugrid_iverts.dat")
+    iverts, xc, yc = load_iverts(fname)
+
+    ncpl = np.array(5 * [len(iverts)])
+    g = UnstructuredGrid(verts, iverts, xc, yc, ncpl=ncpl)
+    assert isinstance(g.grid_lines, list)
+    assert np.allclose(g.ncpl, ncpl)
+    assert g.extent == (0.0, 700.0, 0.0, 700.0)
+    assert g._vertices.shape == (156,)
+    assert g.nnodes == g.ncpl.sum() == 1090
+    return
+
+
+def test_triangle_unstructured_grid():
+    maximum_area = 30000.0
+    extent = (214270.0, 221720.0, 4366610.0, 4373510.0)
+    domainpoly = [
+        (extent[0], extent[2]),
+        (extent[1], extent[2]),
+        (extent[1], extent[3]),
+        (extent[0], extent[3]),
+    ]
+    tri = Triangle(maximum_area=maximum_area, angle=30, model_ws=tpth)
+    tri.add_polygon(domainpoly)
+    tri.build(verbose=False)
+    verts = [[iv, x, y] for iv, (x, y) in enumerate(tri.verts)]
+    iverts = tri.iverts
+    xc, yc = tri.get_xcyc().T
+    ncpl = np.array([len(iverts)])
+    g = UnstructuredGrid(
+        vertices=verts,
+        iverts=iverts,
+        ncpl=ncpl,
+        xcenters=xc,
+        ycenters=yc,
+    )
+    assert len(g.grid_lines) == 8190
+    assert g.nnodes == g.ncpl == 2730
+    return
+
+
+if __name__ == "__main__":
+    test_unstructured_grid_shell()
+    test_unstructured_grid_dimensions()
+    test_unstructured_minimal_grid()
+    test_unstructured_complete_grid()
+    test_loading_argus_meshes()
+    test_create_unstructured_grid_from_verts()
+    test_triangle_unstructured_grid()