feat(grid): add ihc to UnstructuredGrid

flopy/discretization/unstructuredgrid.py

@@ -73,6 +73,11 @@ class UnstructuredGrid(Grid):
         optional number of connections per node array
     ja : list or ndarray
         optional jagged connection array
+    ihc : list or ndarray
+        optional horizontal connection indicator array (for MODFLOW 6 DISU).
+        For each connection in the ja array: ihc = 0 indicates a vertical
+        connection, ihc = 1 or 2 indicates a horizontal connection, with 2
+        indicating that horizontal connections are vertically staggered.
     **kwargs : dict, optional
         Support deprecated keyword options.
 
@@ -132,6 +137,7 @@ def __init__(
         angrot=0.0,
         iac=None,
         ja=None,
+        ihc=None,
         cell2d=None,
         **kwargs,
     ):
@@ -187,6 +193,7 @@ def __init__(
 
         self._iac = iac
         self._ja = ja
+        self._ihc = ihc
 
     def set_ncpl(self, ncpl):
         if isinstance(ncpl, int):
@@ -299,6 +306,10 @@ def iac(self):
     def ja(self):
         return self._ja
 
+    @property
+    def ihc(self):
+        return self._ihc
+
     @property
     def ncpl(self):
         return self._ncpl
@@ -656,6 +667,9 @@ def convert_grid(self, factor):
                 xoff=self.xoffset * factor,
                 yoff=self.yoffset * factor,
                 angrot=self.angrot,
+                iac=self._iac,
+                ja=self._ja,
+                ihc=self._ihc,
             )
         else:
             raise AssertionError("Grid is not complete and cannot be converted")
@@ -716,6 +730,7 @@ def clean_iverts(self, inplace=False):
                     angrot=self.angrot,
                     iac=self._iac,
                     ja=self._ja,
+                    ihc=self._ihc,
                 )
 
     def intersect(self, x, y, z=None, local=False, forgive=False):

flopy/mf6/mfmodel.py

@@ -562,6 +562,7 @@ def modelgrid(self):
                 angrot=self._modelgrid.angrot,
                 iac=dis.iac.array,
                 ja=dis.ja.array,
+                ihc=dis.ihc.array,
             )
         elif self.get_grid_type() == DiscretizationType.DISV1D:
             dis = self.get_package("disv1d")

flopy/plot/plotutil.py

@@ -3064,6 +3064,11 @@ def is_vertical_barrier(mg, cellid1, cellid2) -> bool:
     """
     Determine if a barrier is vertical (between vertically stacked cells).
 
+    For structured (DIS) and vertex (DISV) grids, uses cellid structure to
+    determine orientation. For unstructured (DISU) grids, prefers to use the
+    ihc (horizontal connection indicator) array if available, falling back to
+    geometric analysis of shared face z-coordinates.
+
     Parameters
     ----------
     mg : Grid
@@ -3092,14 +3097,28 @@ def is_vertical_barrier(mg, cellid1, cellid2) -> bool:
         return cellid1[0] != cellid2[0] and cellid1[1] == cellid2[1]
     else:
         # Unstructured grid: (node,)
-        # Infer from geometry: check the orientation of the shared face
-        # If the face is horizontal (all z-coords equal), it's a vertical barrier
-        # If the face is vertical (z-coords differ), it's a horizontal barrier
-        shared_face_3d = get_shared_face_3d(mg, cellid1, cellid2)
-        if shared_face_3d is None:
-            # No shared face found, can't determine orientation
+        if mg.grid_type != "unstructured":
+            raise ValueError("Given 1-element node number, expected unstructured grid")
+
+        # Try to use connectivity data if available
+        if mg.ihc is not None and mg.iac is not None:
+            node1 = cellid1[0]
+            node2 = cellid2[0]
+            if node1 < len(mg.iac) - 1:
+                idx0 = mg.iac[node1]
+                idx1 = mg.iac[node1 + 1]
+                ja = mg.ja
+                ihc = mg.ihc
+                for i in range(idx0 + 1, idx1):
+                    if ja[i] == node2:
+                        # ihc == 0 means vertical, ihc != 0 means horizontal
+                        return ihc[i] == 0
+
+        # Fall back to geometric approach if ihc not available. In this case we
+        # look at the cell geometries. Check the orientation of the shared face.
+        # If the face is horizontal (z coordinates equal), then it's a vertical
+        # barrier. If the face is vertical (z coords vary), horizontal barrier.
+        if (shared_face_3d := get_shared_face_3d(mg, cellid1, cellid2)) is None:
             return False
-
-        # Check if all z-coordinates are the same (horizontal face)
         z_coords = [v[2] for v in shared_face_3d]
         return np.allclose(z_coords, z_coords[0], rtol=1e-5)