catch duplicates

autotest/test_gridgen.py

@@ -4,6 +4,7 @@
 
 import matplotlib.pyplot as plt
 import numpy as np
+import pandas as pd
 import pytest
 from matplotlib.collections import LineCollection, PathCollection, QuadMesh
 from modflow_devtools.markers import requires_exe, requires_pkg
@@ -853,10 +854,14 @@ def test_gridgen(function_tmpdir):
     assert max(ja0) <= disu_vp.nodelay[0], msg
 
 
-@pytest.mark.xfail
 @requires_exe("mf6", "gridgen")
 def test_flopy_issue_1492(function_tmpdir):
-    name = "mymodel"
+    """
+    Submitted by David Brakenhoff in
+    https://github.com/modflowpy/flopy/issues/1492
+    """
+
+    name = "issue1492"
     nlay = 3
     nrow = 10
     ncol = 10
@@ -871,7 +876,6 @@ def test_flopy_issue_1492(function_tmpdir):
         sim_name=name, sim_ws=function_tmpdir, exe_name="mf6"
     )
     gwf = flopy.mf6.ModflowGwf(sim, modelname=name)
-
     dis = flopy.mf6.ModflowGwfdis(
         gwf,
         nlay=nlay,
@@ -882,6 +886,7 @@ def test_flopy_issue_1492(function_tmpdir):
         top=top,
         botm=botm,
     )
+    og_grid = gwf.modelgrid
 
     # Create and build the gridgen model with a refined area in the middle
     g = Gridgen(dis, model_ws=function_tmpdir)
@@ -925,17 +930,30 @@ def test_flopy_issue_1492(function_tmpdir):
         head_filerecord=head_file,
         saverecord=[("HEAD", "ALL"), ("BUDGET", "ALL")],
     )
+
     sim.write_simulation()
     success, _ = sim.run_simulation(silent=False)
     assert success
 
     head = gwf.output.head().get_data()
     bud = gwf.output.budget()
     spdis = bud.get_data(text="DATA-SPDIS")[0]
-
-    pmv = flopy.plot.PlotMapView(gwf)
-    pmv.plot_array(head)
-    pmv.plot_grid(colors="white")
-    pmv.contour_array(head, levels=[0.2, 0.4, 0.6, 0.8], linewidths=3.0)
-    pmv.plot_vector(spdis["qx"], spdis["qy"], color="white")
-    plt.show()
+    grid = gwf.modelgrid
+    og_verts = pd.DataFrame(
+        og_grid.verts, columns=["x", "y"]
+    )  # .round(3).sort_values(by=["x", "y"], ignore_index=True).drop_duplicates(ignore_index=True)
+    mg_verts = pd.DataFrame(
+        grid.verts, columns=["x", "y"]
+    )  # .round(3).sort_values(by=["x", "y"], ignore_index=True).drop_duplicates(ignore_index=True)
+
+    plot_debug = False
+    if plot_debug:
+        pmv = flopy.plot.PlotMapView(gwf)
+        pmv.plot_array(head)
+        pmv.plot_grid(colors="white")
+        ax = plt.gca()
+        verts = grid.verts
+        ax.plot(verts[:, 0], verts[:, 1], "bo", alpha=0.25, ms=5)
+        pmv.contour_array(head, levels=[0.2, 0.4, 0.6, 0.8], linewidths=3.0)
+        pmv.plot_vector(spdis["qx"], spdis["qy"], color="white")
+        plt.show()

flopy/utils/cvfdutil.py

@@ -1,4 +1,5 @@
 import numpy as np
+import pandas as pd
 
 from .utl_import import import_optional_dependency
 
@@ -176,7 +177,7 @@ def to_cvfd(
         xcyc[icell, 1] = yc
         ivertlist = []
         for p in points:
-            pt = tuple(p)
+            pt = (round(p[0], 9), round(p[1], 9))
             if pt in vertexdict:
                 ivert = vertexdict[pt]
             else:
@@ -240,10 +241,21 @@ def to_cvfd(
         if verbose:
             print("Done checking for hanging nodes.")
 
-    verts = np.array(vertexdict_keys)
-    iverts = vertexlist
-
-    return verts, iverts
+    # drop duplicate vertices
+    verts = (
+        pd.DataFrame(np.array(vertexdict_keys), columns=["x", "y"])
+        .round(9)
+        .drop_duplicates(["x", "y"], ignore_index=False)
+        .to_numpy()
+    )
+
+    def get_iverts():
+        vtups = [tuple(v) for v in verts.tolist()]
+        vdict = {k: v for k, v in vertexdict.items() if k in vtups}
+        for v in vertexlist:
+            yield [vv for vv in v if vv in vdict.values()]
+
+    return verts, list(get_iverts())
 
 
 def shapefile_to_cvfd(shp, **kwargs):

flopy/utils/gridgen.py

@@ -1891,8 +1891,9 @@ def _mkvertdict(self):
         idx = attributes.index("nodenumber")
         for i in range(len(shapes)):
             nodenumber = int(records[i][idx]) - 1
-            self._vertdict[nodenumber] = shapes[i].points
-        return
+            points = shapes[i].points
+            # import pdb; pdb.set_trace()
+            self._vertdict[nodenumber] = points
 
     @staticmethod
     def read_qtg_nod(