diff --git a/xugrid/conversion.py b/xugrid/conversion.py index ee42cdda..d6540373 100644 --- a/xugrid/conversion.py +++ b/xugrid/conversion.py @@ -283,23 +283,32 @@ def bounds1d_to_vertices(bounds: np.ndarray): def bounds2d_to_vertices( x_bounds: np.ndarray, y_bounds: np.ndarray ) -> Tuple[FloatArray, BoolArray]: - # Create 2D coords, discard NaNs. + x = x_bounds.reshape(-1, 4) + y = y_bounds.reshape(-1, 4) + # Make sure repeated nodes are consecutive so we can find them later. + # lexsort along axis 1. + sorter = np.lexsort((y, x)) face_node_coordinates = np.stack( ( - x_bounds.reshape(-1, 4), - y_bounds.reshape(-1, 4), + np.take_along_axis(x, sorter, axis=1), + np.take_along_axis(y, sorter, axis=1), ), axis=-1, ) index = np.isfinite(face_node_coordinates.reshape(-1, 8)).all(axis=-1) face_node_coordinates = face_node_coordinates[index] - # Guarantee counterclockwise orientation + # Guarantee counterclockwise orientation. face_centroids = np.mean(face_node_coordinates, axis=1) dx = face_node_coordinates[..., 0] - face_centroids[:, np.newaxis, 0] dy = face_node_coordinates[..., 1] - face_centroids[:, np.newaxis, 1] angle = np.arctan2(dy, dx) + # When nodes are repeated, make sure the repeated node ends up as last so we can + # construct the face_node_connectivity directly from it. We do so by inserting + # an angle of np.inf for the repeated nodes. + angle[:, 1:][angle[:, 1:] == angle[:, :-1]] = np.inf counterclockwise = np.argsort(angle, axis=1) + face_node_coordinates = np.take_along_axis( face_node_coordinates, counterclockwise[..., None], axis=1 ) diff --git a/xugrid/ugrid/ugrid2d.py b/xugrid/ugrid/ugrid2d.py index 6691d88f..bf2b586a 100644 --- a/xugrid/ugrid/ugrid2d.py +++ b/xugrid/ugrid/ugrid2d.py @@ -2251,10 +2251,29 @@ def from_structured_bounds( xy, inverse = np.unique( face_node_coordinates.reshape((-1, 2)), return_inverse=True, axis=0 ) + face_node_connectivity = inverse.reshape((-1, 4)) + + # Check whether all coordinates form valid UGRID topologies. + # We can only maintain triangles and quadrangles. + valid = face_node_connectivity != np.roll(face_node_connectivity, 1, axis=1) + if not valid.all(): + xy = xy[np.unique(face_node_connectivity[valid])] + keep = valid.sum(axis=1) >= 3 + warnings.warn( + "A UGRID2D face requires at least three unique vertices.\n" + f"Your structured bounds contain {len(keep) - keep.sum()} invalid faces.\n" + "These will be omitted from the Ugrid2d topology.", + ) + index = index[keep] + face_node_connectivity[~valid] = -1 + face_node_connectivity = connectivity.renumber( + face_node_connectivity[keep] + ) + grid = Ugrid2d( *xy.T, -1, - face_node_connectivity=inverse.reshape((-1, 4)), + face_node_connectivity=face_node_connectivity, name=name, ) else: