refactor(plot): better shared face finding for HFB plotting (#2682)

Improve the shared face finding algorithm for HFB (Horizontal Flow Barrier) plotting. Factor out @jlarsen-usgs's index-based approach from #2671 into reusable function(s) and use it instead of floating point comparisons as before. Should be faster.

I moved the new face-related functions to a new module faceutil.py instead of plotutil.py as they are used not only for plotting but for export, and may be useful in other cases. I considered putting them in gridutil.py but that is used by the grid classes, while the face functions use the grid classes, so it would have created a circular import.

I also changed the name of is_vertical_barrier -> is_vertical and inverted the semantics so it works on faces, not barriers (a face between horizontally adjacent cells is a vertical polygon, and a face between vertically adjacent cells is a horizontal polygon, but we call the former a "horizontal"  flow barrier and the latter a "vertical" flow barrier). I figure writing these utils in terms of faces for the general case makes sense, and just interpreting the results for the specific case of HFBs.

This PR builds on

- #2671's shared face finding approach
- #2680 which gives all grid types a consistent get_node() method
- #2681 which adds an ihc init parameter/attribute to UnstructuredGrid

Some ad hoc cell ID -> node number conversion utilities are removed and get_node() used instead.

I am classifying all this a refactor though it introduces a new utility module, as the new utils are in service of the refactor.

Author: wpbonelli

autotest/test_faceutil.py

@@ -0,0 +1,182 @@
+import numpy as np
+import pytest
+
+from flopy.discretization import StructuredGrid, UnstructuredGrid, VertexGrid
+from flopy.utils.faceutil import get_shared_face_3d
+
+
+def test_get_shared_face_3d_structured_horizontal():
+    nlay, nrow, ncol = 1, 2, 2
+    delr = np.array([1.0, 1.0])
+    delc = np.array([1.0, 1.0])
+    top = np.array([[10.0, 10.0], [10.0, 10.0]])
+    botm = np.array([[[0.0, 0.0], [0.0, 0.0]]])
+
+    grid = StructuredGrid(delc=delc, delr=delr, top=top, botm=botm)
+
+    # Test horizontally adjacent cells (0, 0, 0) and (0, 0, 1)
+    cellid1 = (0, 0, 0)
+    cellid2 = (0, 0, 1)
+
+    face = get_shared_face_3d(grid, cellid1, cellid2)
+
+    assert face is not None
+    assert len(face) == 4  # Vertical face should have 4 vertices
+
+    # Check that face is vertical (two z-coordinates: top and bottom)
+    z_coords = sorted({v[2] for v in face})
+    assert len(z_coords) == 2
+    assert z_coords[0] == 0.0  # bottom
+    assert z_coords[1] == 10.0  # top
+
+    # Check x-coordinate is at the shared boundary (x=1.0)
+    x_coords = {v[0] for v in face}
+    assert 1.0 in x_coords
+
+
+def test_get_shared_face_3d_structured_vertical():
+    nlay, nrow, ncol = 2, 2, 2
+    delr = np.array([1.0, 1.0])
+    delc = np.array([1.0, 1.0])
+    top = np.array([[10.0, 10.0], [10.0, 10.0]])
+    botm = np.array([[[5.0, 5.0], [5.0, 5.0]], [[0.0, 0.0], [0.0, 0.0]]])
+
+    grid = StructuredGrid(delc=delc, delr=delr, top=top, botm=botm)
+
+    # Test vertically adjacent cells (0, 0, 0) and (1, 0, 0)
+    cellid1 = (0, 0, 0)
+    cellid2 = (1, 0, 0)
+
+    face = get_shared_face_3d(grid, cellid1, cellid2)
+
+    assert face is not None
+    assert len(face) == 2  # Horizontal face should have 2 vertices (the edge)
+
+    # Check that face is horizontal (all z-coordinates the same at interface)
+    z_coords = {v[2] for v in face}
+    assert len(z_coords) == 1
+    assert next(iter(z_coords)) == 5.0  # Interface between layers
+
+
+def test_get_shared_face_3d_vertex():
+    nlay = 2
+    ncpl = 2
+
+    # Two square cells side by side
+    vertices = [
+        [0, 0.0, 0.0],
+        [1, 1.0, 0.0],
+        [2, 2.0, 0.0],
+        [3, 0.0, 1.0],
+        [4, 1.0, 1.0],
+        [5, 2.0, 1.0],
+    ]
+
+    cell2d = [
+        [0, 0.5, 0.5, 4, 0, 1, 4, 3],
+        [1, 1.5, 0.5, 4, 1, 2, 5, 4],
+    ]
+
+    top = np.array([10.0, 10.0])
+    botm = np.array([[5.0, 5.0], [0.0, 0.0]])
+
+    grid = VertexGrid(vertices=vertices, cell2d=cell2d, top=top, botm=botm, nlay=nlay)
+
+    # Test horizontally adjacent cells (0, 0) and (0, 1)
+    cellid1 = (0, 0)
+    cellid2 = (0, 1)
+
+    face = get_shared_face_3d(grid, cellid1, cellid2)
+
+    assert face is not None
+    assert len(face) == 4  # Vertical face
+
+    # Test vertically adjacent cells (0, 0) and (1, 0)
+    cellid1 = (0, 0)
+    cellid2 = (1, 0)
+
+    face = get_shared_face_3d(grid, cellid1, cellid2)
+
+    assert face is not None
+    assert len(face) == 2  # Horizontal face
+
+    # Check all z-coordinates are the same (interface elevation)
+    z_coords = {v[2] for v in face}
+    assert len(z_coords) == 1
+    assert next(iter(z_coords)) == 5.0
+
+
+def test_get_shared_face_3d_unstructured():
+    # Two cells stacked vertically
+    vertices = [
+        # Bottom layer vertices (z=0)
+        [0, 0.0, 0.0, 0.0],
+        [1, 1.0, 0.0, 0.0],
+        [2, 1.0, 1.0, 0.0],
+        [3, 0.0, 1.0, 0.0],
+        # Middle layer vertices (z=5)
+        [4, 0.0, 0.0, 5.0],
+        [5, 1.0, 0.0, 5.0],
+        [6, 1.0, 1.0, 5.0],
+        [7, 0.0, 1.0, 5.0],
+        # Top layer vertices (z=10)
+        [8, 0.0, 0.0, 10.0],
+        [9, 1.0, 0.0, 10.0],
+        [10, 1.0, 1.0, 10.0],
+        [11, 0.0, 1.0, 10.0],
+    ]
+
+    iverts = [
+        [0, 1, 2, 3, 4, 5, 6, 7],  # Bottom cell
+        [4, 5, 6, 7, 8, 9, 10, 11],  # Top cell
+    ]
+
+    xc = np.array([0.5, 0.5])
+    yc = np.array([0.5, 0.5])
+    top = np.array([10.0, 10.0])
+    botm = np.array([0.0, 5.0])
+
+    grid = UnstructuredGrid(
+        vertices=vertices,
+        iverts=iverts,
+        xcenters=xc,
+        ycenters=yc,
+        top=top,
+        botm=botm,
+        ncpl=[2],
+    )
+
+    # Test shared face between two vertically stacked cells
+    cellid1 = (0,)
+    cellid2 = (1,)
+
+    face = get_shared_face_3d(grid, cellid1, cellid2)
+
+    assert face is not None
+    assert len(face) == 4  # Should have 4 shared vertices at the interface
+
+    # Check all shared vertices are at z=5.0 (the interface)
+    z_coords = {v[2] for v in face}
+    assert len(z_coords) == 1
+    assert next(iter(z_coords)) == 5.0
+
+
+def test_get_shared_face_3d_error_cases():
+    nlay, nrow, ncol = 1, 2, 2
+    delr = np.array([1.0, 1.0])
+    delc = np.array([1.0, 1.0])
+    top = np.array([[10.0, 10.0], [10.0, 10.0]])
+    botm = np.array([[[0.0, 0.0], [0.0, 0.0]]])
+
+    grid = StructuredGrid(delc=delc, delr=delr, top=top, botm=botm)
+
+    # Test with same cellid
+    with pytest.raises(ValueError, match="cellid1 and cellid2 must be different"):
+        get_shared_face_3d(grid, (0, 0, 0), (0, 0, 0))
+
+    # Test with non-adjacent cells (should return None)
+    cellid1 = (0, 0, 0)
+    cellid2 = (0, 1, 1)  # Diagonally opposite, not adjacent
+
+    face = get_shared_face_3d(grid, cellid1, cellid2)
+    assert face is None

flopy/mf6/mfpackage.py

@@ -2160,7 +2160,7 @@ def to_geodataframe(self, gdf=None, kper=0, full_grid=True, shorten_attr=False,
                 if modelgrid is not None:
                     if self.package_type == "hfb":
                         gpd = import_optional_dependency("geopandas")
-                        from ..plot.plotutil import hfb_data_to_linework
+                        from ..utils.faceutil import hfb_data_to_linework
 
                         recarray = self.stress_period_data.data[kper]
                         lines = hfb_data_to_linework(recarray, modelgrid)

flopy/modflow/mfhfb.py

@@ -190,7 +190,7 @@ def _get_hfb_lines(self):
         -------
             list : list of hfb lines
         """
-        from ..plot.plotutil import hfb_data_to_linework
+        from ..utils.faceutil import hfb_data_to_linework
 
         return hfb_data_to_linework(self.hfb_data, self.parent.modelgrid)
 

flopy/plot/crosssection.py

@@ -13,7 +13,7 @@
 from . import plotutil
 from .plotutil import (
     get_shared_face_3d,
-    is_vertical_barrier,
+    is_vertical,
     to_mp7_endpoints,
     to_mp7_pathlines,
 )
@@ -781,22 +781,6 @@ def plot_grid(self, **kwargs):
             ax.add_collection(col)
         return col
 
-    def _cellid_to_node(self, cellid):
-        """
-        Convert a cellid tuple to a node number.
-
-        Parameters
-        ----------
-        cellid : tuple
-            Cell identifier
-
-        Returns
-        -------
-        int
-            Node number
-        """
-        return self.mg.get_node([cellid])[0]
-
     def _plot_vertical_hfb_lines(self, color=None, **kwargs):
         """
         Plot vertical HFBs as lines at layer interfaces.
@@ -848,8 +832,8 @@ def _plot_vertical_hfb_lines(self, color=None, **kwargs):
 
             # Get the x-coordinates along the cross section for this cell
             # Need to find this cell in projpts - convert both cellids to nodes
-            node1 = self._cellid_to_node(cellid1)
-            node2 = self._cellid_to_node(cellid2)
+            node1 = self.mg.get_node([cellid1])[0]
+            node2 = self.mg.get_node([cellid2])[0]
 
             # Get x-coordinates from either node's projection
             xcoords = []
@@ -976,11 +960,11 @@ def plot_bc(
                             cellid1 = tuple(entry["cellid1"])
                             cellid2 = tuple(entry["cellid2"])
 
-                            if is_vertical_barrier(self.mg, cellid1, cellid2):
-                                # Store vertical HFBs for line plotting
+                            if not is_vertical(self.mg, cellid1, cellid2):
+                                # horizontal face, vertical barrier
                                 vertical_hfbs.append((cellid1, cellid2))
                             else:
-                                # Horizontal barriers - plot both cells as patches
+                                # vertical face, horizontal barrier
                                 cellids.append(list(entry["cellid1"]))
                                 cellids.append(list(entry["cellid2"]))
 

flopy/plot/map.py

@@ -13,7 +13,7 @@
 from .plotutil import (
     get_shared_face,
     get_shared_face_3d,
-    is_vertical_barrier,
+    is_vertical,
     to_mp7_pathlines,
 )
 
@@ -469,20 +469,18 @@ def _plot_barrier_bc(self, barrier_data, color=None, name=None, **kwargs):
         """
         ax = kwargs.pop("ax", self.ax)
 
-        # Separate horizontal and vertical barriers
         horizontal_line_segments = []
         vertical_cell_indices = []
 
         for cellid1, cellid2 in barrier_data:
             # Only plot barriers on the current layer (for layered grids)
             # For DISU (len==1), plot all barriers since there's no layer filtering
-            if len(cellid1) >= 2:  # DIS or DISV (has layer info)
+            if len(cellid1) >= 2:
                 if cellid1[0] != self.layer and cellid2[0] != self.layer:
                     continue
 
-            # Check if this is a vertical barrier
-            if is_vertical_barrier(self.mg, cellid1, cellid2):
-                # For vertical barriers, plot the cells on the current layer
+            # horizontal face, vertical barrier
+            if not is_vertical(self.mg, cellid1, cellid2):
                 if len(cellid1) >= 2:
                     if cellid1[0] == self.layer:
                         if len(cellid1) == 3:
@@ -499,7 +497,7 @@ def _plot_barrier_bc(self, barrier_data, color=None, name=None, **kwargs):
                         else:
                             vertical_cell_indices.append([self.layer, cellid2[1]])
             else:
-                # Horizontal barrier - plot as line on shared face
+                # vertical face, horizontal barrier
                 shared_face = get_shared_face(self.mg, cellid1, cellid2)
                 if shared_face is not None:
                     horizontal_line_segments.append(shared_face)