fix(voronoi): clean up voronoi examples and add error check (#1323)

* Close #1299
* Voronoi issues in example were caused by closed polygons being passed to triangle
* Revised examples to meet requirements
* Add check for duplicate vertices in triangular mesh
* Updated t075 to check for invalid cells
* Changed several filtering loops in grid methods to use standard list comprehension

.github/workflows/ci.yml

@@ -82,7 +82,7 @@ jobs:
         echo "if black check fails run"
         echo "  black ./flopy"
         echo "and then commit the changes."
-        black --check ./flopy
+        black --check --diff ./flopy
 
     - name: Run flake8
       run: |

autotest/t075_test_ugrid.py

@@ -279,7 +279,7 @@ def test_voronoi_grid0(plot=False):
     test_setup = FlopyTestSetup(verbose=True, test_dirs=model_ws)
 
     name = "vor0"
-    answer_ncpl = 3804
+    answer_ncpl = 3803
     domain = [
         [1831.381546, 6335.543757],
         [4337.733475, 6851.136153],
@@ -298,7 +298,7 @@ def test_voronoi_grid0(plot=False):
         [1330.11116, 1809.788273],
         [399.1804436, 2998.515188],
         [914.7728404, 5132.494831],
-        [1831.381546, 6335.543757],
+#        [1831.381546, 6335.543757],
     ]
     area_max = 100.0 ** 2
     tri = Triangle(maximum_area=area_max, angle=30, model_ws=model_ws)
@@ -308,11 +308,6 @@ def test_voronoi_grid0(plot=False):
 
     vor = VoronoiGrid(tri)
     gridprops = vor.get_gridprops_vertexgrid()
-    ncpl = gridprops["ncpl"]
-    assert (
-        ncpl == answer_ncpl
-    ), f"Number of cells should be {answer_ncpl}. Found {ncpl}"
-
     voronoi_grid = VertexGrid(**gridprops, nlay=1)
 
     if plot:
@@ -324,6 +319,19 @@ def test_voronoi_grid0(plot=False):
         voronoi_grid.plot(ax=ax)
         plt.savefig(os.path.join(model_ws, f"{name}.png"))
 
+    # ensure proper number of cells
+    ncpl = gridprops["ncpl"]
+    errmsg = f"Number of cells should be {answer_ncpl}. Found {ncpl}"
+    assert ncpl == answer_ncpl, errmsg
+
+    # ensure that all cells have 3 or more points
+    ninvalid_cells = []
+    for icell, ivts in enumerate(vor.iverts):
+        if len(ivts) < 3:
+            ninvalid_cells.append(icell)
+    errmsg = f"The following cells do not have 3 or more vertices.\n{ninvalid_cells}"
+    assert len(ninvalid_cells) == 0, errmsg
+
     return
 
 
@@ -346,10 +354,6 @@ def test_voronoi_grid1(plot=False):
     vor = VoronoiGrid(tri)
     gridprops = vor.get_gridprops_vertexgrid()
     voronoi_grid = VertexGrid(**gridprops, nlay=1)
-    ncpl = gridprops["ncpl"]
-    assert (
-        ncpl == answer_ncpl
-    ), f"Number of cells should be {answer_ncpl}. Found {ncpl}"
 
     if plot:
         import matplotlib.pyplot as plt
@@ -360,6 +364,19 @@ def test_voronoi_grid1(plot=False):
         voronoi_grid.plot(ax=ax)
         plt.savefig(os.path.join(model_ws, f"{name}.png"))
 
+    # ensure proper number of cells
+    ncpl = gridprops["ncpl"]
+    errmsg = f"Number of cells should be {answer_ncpl}. Found {ncpl}"
+    assert ncpl == answer_ncpl, errmsg
+
+    # ensure that all cells have 3 or more points
+    ninvalid_cells = []
+    for icell, ivts in enumerate(vor.iverts):
+        if len(ivts) < 3:
+            ninvalid_cells.append(icell)
+    errmsg = f"The following cells do not have 3 or more vertices.\n{ninvalid_cells}"
+    assert len(ninvalid_cells) == 0, errmsg
+
     return
 
 
@@ -381,10 +398,6 @@ def test_voronoi_grid2(plot=False):
     vor = VoronoiGrid(tri)
     gridprops = vor.get_gridprops_vertexgrid()
     voronoi_grid = VertexGrid(**gridprops, nlay=1)
-    ncpl = gridprops["ncpl"]
-    assert (
-        ncpl == answer_ncpl
-    ), f"Number of cells should be {answer_ncpl}. Found {ncpl}"
 
     if plot:
         import matplotlib.pyplot as plt
@@ -395,6 +408,19 @@ def test_voronoi_grid2(plot=False):
         voronoi_grid.plot(ax=ax)
         plt.savefig(os.path.join(model_ws, f"{name}.png"))
 
+    # ensure proper number of cells
+    ncpl = gridprops["ncpl"]
+    errmsg = f"Number of cells should be {answer_ncpl}. Found {ncpl}"
+    assert ncpl == answer_ncpl, errmsg
+
+    # ensure that all cells have 3 or more points
+    ninvalid_cells = []
+    for icell, ivts in enumerate(vor.iverts):
+        if len(ivts) < 3:
+            ninvalid_cells.append(icell)
+    errmsg = f"The following cells do not have 3 or more vertices.\n{ninvalid_cells}"
+    assert len(ninvalid_cells) == 0, errmsg
+
     return
 
 
@@ -426,10 +452,6 @@ def test_voronoi_grid3(plot=False):
     vor = VoronoiGrid(tri)
     gridprops = vor.get_gridprops_vertexgrid()
     voronoi_grid = VertexGrid(**gridprops, nlay=1)
-    ncpl = gridprops["ncpl"]
-    assert (
-        ncpl == answer_ncpl
-    ), f"Number of cells should be {answer_ncpl}. Found {ncpl}"
 
     if plot:
         import matplotlib.pyplot as plt
@@ -440,6 +462,19 @@ def test_voronoi_grid3(plot=False):
         voronoi_grid.plot(ax=ax)
         plt.savefig(os.path.join(model_ws, f"{name}.png"))
 
+    # ensure proper number of cells
+    ncpl = gridprops["ncpl"]
+    errmsg = f"Number of cells should be {answer_ncpl}. Found {ncpl}"
+    assert ncpl == answer_ncpl, errmsg
+
+    # ensure that all cells have 3 or more points
+    ninvalid_cells = []
+    for icell, ivts in enumerate(vor.iverts):
+        if len(ivts) < 3:
+            ninvalid_cells.append(icell)
+    errmsg = f"The following cells do not have 3 or more vertices.\n{ninvalid_cells}"
+    assert len(ninvalid_cells) == 0, errmsg
+
     return
 
 
@@ -464,10 +499,6 @@ def test_voronoi_grid4(plot=False):
     vor = VoronoiGrid(tri)
     gridprops = vor.get_gridprops_vertexgrid()
     voronoi_grid = VertexGrid(**gridprops, nlay=1)
-    ncpl = gridprops["ncpl"]
-    assert (
-        ncpl == answer_ncpl
-    ), f"Number of cells should be {answer_ncpl}. Found {ncpl}"
 
     if plot:
         import matplotlib.pyplot as plt
@@ -478,6 +509,19 @@ def test_voronoi_grid4(plot=False):
         voronoi_grid.plot(ax=ax)
         plt.savefig(os.path.join(model_ws, f"{name}.png"))
 
+    # ensure proper number of cells
+    ncpl = gridprops["ncpl"]
+    errmsg = f"Number of cells should be {answer_ncpl}. Found {ncpl}"
+    assert ncpl == answer_ncpl, errmsg
+
+    # ensure that all cells have 3 or more points
+    ninvalid_cells = []
+    for icell, ivts in enumerate(vor.iverts):
+        if len(ivts) < 3:
+            ninvalid_cells.append(icell)
+    errmsg = f"The following cells do not have 3 or more vertices.\n{ninvalid_cells}"
+    assert len(ninvalid_cells) == 0, errmsg
+
     return
 
 
@@ -486,10 +530,10 @@ def test_voronoi_grid5(plot=False):
     test_setup = FlopyTestSetup(verbose=True, test_dirs=model_ws)
 
     name = "vor5"
-    answer_ncpl = 1067
+    answer_ncpl = 1305
     active_domain = [(0, 0), (100, 0), (100, 100), (0, 100)]
     area1 = [(10, 10), (40, 10), (40, 40), (10, 40)]
-    area2 = [(50, 50), (90, 50), (90, 90), (50, 90)]
+    area2 = [(70, 70), (90, 70), (90, 90), (70, 90)]
 
     tri = Triangle(angle=30, model_ws=model_ws)
 
@@ -515,25 +559,26 @@ def test_voronoi_grid5(plot=False):
     # tri.add_hole((70, 20))
 
     # add line through domain to force conforming cells
-    line = [(x, x) for x in np.linspace(10, 90, 100)]
+    line = [(x, x) for x in np.linspace(11, 89, 100)]
     tri.add_polygon(line)
 
     # then regions and other polygons should follow
     tri.add_polygon(area1)
     tri.add_polygon(area2)
     tri.add_region((1, 1), 0, maximum_area=100)  # point inside active domain
     tri.add_region((11, 11), 1, maximum_area=10)  # point inside area1
-    tri.add_region((70, 70), 2, maximum_area=3)  # point inside area2
+    tri.add_region((70, 70), 2, maximum_area=1)  # point inside area2
 
     tri.build(verbose=False)
 
     vor = VoronoiGrid(tri)
     gridprops = vor.get_gridprops_vertexgrid()
     voronoi_grid = VertexGrid(**gridprops, nlay=1)
-    ncpl = gridprops["ncpl"]
-    assert (
-        ncpl == answer_ncpl
-    ), f"Number of cells should be {answer_ncpl}. Found {ncpl}"
+
+    ninvalid_cells = []
+    for icell, ivts in enumerate(vor.iverts):
+        if len(ivts) < 3:
+            ninvalid_cells.append(icell)
 
     if plot:
         import matplotlib.pyplot as plt
@@ -542,8 +587,21 @@ def test_voronoi_grid5(plot=False):
         ax = fig.add_subplot()
         ax.set_aspect("equal")
         voronoi_grid.plot(ax=ax)
+
+        # plot invalid cells
+        ax.plot(voronoi_grid.xcellcenters[ninvalid_cells], voronoi_grid.ycellcenters[ninvalid_cells], 'ro')
+
         plt.savefig(os.path.join(model_ws, f"{name}.png"))
 
+    # ensure proper number of cells
+    ncpl = gridprops["ncpl"]
+    errmsg = f"Number of cells should be {answer_ncpl}. Found {ncpl}"
+    assert ncpl == answer_ncpl, errmsg
+
+    # ensure that all cells have 3 or more points
+    errmsg = f"The following cells do not have 3 or more vertices.\n{ninvalid_cells}"
+    assert len(ninvalid_cells) == 0, errmsg
+
     return
 
 
@@ -556,9 +614,9 @@ def test_voronoi_grid5(plot=False):
     # test_create_unstructured_grid_from_verts()
     # test_triangle_unstructured_grid()
     # test_voronoi_vertex_grid()
-    test_voronoi_grid0(plot=True)
-    test_voronoi_grid1(plot=True)
-    test_voronoi_grid2(plot=True)
-    test_voronoi_grid3(plot=True)
-    test_voronoi_grid4(plot=True)
+    #test_voronoi_grid0(plot=True)
+    #test_voronoi_grid1(plot=True)
+    #test_voronoi_grid2(plot=True)
+    #test_voronoi_grid3(plot=True)
+    #test_voronoi_grid4(plot=True)
     test_voronoi_grid5(plot=True)

flopy/discretization/grid.py

@@ -353,6 +353,10 @@ def lenuni(self):
     def idomain(self):
         return copy.deepcopy(self._idomain)
 
+    @idomain.setter
+    def idomain(self, idomain):
+        self._idomain = idomain
+
     @property
     def ncpl(self):
         raise NotImplementedError("must define ncpl in child class")

flopy/discretization/unstructuredgrid.py

@@ -201,7 +201,9 @@ def nvert(self):
     @property
     def iverts(self):
         if self._iverts is not None:
-            return [list(filter((None).__ne__, i)) for i in self._iverts]
+            return [
+                [ivt for ivt in t if ivt is not None] for t in self._iverts
+            ]
 
     @property
     def verts(self):

flopy/discretization/vertexgrid.py

@@ -133,12 +133,16 @@ def iverts(self):
     @property
     def cell1d(self):
         if self._cell1d is not None:
-            return [list(filter((None).__ne__, t)) for t in self._cell1d]
+            return [
+                [ivt for ivt in t if ivt is not None] for t in self._cell2d
+            ]
 
     @property
     def cell2d(self):
         if self._cell2d is not None:
-            return [list(filter((None).__ne__, t)) for t in self._cell2d]
+            return [
+                [ivt for ivt in t if ivt is not None] for t in self._cell2d
+            ]
 
     @property
     def verts(self):

flopy/utils/voronoi.py

@@ -82,16 +82,27 @@ def tri2vor(tri, **kwargs):
     tri_verts = tri.verts
     tri_iverts = tri.iverts
     tri_edge = tri.edge
+    npoints = tri_verts.shape[0]
+    ntriangles = len(tri_iverts)
+    nedges = tri_edge.shape[0]
+
+    # check to make sure there are no duplicate points
+    tri_verts_unique = np.unique(tri_verts, axis=0)
+    if tri_verts.shape != tri_verts_unique.shape:
+        npoints_unique = tri_verts_unique.shape[0]
+        errmsg = (
+            f"There are duplicate points in the triangular mesh. "
+            f"These can be caused by overlapping regions, holes, and "
+            f"refinement features.  The triangular mesh has {npoints} "
+            f"points but only {npoints_unique} are unique."
+        )
+        raise Exception(errmsg)
 
     # construct the voronoi grid
     vor = Voronoi(tri_verts, **kwargs)
     ridge_points = vor.ridge_points
     ridge_vertices = vor.ridge_vertices
 
-    npoints = tri_verts.shape[0]
-    ntriangles = len(tri_iverts)
-    nedges = tri_edge.shape[0]
-
     # test the voronoi vertices, and mark those outside of the domain
     nvertices = vor.vertices.shape[0]
     xc = vor.vertices[:, 0].reshape((nvertices, 1))
@@ -115,6 +126,9 @@ def tri2vor(tri, **kwargs):
     # renumber valid vertices consecutively
     idx_vertindex[idx_filtered] = np.arange(nvalid_vertices)
 
+    # Create new lists for the voronoi grid vertices and the
+    # voronoi grid incidence list.  There should be one voronoi
+    # cell for each vertex point in the triangular grid
     vor_verts = [(x, y) for x, y in vor.vertices[idx_filtered]]
     vor_iverts = [[] for i in range(npoints)]
 
@@ -181,9 +195,8 @@ def tri2vor(tri, **kwargs):
     if True:
         vor_verts = np.array(vor_verts)
         for icell in range(len(vor_iverts)):
-            vor_iverts[icell] = get_sorted_vertices(
-                vor_iverts[icell], vor_verts
-            )
+            iverts_cell = vor_iverts[icell]
+            vor_iverts[icell] = get_sorted_vertices(iverts_cell, vor_verts)
 
     return vor_verts, vor_iverts