From ad4e817d3feb76fe72873f54adbf383aaf649bfe Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Tue, 23 Dec 2025 08:33:37 +0100 Subject: [PATCH] Fixing bug in uxgrid index_search This makes sure shapes stay consistent for single-face queries, and adds a small safeguard against zero norms --- src/parcels/_core/index_search.py | 14 ++++++++------ src/parcels/_core/uxgrid.py | 2 +- 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/src/parcels/_core/index_search.py b/src/parcels/_core/index_search.py index a225621fc..c2dd18351 100644 --- a/src/parcels/_core/index_search.py +++ b/src/parcels/_core/index_search.py @@ -239,18 +239,20 @@ def uxgrid_point_in_cell(grid, y: np.ndarray, x: np.ndarray, yi: np.ndarray, xi: axis=-1, ) - # Get projection points onto element plane - # for the projection, all points are computed relative to v0 - r1 = np.squeeze(face_vertices[:, 1, :] - face_vertices[:, 0, :]) # (M,3) - r2 = np.squeeze(face_vertices[:, 2, :] - face_vertices[:, 0, :]) # (M,3) + # Get projection points onto element plane. Keep the leading + # dimension even for single-face queries so shapes remain (M,3). + r1 = face_vertices[:, 1, :] - face_vertices[:, 0, :] + r2 = face_vertices[:, 2, :] - face_vertices[:, 0, :] nhat = np.cross(r1, r2) norm = np.linalg.norm(nhat, axis=-1) + # Avoid division by zero for degenerate faces + norm = np.where(norm == 0.0, 1.0, norm) nhat = nhat / norm[:, None] # Calculate the component of the points in the direction of nhat - ptilde = points - np.squeeze(face_vertices[:, 0, :]) + ptilde = points - face_vertices[:, 0, :] pdotnhat = np.sum(ptilde * nhat, axis=-1) # Reconstruct points with normal component removed. - points = ptilde - pdotnhat[:, None] * nhat + np.squeeze(face_vertices[:, 0, :]) + points = ptilde - pdotnhat[:, None] * nhat + face_vertices[:, 0, :] else: nids = grid.uxgrid.face_node_connectivity[xi].values diff --git a/src/parcels/_core/uxgrid.py b/src/parcels/_core/uxgrid.py index 5e2f628a8..a1d45e796 100644 --- a/src/parcels/_core/uxgrid.py +++ b/src/parcels/_core/uxgrid.py @@ -103,7 +103,7 @@ def search(self, z, y, x, ei=None, tol=1e-6): if np.any(ei): indices = self.unravel_index(ei) fi = indices.get("FACE") - is_in_cell, coords = uxgrid_point_in_cell(self.uxgrid, y, x, fi, fi) + is_in_cell, coords = uxgrid_point_in_cell(self, y, x, fi, fi) y_check = y[is_in_cell == 0] x_check = x[is_in_cell == 0] zero_indices = np.where(is_in_cell == 0)[0]