diff --git a/autotest/test_util_geometry.py b/autotest/test_util_geometry.py
index 79135f326..2c70c730f 100644
--- a/autotest/test_util_geometry.py
+++ b/autotest/test_util_geometry.py
@@ -30,11 +30,60 @@ def test_does_isclockwise_work():
     assert bool(rslt) is False, "is_clockwise() failed"
 
 
-def test_point_in_polygon_basic():
+def debug_plot(grid, cell, xpts, ypts, mask):
+    import matplotlib.pyplot as plt
+
+    grid.plot()
+    plt.plot(*list(zip(*cell)))
+    plt.scatter(xpts[mask], ypts[mask], c="red")
+    # plt.show()
+
+
+def test_point_in_polygon_interior():
     grid = GridCases().structured_small()
+    cell = grid.verts[grid.iverts[0]].tolist()
     xpts = grid.xcellcenters
     ypts = grid.ycellcenters
-    poly = grid.verts[grid.iverts[0]].tolist()
-    mask = point_in_polygon(xpts, ypts, poly)
+    mask = point_in_polygon(xpts, ypts, cell)
     assert mask.sum() == 1
     assert mask[0, 0] == True
+    debug_plot(grid, cell, xpts, ypts, mask)
+
+
+def test_point_in_polygon_vertices():
+    grid = GridCases().structured_small()
+    cell = grid.verts[grid.iverts[0]].tolist()
+    xpts, ypts = list(zip(*cell))
+    xpts = np.array([xpts])
+    ypts = np.array([ypts])
+    mask = point_in_polygon(xpts, ypts, cell)
+    assert mask.sum() == 1  # only bottom left (due to axis-alignment)
+    debug_plot(grid, cell, xpts, ypts, mask)
+
+    # move points inside boundary
+    xpts[0, 2] = xpts[0, 2] - 0.001
+    ypts[0, 2] = ypts[0, 2] + 0.001
+    xpts[0, 1] = xpts[0, 1] - 0.001
+    ypts[0, 1] = ypts[0, 1] - 0.001
+    xpts[0, 0] = xpts[0, 0] + 0.001
+    ypts[0, 0] = ypts[0, 0] - 0.001
+    mask = point_in_polygon(xpts, ypts, cell)
+    assert mask.sum() == 4  # expect all points
+    debug_plot(grid, cell, xpts, ypts, mask)
+
+
+def test_point_in_polygon_faces():
+    grid = GridCases().structured_small()
+    cell = grid.verts[grid.iverts[0]].tolist()
+    xpts_v, ypts_v = list(zip(*cell))
+    xpts_v = np.array([xpts_v])
+    ypts_v = np.array([ypts_v])
+    xpts = np.array(
+        [[xpts_v[0, 0], xpts_v[0, 2], np.mean(xpts_v), np.mean(xpts_v)]]
+    )
+    ypts = np.array(
+        [[np.mean(ypts_v), np.mean(ypts_v), ypts_v[0, 0], ypts_v[0, 2]]]
+    )
+    mask = point_in_polygon(xpts, ypts, cell)
+    assert mask.sum() == 2  # only inner faces
+    debug_plot(grid, cell, xpts, ypts, mask)