handle irregular yz cross sections

CHANGES.rst

@@ -6,6 +6,7 @@ Changelog of threedigrid-builder
 -------------------
 
 - Define new boundary condition types.
+- Handle open YZ profiles correctly when the middle is higher than the sides.
 
 
 1.19.0 (2024-09-10)

threedigrid_builder/grid/cross_section_definitions.py

@@ -319,15 +319,17 @@ def tabulate_yz(shape, width, height):
         seen.add(x)
         zs[i] = x
 
-    # For open profiles, if left and right sides are not equal in Z, we envision a
-    # vertical line from the left or right most point depending on which is lower.
-    # This is to ensure that the width is non-decreasing for the remaining profile.
+    # For open profiles where the left or right sides are not the maximum Z,
+    # we add a vertical line at each end such that the left and right sides are
+    # now the maximum Z. This ensures that the profile is not inverted by
+    # shapely when calculating the widths.
     if not is_closed:
-        if zs[0] < zs[-1]:
-            zs = np.concatenate(([zs[-1]], zs))
+        highest_z = np.max(zs)
+        if zs[0] < highest_z:
+            zs = np.concatenate(([highest_z], zs))
             ys = np.concatenate(([ys[0]], ys))
-        if zs[0] > zs[-1]:
-            zs = np.concatenate((zs, [zs[0]]))
+        if zs[-1] < highest_z:
+            zs = np.concatenate((zs, [highest_z]))
             ys = np.concatenate((ys, [ys[-1]]))
 
     # shapely will automatically close an open profile

threedigrid_builder/tests/test_cross_section_definitions.py

@@ -381,6 +381,32 @@ def test_tabulate_inverted_egg():
                 )
             ),
         ),
+        (  # Open profile, left side rises then falls
+            "0 5 10 15 20 25 30 35 40",
+            "5.2 5.26 5.25 5.2 5.1 5.1 0 5 5",
+            None,
+            None,
+            None,
+            None,
+            None,
+            40.0,
+            5.26,
+            np.column_stack(
+                (
+                    [0.0, 5.0, 5.0, 5.1, 5.1, 5.2, 5.2, 5.25, 5.26],
+                    [0.0, 9.902, 14.902, 15.0, 20.005, 24.995, 25.008, 34.167, 40.0],
+                )
+            ),
+            np.column_stack(
+                (
+                    [0.0, 5.0, 10.0, 15.0, 20.0, 25.0, 30.0, 35.0, 40.0],
+                    [5.2, 5.26, 5.25, 5.2, 5.1, 5.1, 0.0, 5.0, 5.0],
+                    np.zeros(9, dtype=float),
+                    np.zeros(9, dtype=float),
+                    np.zeros(9, dtype=float),
+                )
+            ),
+        ),
     ],
 )
 def test_tabulate_yz(