From 6611d12737bdf3021542d4bddcbc6a49dc36357c Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Tue, 20 Aug 2024 22:35:00 +0200
Subject: [PATCH 1/8] Add an option to skip concave exterior vertices

---
 xugrid/ugrid/connectivity.py | 18 ++++++++----
 xugrid/ugrid/voronoi.py      | 55 ++++++++++++++++++++++++++++++++++--
 2 files changed, 65 insertions(+), 8 deletions(-)

diff --git a/xugrid/ugrid/connectivity.py b/xugrid/ugrid/connectivity.py
index 42e8db28b..82d66c38b 100644
--- a/xugrid/ugrid/connectivity.py
+++ b/xugrid/ugrid/connectivity.py
@@ -571,20 +571,26 @@ def perimeter(
     return np.linalg.norm(dxy, axis=-1).sum(axis=1)
 
 
+def area_from_coordinates(
+    coordinates: FloatArray,
+) -> FloatArray:
+    xy0 = coordinates[:, 0]
+    a = coordinates[:, :-1] - xy0[:, np.newaxis]
+    b = coordinates[:, 1:] - xy0[:, np.newaxis]
+    determinant = np.cross(a, b)
+    return 0.5 * abs(determinant.sum(axis=1))
+
+
 def area(
     face_node_connectivity: IntArray,
     fill_value: int,
     node_x: FloatArray,
     node_y: FloatArray,
-):
+) -> FloatArray:
     nodes = np.column_stack([node_x, node_y])
     closed, _ = close_polygons(face_node_connectivity, fill_value)
     coordinates = nodes[closed]
-    xy0 = coordinates[:, 0]
-    a = coordinates[:, :-1] - xy0[:, np.newaxis]
-    b = coordinates[:, 1:] - xy0[:, np.newaxis]
-    determinant = np.cross(a, b)
-    return 0.5 * abs(determinant.sum(axis=1))
+    return area_from_coordinates(coordinates)
 
 
 def centroids(
diff --git a/xugrid/ugrid/voronoi.py b/xugrid/ugrid/voronoi.py
index 7e7d41bd5..263b792f5 100644
--- a/xugrid/ugrid/voronoi.py
+++ b/xugrid/ugrid/voronoi.py
@@ -21,7 +21,12 @@
 from scipy import sparse
 
 from xugrid.constants import X_EPSILON, FloatArray, IntArray
-from xugrid.ugrid.connectivity import renumber
+from xugrid.ugrid.connectivity import (
+    area_from_coordinates,
+    close_polygons,
+    ragged_index,
+    renumber,
+)
 
 
 def dot_product2d(U: FloatArray, V: FloatArray):
@@ -134,6 +139,42 @@ def exterior_vertices(
     return i_keep, j_keep, vertices_keep, face_i, n_interpolated
 
 
+def choose_convex(
+    i: IntArray,
+    j: IntArray,
+    nodes: FloatArray,
+    original_vertices: FloatArray,
+    n_interpolated: int,
+) -> None:
+    # Determine whether the original vertex or the interpolated vertex
+    # generates the largest area (since the concave face will be smaller than
+    # the convex face.)
+    # Create a face_node_connectivity array.
+    n_vertex = np.bincount(i)
+    n_vertex = n_vertex[n_vertex > 0]
+    n = len(n_vertex)
+    m = n_vertex.max()
+    index = ragged_index(n, m, n_vertex)
+    faces = np.full((n, m), -1)
+    faces[index] = j
+    # Close the polygons so we can easily compute areas.
+    closed, _ = close_polygons(faces, -1)
+    # Make a copy and insert the original vertices.
+    modified_nodes = nodes.copy()
+    modified_nodes[-n_interpolated:] = original_vertices
+
+    # Compare areas of faces
+    convex_area = area_from_coordinates(nodes[closed])
+    modified_area = area_from_coordinates(modified_nodes[closed])
+    original_is_convex = (modified_area >= convex_area)[:, np.newaxis]
+    # All the newly created vertices are found at the end.
+    is_interpolated = faces >= len(nodes) - n_interpolated
+    # No need for unique: every exterior vertex is featured exactly once.
+    use_original = faces[original_is_convex & is_interpolated]
+    nodes[use_original] = modified_nodes[use_original]
+    return
+
+
 def exterior_topology(
     edge_face_connectivity: IntArray,
     edge_node_connectivity: IntArray,
@@ -142,6 +183,7 @@ def exterior_topology(
     vertices: FloatArray,
     centroids: FloatArray,
     add_vertices: bool,
+    skip_concave: bool,
 ):
     """
     Create the exterior topology of the voronoi tesselation.
@@ -220,7 +262,10 @@ def exterior_topology(
     # to generate the proper ordering. We overwrite those substituted points
     # here by their possibly concave true vertices.
     if add_vertices:
-        vor_vertices[-n_interpolated:] = orig_vertices
+        if skip_concave:
+            choose_convex(i, j, vor_vertices, orig_vertices, n_interpolated)
+        else:
+            vor_vertices[-n_interpolated:] = orig_vertices
 
     return vor_vertices, i, j, face_i
 
@@ -234,6 +279,7 @@ def voronoi_topology(
     fill_value: int = None,
     add_exterior: bool = False,
     add_vertices: bool = False,
+    skip_concave: bool = False,
 ) -> Tuple[FloatArray, sparse.csr_matrix]:
     """
     Compute the centroidal voronoi tesslation (CVT) of an existing mesh of
@@ -277,6 +323,9 @@ def voronoi_topology(
         exclusively centroids.
     add_vertices: bool, optional
         Whether to use existing exterior vertices.
+    skip_concave: bool, optional
+        Whether to skip existing exterior vertices if they generate concave
+        faces.
 
     Returns
     -------
@@ -334,6 +383,7 @@ def voronoi_topology(
             vertices,
             centroids,
             add_vertices,
+            skip_concave,
         )
         offset = node_i.max() + 1 if len(node_i > 0) else 0
         i = np.concatenate([node_i, exterior_i + offset])
@@ -347,4 +397,5 @@ def voronoi_topology(
     i = renumber(i)
     coo_content = (j, (i, j))
     face_node_connectivity = sparse.coo_matrix(coo_content)
+
     return vor_vertices, face_node_connectivity, face_i

From 84faaf775e8caad0e485dbd21950728c404f668a Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Wed, 21 Aug 2024 16:28:07 +0200
Subject: [PATCH 2/8] Implement smooth interpolation along exterior edges

---
 xugrid/regrid/unstructured.py | 94 +++++++++++++++++++++++++++--------
 xugrid/ugrid/voronoi.py       | 34 ++++++++++---
 2 files changed, 100 insertions(+), 28 deletions(-)

diff --git a/xugrid/regrid/unstructured.py b/xugrid/regrid/unstructured.py
index 5fc514dee..a59525302 100644
--- a/xugrid/regrid/unstructured.py
+++ b/xugrid/regrid/unstructured.py
@@ -3,10 +3,51 @@
 
 import xugrid as xu
 from xugrid.constants import FloatDType
-from xugrid.ugrid import voronoi
+from xugrid.ugrid import connectivity, voronoi
 from xugrid.ugrid.ugrid2d import Ugrid2d
 
 
+def replace_interpolated_weights(
+    vertices,
+    faces,
+    face_index,
+    weights,
+    node_to_node_map,
+    node_index_threshold,
+):
+    n, m = weights.shape
+    for i in range(n):
+        face = faces[face_index[i]]
+        weights_row = weights[i]
+        for j in range(m):
+            p = face[j]
+            w = weights_row[j]
+            if (p < node_index_threshold) or (w <= 0):
+                continue
+            # Find the two surrounding nodes (q and r)
+            q, r = node_to_node_map[p - node_index_threshold]
+            px, py = vertices[p]
+            qx, qy = vertices[q]
+            rx, ry = vertices[r]
+            # Compute the euclidian distance to both
+            p_q = np.sqrt((qx - px) ** 2 + (qy - py) ** 2)
+            p_r = np.sqrt((rx - px) ** 2 + (ry - py) ** 2)
+            total = p_q + p_r
+            # Redistribute weight according to inverse distance.
+            weight_q = (p_r / total) * w
+            weight_r = (p_q / total) * w
+            # Set weights to zero for p, and add to r and q.
+            weights[i, j] = 0.0
+            # Search for p and q
+            for jj in range(m):
+                node = face[jj]
+                if node == q:
+                    weights[i, jj] += weight_q
+                if node == r:
+                    weights[i, jj] += weight_r
+    return
+
+
 class UnstructuredGrid2d:
     """
     Stores only the grid topology.
@@ -97,16 +138,23 @@ def barycentric(self, other):
         grid = self.ugrid_topology
 
         # Create a voronoi grid to get surrounding nodes as vertices
-        vertices, faces, node_to_face_index = voronoi.voronoi_topology(
+        (
+            vertices,
+            faces,
+            node_to_face_index,
+            node_to_node_map,
+        ) = voronoi.voronoi_topology(
             grid.node_face_connectivity,
             grid.node_coordinates,
             grid.centroids,
-            #    edge_face_connectivity=grid.edge_face_connectivity,
-            #    edge_node_connectivity=grid.edge_node_connectivity,
-            #    fill_value=grid.fill_value,
-            #    add_exterior=True,
-            #    add_vertices=False,
+            edge_face_connectivity=grid.edge_face_connectivity,
+            edge_node_connectivity=grid.edge_node_connectivity,
+            fill_value=grid.fill_value,
+            add_exterior=True,
+            add_vertices=True,
+            skip_concave=True,
         )
+        faces = connectivity.to_dense(faces, fill_value=-1)
 
         voronoi_grid = Ugrid2d(
             vertices[:, 0],
@@ -114,29 +162,31 @@ def barycentric(self, other):
             -1,
             faces,
         )
-
         face_index, weights = voronoi_grid.compute_barycentric_weights(points)
-        n_points, n_max_node = weights.shape
+
+        # Find which nodes are interpolated. Redistribute their weights
+        # according to distance to projection vertex.
+        replace_interpolated_weights(
+            vertices=vertices,
+            faces=faces,
+            face_index=face_index,
+            weights=weights,
+            node_to_node_map=node_to_node_map,
+            node_index_threshold=len(vertices) - len(node_to_node_map),
+        )
+
+        # Discards 0 weights and points that fall outside of the grid.
+        outside = grid.locate_points(points) == -1
+        weights[outside] = 0
         keep = weights.ravel() > 0
         source_index = node_to_face_index[
             voronoi_grid.face_node_connectivity[face_index]
         ].ravel()[keep]
+
+        n_points, n_max_node = weights.shape
         target_index = np.repeat(np.arange(n_points), n_max_node)[keep]
         weights = weights.ravel()[keep]
 
-        # Look for points falling outside.
-        outside = face_index == -1
-        other_points = points[outside]
-        sampled_index = grid.locate_points(other_points)
-        sampled_inside = sampled_index != grid.fill_value
-        other_source = sampled_index[sampled_inside]
-        other_target = np.arange(n_points)[outside][sampled_inside]
-
-        # Combine first and second
-        source_index = np.concatenate((source_index, other_source))
-        target_index = np.concatenate((target_index, other_target))
-        weights = np.concatenate((weights, np.ones(other_target.size, dtype=float)))
-
         order = np.argsort(target_index)
         return source_index[order], target_index[order], weights[order]
 
diff --git a/xugrid/ugrid/voronoi.py b/xugrid/ugrid/voronoi.py
index 263b792f5..4cf159826 100644
--- a/xugrid/ugrid/voronoi.py
+++ b/xugrid/ugrid/voronoi.py
@@ -118,6 +118,7 @@ def exterior_vertices(
     n = n_centroid + len(vertices_keep)
     j_keep = np.repeat(np.arange(n_centroid, n), 2)
     n_interpolated = 0
+    interpolation_map = None
 
     # We add a substitution value for the actual vertex
     if add_vertices:
@@ -136,7 +137,11 @@ def exterior_vertices(
         # this: these are associated with two faces. So we set a value of -1 here.
         face_i = np.concatenate([face_i, np.full(n_interpolated, -1)])
 
-    return i_keep, j_keep, vertices_keep, face_i, n_interpolated
+        # For the interpolated vertices, it depends on the two other nodes (which
+        # depend on centroids). Store for each interpolation on which two nodes it relies.
+        interpolation_map = jj.reshape((-1, 2))
+
+    return i_keep, j_keep, vertices_keep, face_i, n_interpolated, interpolation_map
 
 
 def choose_convex(
@@ -227,7 +232,14 @@ def exterior_topology(
         fill_value,
     )
     i1, j1 = exterior_centroids(node_face_connectivity)
-    i2, j2, projected_vertices, face_i, n_interpolated = exterior_vertices(
+    (
+        i2,
+        j2,
+        projected_vertices,
+        face_i,
+        n_interpolated,
+        interpolation_map,
+    ) = exterior_vertices(
         edge_face_connectivity,
         edge_node_connectivity,
         fill_value,
@@ -267,7 +279,7 @@ def exterior_topology(
         else:
             vor_vertices[-n_interpolated:] = orig_vertices
 
-    return vor_vertices, i, j, face_i
+    return vor_vertices, i, j, face_i, interpolation_map
 
 
 def voronoi_topology(
@@ -280,7 +292,7 @@ def voronoi_topology(
     add_exterior: bool = False,
     add_vertices: bool = False,
     skip_concave: bool = False,
-) -> Tuple[FloatArray, sparse.csr_matrix]:
+) -> Tuple[FloatArray, sparse.csr_matrix, IntArray, IntArray]:
     """
     Compute the centroidal voronoi tesslation (CVT) of an existing mesh of
     (convex!) cells using connectivity index arrays.
@@ -335,6 +347,9 @@ def voronoi_topology(
         Connects the nodes of the voronoi topology to the faces of the original
         grid. Exterior vertices (when ``add_vertices=True``) are given an index
         of -1.
+    interpolation_map: ndarray of ints with shape ``(n_interpolated, 2)``
+        Marks for each interpolated point from which nodes it has been
+        interpolated.
     """
     if add_exterior:
         if any(
@@ -375,7 +390,13 @@ def voronoi_topology(
     j = face_i[order]
 
     if add_exterior:
-        vor_vertices, exterior_i, exterior_j, face_i = exterior_topology(
+        (
+            vor_vertices,
+            exterior_i,
+            exterior_j,
+            face_i,
+            interpolation_map,
+        ) = exterior_topology(
             edge_face_connectivity,
             edge_node_connectivity,
             node_face_connectivity,
@@ -389,6 +410,7 @@ def voronoi_topology(
         i = np.concatenate([node_i, exterior_i + offset])
         j = np.concatenate([j, exterior_j])
     else:
+        interpolation_map = None
         vor_vertices = centroids[np.unique(face_i)]
         face_i = np.arange(face_i.max() + 1)
         i = node_i
@@ -398,4 +420,4 @@ def voronoi_topology(
     coo_content = (j, (i, j))
     face_node_connectivity = sparse.coo_matrix(coo_content)
 
-    return vor_vertices, face_node_connectivity, face_i
+    return vor_vertices, face_node_connectivity, face_i, interpolation_map

From 2f924f2fcc4b2014863aad50ee12bfe9b1e3c8e3 Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Wed, 21 Aug 2024 20:57:37 +0200
Subject: [PATCH 3/8] Fix structured regridding: do linear interpolation within
 bounds, not just between centroids.

---
 tests/test_regrid/test_structured.py | 83 +++++++++++++++++++---------
 xugrid/regrid/structured.py          | 58 ++++++++++++-------
 2 files changed, 95 insertions(+), 46 deletions(-)

diff --git a/tests/test_regrid/test_structured.py b/tests/test_regrid/test_structured.py
index ca6e8db9b..cd812c2d8 100644
--- a/tests/test_regrid/test_structured.py
+++ b/tests/test_regrid/test_structured.py
@@ -308,9 +308,9 @@ def test_linear_weights_1d(
     # --------
     assert_expected_overlap(
         *grid_data_a_1d.linear_weights(grid_data_c_1d),
-        np.array([1, 0, 2, 1]),
-        np.array([1, 1, 2, 2]),
-        np.array([0.8, 0.2, 0.8, 0.2]),
+        np.array([0, 0, 1, 0, 2, 1]),
+        np.array([0, 0, 1, 1, 2, 2]),
+        np.array([0.0, 1.0, 0.8, 0.2, 0.8, 0.2]),
     )
 
     # --------
@@ -324,9 +324,9 @@ def test_linear_weights_1d(
     # --------
     assert_expected_overlap(
         *grid_data_a_1d.linear_weights(grid_data_d_1d),
-        np.array([0, 1, 1, 0, 1, 2]),
-        np.array([1, 1, 2, 2, 3, 3]),
-        np.array([0.9, 0.1, 0.6, 0.4, 0.9, 0.1]),
+        np.array([0, 0, 0, 1, 1, 0, 1, 2]),
+        np.array([0, 0, 1, 1, 2, 2, 3, 3]),
+        np.array([0.0, 0.1, 0.9, 0.1, 0.6, 0.4, 0.9, 0.1]),
     )
 
     # non-equidistant
@@ -339,9 +339,9 @@ def test_linear_weights_1d(
     # --------
     assert_expected_overlap(
         *grid_data_a_1d.linear_weights(grid_data_e_1d),
-        np.array([0, 1, 1, 2]),
-        np.array([1, 1, 2, 2]),
-        np.array([0.65, 0.35, 0.9, 0.1]),
+        np.array([0, 0, 0, 1, 1, 2]),
+        np.array([0, 0, 1, 1, 2, 2]),
+        np.array([0.0, 1.0, 0.65, 0.35, 0.9, 0.1]),
     )
 
     # 1-1 grid
@@ -354,9 +354,9 @@ def test_linear_weights_1d(
     # --------
     assert_expected_overlap(
         *grid_data_b_1d.linear_weights(grid_data_b_1d),
-        np.array([1, 0, 2, 1]),
-        np.array([1, 1, 2, 2]),
-        np.array([1.0, 0.0, 1.0, 0.0]),
+        np.array([0, 0, 1, 0, 2, 1, 3, 2]),
+        np.array([0, 0, 1, 1, 2, 2, 3, 3]),
+        np.array([0.0, 1.0, 1.0, 0.0, 1.0, 0.0, 1.0, 0.0]),
     )
 
 
@@ -379,32 +379,63 @@ def test_linear_weights_2d(
 
     # --------
     # source   targets      weight
+    # 4        0, 0, 3, 3   0%  50% 0%  50%
     # 5        0, 1, 3, 4   10%	40%	10%	40%
     # 6        1, 2, 4, 5   10%	40%	10%	40%
+    # 8        3, 3, 6, 6   0%  50% 0%  50%
     # 9        3, 4, 6, 7   10%	40%	10%	40%
     # 10       4, 5, 7, 8   10%	40%	10%	40%
     # --------
     assert_expected_overlap(
         *grid_data_a_layered_2d.linear_weights(grid_data_c_2d),
-        np.array([1, 0, 4, 3, 2, 1, 5, 4, 4, 3, 7, 6, 5, 4, 8, 7]),
-        np.array([5, 5, 5, 5, 6, 6, 6, 6, 9, 9, 9, 9, 10, 10, 10, 10]),
-        np.array([0.4, 0.1, 0.4, 0.1] * 4),
+        np.array(
+            [0, 0, 3, 3, 1, 0, 3, 4, 5, 4, 2, 1, 3, 3, 6, 6, 4, 3, 7, 6, 8, 5, 4, 7]
+        ),
+        np.array(
+            [4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6, 8, 8, 8, 8, 9, 9, 9, 9, 10, 10, 10, 10]
+        ),
+        np.array(
+            [
+                0.0,
+                0.5,
+                0.0,
+                0.5,
+                0.4,
+                0.1,
+                0.1,
+                0.4,
+                0.4,
+                0.1,
+                0.4,
+                0.1,
+                0.0,
+                0.5,
+                0.0,
+                0.5,
+                0.4,
+                0.1,
+                0.4,
+                0.1,
+                0.4,
+                0.4,
+                0.1,
+                0.1,
+            ]
+        ),
     )
 
     # 1-1
     # --------
     # source   targets      weight
-    # 5        4, 5, 8, 9   0% 100% 0% 0%
-    # 6        5, 6, 9,10   0% 100% 0% 0%
-    # 9        8, 9,12,13   0% 100% 0% 0%
-    # 10       9,10,13,14   0% 100% 0% 0%
-    # --------
-    assert_expected_overlap(
-        *grid_data_b_2d.linear_weights(grid_data_b_2d),
-        np.array([5, 4, 9, 8, 6, 5, 10, 9, 9, 8, 13, 12, 10, 9, 14, 13]),
-        np.array([5, 5, 5, 5, 6, 6, 6, 6, 9, 9, 9, 9, 10, 10, 10, 10]),
-        np.array([1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 0]),
-    )
+    # 0-15     0-15         0% 0% 0% 100% (or shuffled)
+    # result should be 1:1 mapping
+    # --------
+    source, target, weights = grid_data_b_2d.linear_weights(grid_data_b_2d)
+    expected_target = np.repeat(np.arange(16), 4)
+    assert np.array_equal(target, expected_target)
+    assert np.array_equal(np.unique(weights), [0, 1])
+    check_source = source[weights != 0]
+    assert np.array_equal(check_source, np.arange(16))
 
 
 def test_nonscalar_dx():
diff --git a/xugrid/regrid/structured.py b/xugrid/regrid/structured.py
index 8b14044f8..7ebf54ebe 100644
--- a/xugrid/regrid/structured.py
+++ b/xugrid/regrid/structured.py
@@ -6,6 +6,7 @@
 While the unstructured logic would work for structured data as well, it is much
 less efficient than utilizing the structure of the coordinates.
 """
+
 from typing import Any, Tuple, Union
 
 import numpy as np
@@ -176,6 +177,7 @@ def valid_nodes_within_bounds_and_extend(
             valid target indexes
         """
         source_index, target_index = self.valid_nodes_within_bounds(other)
+        return source_index, target_index
         valid = (other.midpoints[target_index] > self.midpoints[0]) & (
             other.midpoints[target_index] < self.midpoints[-1]
         )
@@ -283,31 +285,47 @@ def compute_linear_weights_to_centroids(
         weights: np.array
         neighbor: np.narray
         """
-        source_midpoint_index = self.maybe_reverse_index(source_index)
-        target_midpoints_index = other.maybe_reverse_index(target_index)
-        neighbor = np.ones(target_midpoints_index.size, dtype=int)
-        # cases where midpoint target <= midpoint source
-        condition = (
-            other.midpoints[target_midpoints_index]
-            <= self.midpoints[source_midpoint_index]
-        )
-        neighbor[condition] = -neighbor[condition]
-
         if self.midpoints.size < 2:
             raise ValueError(
                 f"Coordinate {self.name} has size: {self.midpoints.size}. "
                 "At least two points are required for interpolation."
             )
-        weights = 1 - (
-            (
-                other.midpoints[target_midpoints_index]
-                - self.midpoints[source_midpoint_index]
-            )
-            / (
-                self.midpoints[source_midpoint_index + neighbor]
-                - self.midpoints[source_midpoint_index]
-            )
-        )
+
+        source_midpoint_index = self.maybe_reverse_index(source_index)
+        target_midpoint_index = other.maybe_reverse_index(target_index)
+        # cases where midpoint target <= midpoint source: set neighbor to -1
+        neighbor = np.where(
+            other.midpoints[target_midpoint_index]
+            <= self.midpoints[source_midpoint_index],
+            -1,
+            1,
+        )
+
+        # When we exceed the original domain, it should still interpolate
+        # within the bounds.
+        # Make sure neighbor falls in [0, n)
+        neighbor_index = np.clip(
+            source_midpoint_index + neighbor, 0, self.midpoints.size - 1
+        )
+        # Update neighbor since we return it
+        neighbor = neighbor_index - source_midpoint_index
+
+        # If neighbor is 0, we end up computing zero distance, since we're
+        # comparing a midpoint to iself. Instead, set a weight of 1.0 on one,
+        # (and impliclity 0 in the other). Similarly, if source and target
+        # midpoints coincide, the distance may end up 0.
+        length = (
+            other.midpoints[target_midpoint_index]
+            - self.midpoints[source_midpoint_index]
+        )
+        total_length = (
+            self.midpoints[neighbor_index] - self.midpoints[source_midpoint_index]
+        )
+        # Do not divide by zero.
+        # We will overwrite the value anyway at neighbor == 0.
+        total_length[total_length == 0] = 1
+        weights = 1 - (length / total_length)
+        weights[neighbor == 0] = 0.0
         condition = np.logical_and(weights < 0.0, weights > 1.0)
         if condition.any():
             raise ValueError(

From 9dd26540486e278fdf5e2b230cd438cf245e10a2 Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Wed, 21 Aug 2024 21:00:09 +0200
Subject: [PATCH 4/8] Add support to skip concave vertices in voronoi
 tesselation. Use this in barycentric interpolation to interpolate beyond the
 centroids.

TODO: fix edge cases in barycentric interpolation. Might be numba_celltree?
---
 examples-dev/voronoi.py             |  74 +++++++++++------
 tests/test_regrid/test_regridder.py |  36 +++++++++
 tests/test_ugrid2d.py               |  20 +++--
 tests/test_voronoi.py               |  73 +++++++++++++----
 xugrid/regrid/unstructured.py       |  10 ++-
 xugrid/ugrid/ugrid2d.py             |  26 ++++--
 xugrid/ugrid/voronoi.py             | 120 +++++++++++++++++++---------
 7 files changed, 264 insertions(+), 95 deletions(-)

diff --git a/examples-dev/voronoi.py b/examples-dev/voronoi.py
index 1153dc333..758ea5a70 100644
--- a/examples-dev/voronoi.py
+++ b/examples-dev/voronoi.py
@@ -133,14 +133,13 @@ def comparison_plot(
 centroids = vertices[faces].mean(axis=1)
 
 node_face_connectivity = connectivity.invert_dense_to_sparse(faces, -1)
-voronoi_vertices, voronoi_faces, face_index = voronoi.voronoi_topology(
+voronoi_vertices, voronoi_faces, face_index, _ = voronoi.voronoi_topology(
     node_face_connectivity,
     vertices,
     centroids,
     add_exterior=False,
     add_vertices=False,
 )
-voronoi_faces = connectivity.to_dense(voronoi_faces, -1)
 
 # %%
 # We can compare the two meshes:
@@ -163,7 +162,7 @@ def comparison_plot(
     faces, -1
 )
 edge_face_connectivity = connectivity.invert_dense(face_edge_connectivity, -1)
-voronoi_vertices, voronoi_faces, face_index = voronoi.voronoi_topology(
+voronoi_vertices, voronoi_faces, face_index, _ = voronoi.voronoi_topology(
     node_face_connectivity,
     vertices,
     centroids,
@@ -173,7 +172,6 @@ def comparison_plot(
     add_exterior=True,
     add_vertices=True,
 )
-voronoi_faces = connectivity.to_dense(voronoi_faces, -1)
 
 comparison_plot(vertices, faces, centroids, voronoi_vertices, voronoi_faces)
 
@@ -197,7 +195,7 @@ def comparison_plot(
     faces, -1
 )
 edge_face_connectivity = connectivity.invert_dense(face_edge_connectivity, -1)
-voronoi_vertices, voronoi_faces, face_index = voronoi.voronoi_topology(
+voronoi_vertices, voronoi_faces, face_index, _ = voronoi.voronoi_topology(
     node_face_connectivity,
     vertices,
     centroids,
@@ -207,7 +205,6 @@ def comparison_plot(
     add_exterior=True,
     add_vertices=True,
 )
-voronoi_faces = connectivity.to_dense(voronoi_faces, -1)
 
 comparison_plot(vertices, faces, centroids, voronoi_vertices, voronoi_faces)
 
@@ -225,7 +222,7 @@ def comparison_plot(
     faces, -1
 )
 edge_face_connectivity = connectivity.invert_dense(face_edge_connectivity, -1)
-voronoi_vertices, voronoi_faces, face_index = voronoi.voronoi_topology(
+voronoi_vertices, voronoi_faces, face_index, _ = voronoi.voronoi_topology(
     node_face_connectivity,
     vertices,
     centroids,
@@ -235,25 +232,44 @@ def comparison_plot(
     add_exterior=True,
     add_vertices=False,
 )
-voronoi_faces = connectivity.to_dense(voronoi_faces, -1)
 
 comparison_plot(vertices, faces, centroids, voronoi_vertices, voronoi_faces)
 
 # %%
 # This will (obviously) result in a mesh that does not preserve the exterior
-# exactly.
-#
-# These are the three options, side by side:
+# exactly. Alternatively, we can choose to skip the exterior vertex if it
+# creates a concave face:
 
-nodes0, faces0, face_index0 = voronoi.voronoi_topology(
+node_face_connectivity = connectivity.invert_dense_to_sparse(faces, -1)
+edge_node_connectivity, face_edge_connectivity = connectivity.edge_connectivity(
+    faces, -1
+)
+edge_face_connectivity = connectivity.invert_dense(face_edge_connectivity, -1)
+voronoi_vertices, voronoi_faces, face_index, _ = voronoi.voronoi_topology(
+    node_face_connectivity,
+    vertices,
+    centroids,
+    edge_face_connectivity=edge_face_connectivity,
+    edge_node_connectivity=edge_node_connectivity,
+    fill_value=-1,
+    add_exterior=True,
+    add_vertices=True,
+    skip_concave=True,
+)
+
+comparison_plot(vertices, faces, centroids, voronoi_vertices, voronoi_faces)
+
+# %%
+# These are the four options, side by side:
+
+nodes0, faces0, face_index0, _ = voronoi.voronoi_topology(
     node_face_connectivity,
     vertices,
     centroids,
 )
-faces0 = connectivity.to_dense(faces0, -1)
 edges0, _ = connectivity.edge_connectivity(faces0, -1)
 
-nodes1, faces1, face_index1 = voronoi.voronoi_topology(
+nodes1, faces1, face_index1, _ = voronoi.voronoi_topology(
     node_face_connectivity,
     vertices,
     centroids,
@@ -263,10 +279,9 @@ def comparison_plot(
     add_exterior=True,
     add_vertices=False,
 )
-faces1 = connectivity.to_dense(faces1, -1)
 edges1, _ = connectivity.edge_connectivity(faces1, -1)
 
-nodes2, faces2, _ = voronoi.voronoi_topology(
+nodes2, faces2, _, _ = voronoi.voronoi_topology(
     node_face_connectivity,
     vertices,
     centroids,
@@ -276,20 +291,31 @@ def comparison_plot(
     add_exterior=True,
     add_vertices=True,
 )
-faces2 = connectivity.to_dense(faces2, -1)
 edges2, _ = connectivity.edge_connectivity(faces2, -1)
 
+nodes3, faces3, face_index3, node_map3 = voronoi.voronoi_topology(
+    node_face_connectivity,
+    vertices,
+    centroids,
+    edge_face_connectivity=edge_face_connectivity,
+    edge_node_connectivity=edge_node_connectivity,
+    fill_value=-1,
+    add_exterior=True,
+    add_vertices=True,
+    skip_concave=True,
+)
+edges3, _ = connectivity.edge_connectivity(faces3, -1)
 
 fig, axes = plt.subplots(
     nrows=1,
-    ncols=3,
-    figsize=(15, 5),
+    ncols=4,
+    figsize=(20, 5),
     subplot_kw={"box_aspect": 1},
     sharey=True,
     sharex=True,
 )
-all_edges = [edges0, edges1, edges2]
-all_nodes = [nodes0, nodes1, nodes2]
+all_edges = [edges0, edges1, edges2, edges3]
+all_nodes = [nodes0, nodes1, nodes2, nodes3]
 for ax, e, v in zip(axes, all_edges, all_nodes):
     edge_plot(v, e, ax, colors="red")
     ax.scatter(*centroids.T, color="red")
@@ -312,9 +338,9 @@ def comparison_plot(
 # Before we can send the data of an unstructured mesh off to a plotting library
 # such as ``matplotlib``, we'll generally need to triangulate the mesh. We can
 # directly use the first two options, since the generated voronoi vertices
-# correspond directly to a cell face. This is not the case for the third option,
-# since it includes some vertices of the original mesh, which are connected to
-# two faces.
+# correspond directly to a cell face. This is not the case for the third or
+# fourth option, since it includes some vertices of the original mesh, which
+# are connected to multiple faces.
 
 triangles0, face_triangles0 = connectivity.triangulate(faces0, -1)
 triangulation0 = mtri.Triangulation(nodes0[:, 0], nodes0[:, 1], triangles0)
diff --git a/tests/test_regrid/test_regridder.py b/tests/test_regrid/test_regridder.py
index bd8c1046c..4c79c152f 100644
--- a/tests/test_regrid/test_regridder.py
+++ b/tests/test_regrid/test_regridder.py
@@ -323,3 +323,39 @@ def test_directional_dependence():
         result.append(regridder.regrid(source))
     first = result.pop(0)
     assert all(first.identical(item) for item in result)
+
+
+def test_barycentric_concave():
+    vertices = np.array(
+        [
+            [0.0, 0.0],
+            [3.0, 0.0],
+            [1.0, 1.0],
+            [0.0, 2.0],
+            [3.0, 2.0],
+        ]
+    )
+    faces = np.array(
+        [
+            [0, 1, 2],
+            [0, 2, 3],
+            [2, 4, 3],
+        ]
+    )
+    grid = xu.Ugrid2d(*vertices.T, -1, faces)
+
+    dx = 0.1
+    x = np.arange(0.0, 3.0, dx) + 0.5 * dx
+    y = np.arange(0.0, 2.0, dx) + 0.5 * dx
+    other = xr.DataArray(
+        data=np.ones((y.size, x.size)), coords={"y": y, "x": x}, dims=("y", "x")
+    )
+
+    uda = xu.UgridDataArray(
+        obj=xr.DataArray([2.0, 0.5, 2.0], dims=[grid.face_dimension]),
+        grid=grid,
+    )
+    regridder = xu.BarycentricInterpolator(source=uda, target=other)
+    result = regridder.regrid(uda)
+    assert result.min() >= 0.5
+    assert result.max() <= 2.0
diff --git a/tests/test_ugrid2d.py b/tests/test_ugrid2d.py
index 07a6e0d61..d2e7f24c7 100644
--- a/tests/test_ugrid2d.py
+++ b/tests/test_ugrid2d.py
@@ -525,15 +525,19 @@ def test_voronoi_topology():
     )
     expected_vertices = np.vstack([CENTROIDS, expected_exterior])
     assert np.allclose(vertices, expected_vertices)
-    assert isinstance(faces, sparse.coo_matrix)
-    expected_row = np.array(
-        [0, 0, 0, 0, 1, 1, 1, 2, 2, 2, 2, 3, 3, 3, 4, 4, 4, 4, 5, 5, 5, 5, 6, 6, 6, 6]
-    )
-    expected_col = np.array(
-        [0, 1, 3, 2, 4, 0, 5, 4, 6, 1, 0, 6, 7, 1, 5, 0, 2, 8, 1, 7, 9, 3, 2, 3, 9, 8]
+    assert isinstance(faces, np.ndarray)
+    expected_faces = np.array(
+        [
+            [0, 1, 3, 2],
+            [4, 0, 5, -1],
+            [4, 6, 1, 0],
+            [6, 7, 1, -1],
+            [5, 0, 2, 8],
+            [1, 7, 9, 3],
+            [2, 3, 9, 8],
+        ]
     )
-    assert np.array_equal(faces.row, expected_row)
-    assert np.array_equal(faces.col, expected_col)
+    assert np.array_equal(faces, expected_faces)
     assert np.array_equal(face_index, [0, 1, 2, 3, 0, 0, 1, 1, 2, 3])
 
 
diff --git a/tests/test_voronoi.py b/tests/test_voronoi.py
index 5d9e1d5ce..9c68c92f6 100644
--- a/tests/test_voronoi.py
+++ b/tests/test_voronoi.py
@@ -168,7 +168,14 @@ def test_interior_centroids(self):
         assert np.array_equal(actual_j, expected_j)
 
     def test_exterior_vertices(self):
-        _, _, actual_vertices, actual_face, n = voronoi.exterior_vertices(
+        (
+            _,
+            _,
+            actual_vertices,
+            actual_face,
+            n,
+            interpolation_map,
+        ) = voronoi.exterior_vertices(
             self.edge_face_connectivity,
             self.edge_node_connectivity,
             self.fill_value,
@@ -179,25 +186,25 @@ def test_exterior_vertices(self):
         assert n == 0
         assert np.allclose(rowsort(actual_vertices), self.exterior_vertices)
         assert np.isin(np.arange(6), actual_face).all()
+        assert interpolation_map is None
 
     def test_voronoi_topology(self):
-        vertices, faces, face_i = voronoi.voronoi_topology(
+        vertices, faces, face_i, _ = voronoi.voronoi_topology(
             self.node_face_connectivity,
             self.vertices,
             self.centroids,
         )
-        actual_faces = connectivity.to_dense(faces, fill_value=-1)
         expected_faces = np.array(
             [
                 [0, 1, 4, 3],
                 [1, 2, 5, 4],
             ]
         )
-        assert actual_faces.shape == (2, 4)
+        assert faces.shape == (2, 4)
         assert np.allclose(rowsort(vertices), self.expected_vertices)
         assert np.array_equal(face_i, [0, 1, 2, 3, 4, 5])
-        assert np.array_equal(actual_faces, expected_faces)
-        assert np.allclose(mesh_area(vertices, actual_faces), 2.0)
+        assert np.array_equal(faces, expected_faces)
+        assert np.allclose(mesh_area(vertices, faces), 2.0)
 
     def test_voronoi_topology__add_exterior(self):
         with pytest.raises(
@@ -210,7 +217,7 @@ def test_voronoi_topology__add_exterior(self):
                 add_exterior=True,
             )
 
-        vertices, faces, face_i = voronoi.voronoi_topology(
+        vertices, faces, face_i, _ = voronoi.voronoi_topology(
             self.node_face_connectivity,
             self.vertices,
             self.centroids,
@@ -219,16 +226,15 @@ def test_voronoi_topology__add_exterior(self):
             self.fill_value,
             add_exterior=True,
         )
-        actual_faces = connectivity.to_dense(faces, fill_value=-1)
         expected_vertices = rowsort(
             np.concatenate([self.expected_vertices, self.exterior_vertices])
         )
-        assert actual_faces.shape == (12, 4)
+        assert faces.shape == (12, 4)
         assert np.allclose(rowsort(vertices), expected_vertices)
         assert (face_i != -1).all()
-        assert np.allclose(mesh_area(vertices, actual_faces), 5.5)
+        assert np.allclose(mesh_area(vertices, faces), 5.5)
 
-        vertices, faces, face_i = voronoi.voronoi_topology(
+        vertices, faces, face_i, node_interpolation = voronoi.voronoi_topology(
             self.node_face_connectivity,
             self.vertices,
             self.centroids,
@@ -238,15 +244,14 @@ def test_voronoi_topology__add_exterior(self):
             add_exterior=True,
             add_vertices=True,
         )
-        actual_faces = connectivity.to_dense(faces, fill_value=-1)
         expected_vertices = rowsort(
             np.concatenate([expected_vertices, self.additional_vertices])
         )
         # This introduces hanging nodes
-        assert actual_faces.shape == (12, 5)
+        assert faces.shape == (12, 5)
         assert np.allclose(rowsort(vertices), expected_vertices)
         assert (face_i == -1).sum() == 10
-        assert np.allclose(mesh_area(vertices, actual_faces), 6.0)
+        assert np.allclose(mesh_area(vertices, faces), 6.0)
 
 
 def test_projected_vertices_on_edge():
@@ -334,3 +339,43 @@ def test_isolated_face():
     voronoi_grid = grid.tesselate_centroidal_voronoi(False, False)
     assert voronoi_grid.n_face == 2
     assert voronoi_grid.n_node == 4
+
+
+def test_concave_voronoi():
+    r"""
+    Three triangles.
+    The lack of a fourth on the creates
+
+    4 (0,2)      5 (3,2)
+         *------------*
+         | \   2    / |
+         |   \    /   |
+         |  1  \/     |
+         |     /\     |
+         |   /    \   |
+         | /   0    \ |
+         *------------*
+    1 (0,0)  3 (1,1)  2 (3,0)
+    """
+    vertices = np.array(
+        [
+            [0.0, 0.0],
+            [3.0, 0.0],
+            [1.0, 1.0],
+            [0.0, 2.0],
+            [3.0, 2.0],
+        ]
+    )
+    faces = np.array(
+        [
+            [0, 1, 2],
+            [0, 2, 3],
+            [2, 4, 3],
+        ]
+    )
+    grid = xu.Ugrid2d(*vertices.T, -1, faces)
+    voronoi_0 = grid.tesselate_centroidal_voronoi(skip_concave=False)
+    voronoi_1 = grid.tesselate_centroidal_voronoi(skip_concave=True)
+    assert voronoi_0.n_face == voronoi_1.n_face
+    assert voronoi_0.n_node == voronoi_1.n_node
+    assert voronoi_0.area.sum() < voronoi_1.area.sum()
diff --git a/xugrid/regrid/unstructured.py b/xugrid/regrid/unstructured.py
index a59525302..ab8a6a3e6 100644
--- a/xugrid/regrid/unstructured.py
+++ b/xugrid/regrid/unstructured.py
@@ -1,12 +1,14 @@
+import numba as nb
 import numpy as np
 import xarray as xr
 
 import xugrid as xu
 from xugrid.constants import FloatDType
-from xugrid.ugrid import connectivity, voronoi
+from xugrid.ugrid import voronoi
 from xugrid.ugrid.ugrid2d import Ugrid2d
 
 
+@nb.njit(cache=True)
 def replace_interpolated_weights(
     vertices,
     faces,
@@ -25,7 +27,8 @@ def replace_interpolated_weights(
             if (p < node_index_threshold) or (w <= 0):
                 continue
             # Find the two surrounding nodes (q and r)
-            q, r = node_to_node_map[p - node_index_threshold]
+            index = p - node_index_threshold
+            q, r = node_to_node_map[index]
             px, py = vertices[p]
             qx, qy = vertices[q]
             rx, ry = vertices[r]
@@ -36,7 +39,7 @@ def replace_interpolated_weights(
             # Redistribute weight according to inverse distance.
             weight_q = (p_r / total) * w
             weight_r = (p_q / total) * w
-            # Set weights to zero for p, and add to r and q.
+            # Set weights to zero for p, and add to r and q weights.
             weights[i, j] = 0.0
             # Search for p and q
             for jj in range(m):
@@ -154,7 +157,6 @@ def barycentric(self, other):
             add_vertices=True,
             skip_concave=True,
         )
-        faces = connectivity.to_dense(faces, fill_value=-1)
 
         voronoi_grid = Ugrid2d(
             vertices[:, 0],
diff --git a/xugrid/ugrid/ugrid2d.py b/xugrid/ugrid/ugrid2d.py
index 319cef5de..f9ae7c5fe 100644
--- a/xugrid/ugrid/ugrid2d.py
+++ b/xugrid/ugrid/ugrid2d.py
@@ -753,7 +753,7 @@ def voronoi_topology(self):
         face_index: 1d array of integers
         """
         if self._voronoi_topology is None:
-            vertices, faces, face_index = voronoi_topology(
+            vertices, faces, face_index, _ = voronoi_topology(
                 self.node_face_connectivity,
                 self.node_coordinates,
                 self.centroids,
@@ -1833,7 +1833,7 @@ def triangulate(self):
         )
         return Ugrid2d(self.node_x, self.node_y, self.fill_value, triangles)
 
-    def _tesselate_voronoi(self, centroids, add_exterior, add_vertices):
+    def _tesselate_voronoi(self, centroids, add_exterior, add_vertices, skip_concave):
         if add_exterior:
             edge_face_connectivity = self.edge_face_connectivity
             edge_node_connectivity = self.edge_node_connectivity
@@ -1841,7 +1841,7 @@ def _tesselate_voronoi(self, centroids, add_exterior, add_vertices):
             edge_face_connectivity = None
             edge_node_connectivity = None
 
-        vertices, faces, _ = voronoi_topology(
+        vertices, faces, _, _ = voronoi_topology(
             self.node_face_connectivity,
             self.node_coordinates,
             centroids,
@@ -1850,11 +1850,13 @@ def _tesselate_voronoi(self, centroids, add_exterior, add_vertices):
             self.fill_value,
             add_exterior,
             add_vertices,
+            skip_concave,
         )
-        faces = connectivity.to_dense(faces, self.fill_value)
         return Ugrid2d(vertices[:, 0], vertices[:, 1], self.fill_value, faces)
 
-    def tesselate_centroidal_voronoi(self, add_exterior=True, add_vertices=True):
+    def tesselate_centroidal_voronoi(
+        self, add_exterior=True, add_vertices=True, skip_concave=False
+    ):
         """
         Create a centroidal Voronoi tesselation of this UGRID2D topology.
 
@@ -1866,14 +1868,19 @@ def tesselate_centroidal_voronoi(self, add_exterior=True, add_vertices=True):
         ----------
         add_exterior: bool, default: True
         add_vertices: bool, default: True
+        skip_concave: bool, default: False
 
         Returns
         -------
         tesselation: Ugrid2d
         """
-        return self._tesselate_voronoi(self.centroids, add_exterior, add_vertices)
+        return self._tesselate_voronoi(
+            self.centroids, add_exterior, add_vertices, skip_concave
+        )
 
-    def tesselate_circumcenter_voronoi(self, add_exterior=True, add_vertices=True):
+    def tesselate_circumcenter_voronoi(
+        self, add_exterior=True, add_vertices=True, skip_concave=False
+    ):
         """
         Create a circumcenter Voronoi tesselation of this UGRID2D topology.
 
@@ -1885,12 +1892,15 @@ def tesselate_circumcenter_voronoi(self, add_exterior=True, add_vertices=True):
         ----------
         add_exterior: bool, default: True
         add_vertices: bool, default: True
+        skip_concave: bool, default: False
 
         Returns
         -------
         tesselation: Ugrid2d
         """
-        return self._tesselate_voronoi(self.circumcenters, add_exterior, add_vertices)
+        return self._tesselate_voronoi(
+            self.circumcenters, add_exterior, add_vertices, skip_concave
+        )
 
     def reverse_cuthill_mckee(self, dimension=None):
         """
diff --git a/xugrid/ugrid/voronoi.py b/xugrid/ugrid/voronoi.py
index 4cf159826..84dee3d7d 100644
--- a/xugrid/ugrid/voronoi.py
+++ b/xugrid/ugrid/voronoi.py
@@ -55,6 +55,17 @@ def compute_centroid(i: IntArray, x: FloatArray, y: FloatArray):
     return _centroid_pandas(i, x, y)
 
 
+def _create_face_node_connectivity(i: IntArray, j: IntArray) -> IntArray:
+    n_vertex = np.bincount(i)
+    n_vertex = n_vertex[n_vertex > 0]
+    n = len(n_vertex)
+    m = n_vertex.max()
+    index = ragged_index(n, m, n_vertex)
+    face_node_connectivity = np.full((n, m), -1)
+    face_node_connectivity[index] = j
+    return face_node_connectivity
+
+
 def exterior_centroids(node_face_connectivity: sparse.csr_matrix):
     n, _ = node_face_connectivity.shape
     # Find exterior nodes NOT associated with any interior edge
@@ -85,6 +96,53 @@ def interior_centroids(
     return i, j
 
 
+def _project_centroid_on_edge(
+    edge_vertices: FloatArray, centroid_vertices: FloatArray
+) -> FloatArray:
+    a = edge_vertices[:, 0, :]
+    b = edge_vertices[:, 1, :]
+    V = b - a
+    U = centroid_vertices - a
+    # np.dot(U, V) doesn't do the desired thing here
+    projected_vertices = a + ((dot_product2d(U, V) / dot_product2d(V, V)) * V.T).T
+    return projected_vertices
+
+
+def _filter_projections(
+    projected_vertices: FloatArray,
+    centroid_vertices: FloatArray,
+    exterior_nodes: IntArray,
+    face_i: IntArray,
+) -> Tuple[IntArray, FloatArray, IntArray]:
+    # Create the numbering pointing to these new vertices.
+    # Discard vertices that overlap with e.g. circumcenters.
+    keep = np.linalg.norm(projected_vertices - centroid_vertices, axis=1) > (
+        X_EPSILON * X_EPSILON
+    )
+    face_i = face_i[keep]
+    vertices = projected_vertices[keep]
+    i = exterior_nodes[keep].ravel()
+    return i, vertices, face_i
+
+
+def _interpolate_between_projections(
+    projected_vertices: FloatArray,
+    exterior_nodes: IntArray,
+    n: int,
+) -> Tuple[IntArray, IntArray, FloatArray, IntArray]:
+    n_new = len(projected_vertices)
+    i = exterior_nodes.ravel()
+    order = np.argsort(i)
+    jj = np.repeat(np.arange(n_new), 2)[order]
+    to_interpolate = projected_vertices[jj]
+    interpolated_vertices = 0.5 * (to_interpolate[::2] + to_interpolate[1::2])
+    # For the interpolated vertices, it depends on the two other nodes (which
+    # depend on centroids). Store for each interpolation on which two nodes it relies.
+    j = np.arange(n, n + len(interpolated_vertices))
+    interpolation_map = jj.reshape((-1, 2))
+    return i[order][::2], j, interpolated_vertices, interpolation_map
+
+
 def exterior_vertices(
     edge_face_connectivity: IntArray,
     edge_node_connectivity: IntArray,
@@ -99,49 +157,39 @@ def exterior_vertices(
     edge_vertices = vertices[exterior_nodes]
     face_i = edge_face_connectivity[is_exterior, 0]
     centroid_vertices = centroids[face_i]
-    a = edge_vertices[:, 0, :]
-    b = edge_vertices[:, 1, :]
-    V = b - a
-    U = centroid_vertices - a
-    # np.dot(U, V) doesn't do the desired thing here
-    projected_vertices = a + ((dot_product2d(U, V) / dot_product2d(V, V)) * V.T).T
-
-    # Create the numbering pointing to these new vertices.
-    # Discard vertices that overlap with e.g. circumcenters.
-    keep = np.linalg.norm(projected_vertices - centroid_vertices, axis=1) > (
-        X_EPSILON * X_EPSILON
+    projected_vertices = _project_centroid_on_edge(edge_vertices, centroid_vertices)
+    i, vertices, face_i = _filter_projections(
+        projected_vertices,
+        centroid_vertices,
+        exterior_nodes,
+        face_i,
     )
-    face_i = face_i[keep]
-    vertices_keep = projected_vertices[keep]
+
     n_centroid = len(centroids)
-    i_keep = exterior_nodes[keep].ravel()
-    n = n_centroid + len(vertices_keep)
-    j_keep = np.repeat(np.arange(n_centroid, n), 2)
+    n = n_centroid + len(vertices)
+    j = np.repeat(np.arange(n_centroid, n), 2)
     n_interpolated = 0
     interpolation_map = None
 
-    # We add a substitution value for the actual vertex
+    # We add a guaranteed convex substitution vertex instead of the actual vertex
     if add_vertices:
-        n_new = len(projected_vertices)
-        i = exterior_nodes.ravel()
-        order = np.argsort(i)
-        jj = np.repeat(np.arange(n_new), 2)[order]
-        to_interpolate = projected_vertices[jj]
-        interpolated = 0.5 * (to_interpolate[::2] + to_interpolate[1::2])
-        n_interpolated = len(interpolated)
-        i_keep = np.concatenate([i_keep, i[order][::2]])
-        j_keep = np.concatenate([j_keep, np.arange(n, n + n_interpolated)])
-        vertices_keep = np.concatenate([vertices_keep, interpolated])
+        (
+            i_new,
+            j_new,
+            interpolated_vertices,
+            interpolation_map,
+        ) = _interpolate_between_projections(projected_vertices, exterior_nodes, n)
+        interpolation_map += n_centroid
+        n_interpolated = len(interpolated_vertices)
+        i = np.concatenate([i, i_new])
+        j = np.concatenate([j, j_new])
+        vertices = np.concatenate([vertices, interpolated_vertices])
         # Create face index: the face of the original mesh associated with every
         # voronoi vertex is a centroid. The exterior vertices are an exception to
         # this: these are associated with two faces. So we set a value of -1 here.
         face_i = np.concatenate([face_i, np.full(n_interpolated, -1)])
 
-        # For the interpolated vertices, it depends on the two other nodes (which
-        # depend on centroids). Store for each interpolation on which two nodes it relies.
-        interpolation_map = jj.reshape((-1, 2))
-
-    return i_keep, j_keep, vertices_keep, face_i, n_interpolated, interpolation_map
+    return i, j, vertices, face_i, n_interpolated, interpolation_map
 
 
 def choose_convex(
@@ -220,7 +268,8 @@ def exterior_topology(
 
     * i refers to the node number of the original mesh.
     * In principle, j refers to the face number (and associated centroids) of
-    the original mesh; exterior vertices require new numbers.
+    the original mesh; exterior vertices require new numbers as they are not a
+    centroid.
 
     Finally, in the resulting new (voronoi) face_node_connectivity, i becomes
     the face number, and j becomes the node number.
@@ -416,8 +465,5 @@ def voronoi_topology(
         i = node_i
         j = renumber(j)
 
-    i = renumber(i)
-    coo_content = (j, (i, j))
-    face_node_connectivity = sparse.coo_matrix(coo_content)
-
+    face_node_connectivity = _create_face_node_connectivity(i, j)
     return vor_vertices, face_node_connectivity, face_i, interpolation_map

From a14e4501925831754628fd774661ed0895908e8b Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Thu, 22 Aug 2024 16:36:22 +0200
Subject: [PATCH 5/8] Update tests for changes in numba_celltree, add test to
 check identical topology in structured vs unstructured form

---
 tests/test_regrid/test_regridder.py | 39 +++++++++++++++++++++++++++++
 tests/test_ugrid2d.py               |  8 +++---
 2 files changed, 43 insertions(+), 4 deletions(-)

diff --git a/tests/test_regrid/test_regridder.py b/tests/test_regrid/test_regridder.py
index 4c79c152f..dfe717535 100644
--- a/tests/test_regrid/test_regridder.py
+++ b/tests/test_regrid/test_regridder.py
@@ -359,3 +359,42 @@ def test_barycentric_concave():
     result = regridder.regrid(uda)
     assert result.min() >= 0.5
     assert result.max() <= 2.0
+    # There should be a large enough "gap" of NaNs on the right side.
+    assert result.isnull().sum() == 200
+
+
+def test_barycentric_structured():
+    da = xr.DataArray(
+        data=[[1, 2], [3, 4]],
+        coords={"y": [1.5, 0.5], "x": [0.5, 1.5]},
+        dims=("y", "x"),
+    )
+    x = np.arange(0.0, 2.0, 0.25) + 0.125
+    y = np.arange(0.0, 2.0, 0.25) + 0.125
+    target = xr.DataArray(
+        data=np.empty((y.size, x.size)), coords={"y": y, "x": x}, dims=("y", "x")
+    )
+
+    regridder = xu.BarycentricInterpolator(source=da, target=target)
+    out_structured = regridder.regrid(da)
+
+    target_uda = xu.UgridDataArray.from_structured(target)
+    regridder = xu.BarycentricInterpolator(source=da, target=target_uda)
+    out_unstructured = regridder.regrid(da)
+
+    # Answers should be same regardless whether the topology is stored as
+    # structured or unstructured.
+    values_structured = out_structured.to_numpy()
+    values_unstructured = out_unstructured.to_numpy()
+    assert np.allclose(values_structured.ravel(), values_unstructured)
+
+    # Check values. Original values should be present at exterior (repeated).
+    assert np.allclose(
+        values_structured[0], [3.0, 3.0, 3.125, 3.375, 3.625, 3.875, 4.0, 4.0]
+    )
+    assert np.allclose(
+        values_structured[-1], [1.0, 1.0, 1.125, 1.375, 1.625, 1.875, 2.0, 2.0]
+    )
+    assert np.allclose(
+        values_structured[:, 0], [3.0, 3.0, 2.75, 2.25, 1.75, 1.25, 1.0, 1.0]
+    )
diff --git a/tests/test_ugrid2d.py b/tests/test_ugrid2d.py
index d2e7f24c7..1809f01b2 100644
--- a/tests/test_ugrid2d.py
+++ b/tests/test_ugrid2d.py
@@ -643,8 +643,8 @@ def test_rasterize():
     x, y, index = grid.rasterize(resolution=0.5)
     expected_index = np.array(
         [
-            [-1, 2, -1, -1],
-            [2, 2, 3, -1],
+            [-1, 2, 3, -1],
+            [2, 2, 3, 3],
             [0, 0, 1, 1],
             [0, 0, 1, 1],
         ]
@@ -658,8 +658,8 @@ def test_rasterize():
     x, y, index = grid.rasterize(resolution=0.5, bounds=bounds)
     expected_index = np.array(
         [
-            [-1, -1, -1, 2, -1, -1],
-            [-1, -1, 2, 2, 3, -1],
+            [-1, -1, -1, 2, 3, -1],
+            [-1, -1, 2, 2, 3, 3],
             [-1, -1, 0, 0, 1, 1],
             [-1, -1, 0, 0, 1, 1],
             [-1, -1, -1, -1, -1, -1],

From d8185cc08ec01880d1cc87451939eec9d5491f52 Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Thu, 22 Aug 2024 16:52:06 +0200
Subject: [PATCH 6/8] Update changelog

---
 docs/changelog.rst | 11 +++++++++++
 1 file changed, 11 insertions(+)

diff --git a/docs/changelog.rst b/docs/changelog.rst
index 57b79a66d..f0ee42371 100644
--- a/docs/changelog.rst
+++ b/docs/changelog.rst
@@ -6,6 +6,17 @@ All notable changes to this project will be documented in this file.
 The format is based on `Keep a Changelog`_, and this project adheres to
 `Semantic Versioning`_.
 
+Unreleased
+----------
+
+Fixed
+~~~~~
+
+- The :class:`xugrid.BarycentricInterpolator` now interpolates according to
+  linear weights within the full bounds of the source grid, rather than only
+  within the centroids of the source grid. Previously, it would give no results
+  beyond the centroids for structured to structured regridding, and it would
+  give nearest results (equal to :class:`CentroidLocatorRegridder`) otherwise.
 
 [0.11.2] 2024-08-16
 -------------------

From 38467a0444fd0fb7fbc8ee0ac0e2934f08313a9d Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Fri, 23 Aug 2024 11:59:09 +0200
Subject: [PATCH 7/8] Latest numba_celltree release

---
 pixi.lock      | 56 +++++++++++++++++++++++++-------------------------
 pyproject.toml |  2 +-
 2 files changed, 29 insertions(+), 29 deletions(-)

diff --git a/pixi.lock b/pixi.lock
index ad2858425..3c0f3591b 100644
--- a/pixi.lock
+++ b/pixi.lock
@@ -231,7 +231,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/linux-64/nspr-4.35-h27087fc_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/nss-3.103-h593d115_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/numba-0.60.0-py312h83e6fd3_0.conda
-      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.1.8-pyhd8ed1ab_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.2.0-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/numcodecs-0.12.1-py312h7070661_1.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/numpy-1.26.4-py312heda63a1_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/openjpeg-2.5.2-h488ebb8_0.conda
@@ -560,7 +560,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/osx-64/nspr-4.35-hea0b92c_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-64/nss-3.103-he7eb89d_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-64/numba-0.60.0-py312hc3b515d_0.conda
-      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.1.8-pyhd8ed1ab_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.2.0-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-64/numcodecs-0.12.1-py312h28f332c_1.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-64/numpy-1.26.4-py312he3a82b2_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-64/openjpeg-2.5.2-h7310d3a_0.conda
@@ -878,7 +878,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/nspr-4.35-hb7217d7_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/nss-3.103-hc42bcbf_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/numba-0.60.0-py312h41cea2d_0.conda
-      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.1.8-pyhd8ed1ab_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.2.0-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/numcodecs-0.12.1-py312h5c2e7bc_1.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/numpy-1.26.4-py312h8442bc7_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/openjpeg-2.5.2-h9f1df11_0.conda
@@ -1200,7 +1200,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/win-64/nh3-0.2.18-py312h68c23d6_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/nodeenv-1.9.1-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/numba-0.60.0-py312hcccf92d_0.conda
-      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.1.8-pyhd8ed1ab_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.2.0-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/numcodecs-0.12.1-py312h275cf98_1.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/numpy-1.26.4-py312h8753938_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/openjpeg-2.5.2-h3d672ee_0.conda
@@ -1561,7 +1561,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/linux-64/nspr-4.35-h27087fc_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/nss-3.103-h593d115_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/numba-0.60.0-py39h0320e7d_0.conda
-      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.1.8-pyhd8ed1ab_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.2.0-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/numcodecs-0.12.1-py39h84cc369_1.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/numpy-1.26.4-py39h474f0d3_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/openjpeg-2.5.2-h488ebb8_0.conda
@@ -1892,7 +1892,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/osx-64/nspr-4.35-hea0b92c_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-64/nss-3.103-he7eb89d_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-64/numba-0.60.0-py39h90d9702_0.conda
-      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.1.8-pyhd8ed1ab_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.2.0-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-64/numcodecs-0.12.1-py39h09c4c31_1.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-64/numpy-1.26.4-py39h28c39a1_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-64/openjpeg-2.5.2-h7310d3a_0.conda
@@ -2212,7 +2212,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/nspr-4.35-hb7217d7_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/nss-3.103-hc42bcbf_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/numba-0.60.0-py39h2d4ef1e_0.conda
-      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.1.8-pyhd8ed1ab_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.2.0-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/numcodecs-0.12.1-py39hbf7db11_1.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/numpy-1.26.4-py39h7aa2656_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/openjpeg-2.5.2-h9f1df11_0.conda
@@ -2536,7 +2536,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/win-64/nh3-0.2.18-py39h7d6e1d9_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/nodeenv-1.9.1-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/numba-0.60.0-py39h5dcb127_0.conda
-      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.1.8-pyhd8ed1ab_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.2.0-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/numcodecs-0.12.1-py39ha51f57c_1.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/numpy-1.26.4-py39hddb5d58_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/openjpeg-2.5.2-h3d672ee_0.conda
@@ -2897,7 +2897,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/linux-64/nspr-4.35-h27087fc_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/nss-3.103-h593d115_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/numba-0.60.0-py310h5dc88bb_0.conda
-      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.1.8-pyhd8ed1ab_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.2.0-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/numcodecs-0.12.1-py310h76e45a6_1.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/numpy-1.26.4-py310hb13e2d6_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/openjpeg-2.5.2-h488ebb8_0.conda
@@ -3227,7 +3227,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/osx-64/nspr-4.35-hea0b92c_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-64/nss-3.103-he7eb89d_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-64/numba-0.60.0-py310h89a1501_0.conda
-      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.1.8-pyhd8ed1ab_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.2.0-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-64/numcodecs-0.12.1-py310he0a0c5d_1.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-64/numpy-1.26.4-py310h4bfa8fc_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-64/openjpeg-2.5.2-h7310d3a_0.conda
@@ -3546,7 +3546,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/nspr-4.35-hb7217d7_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/nss-3.103-hc42bcbf_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/numba-0.60.0-py310h0628f0e_0.conda
-      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.1.8-pyhd8ed1ab_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.2.0-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/numcodecs-0.12.1-py310hcf9f62a_1.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/numpy-1.26.4-py310hd45542a_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/openjpeg-2.5.2-h9f1df11_0.conda
@@ -3869,7 +3869,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/win-64/nh3-0.2.18-py310h4733a37_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/nodeenv-1.9.1-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/numba-0.60.0-py310h7793332_0.conda
-      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.1.8-pyhd8ed1ab_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.2.0-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/numcodecs-0.12.1-py310h9e98ed7_1.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/numpy-1.26.4-py310hf667824_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/openjpeg-2.5.2-h3d672ee_0.conda
@@ -4230,7 +4230,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/linux-64/nspr-4.35-h27087fc_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/nss-3.103-h593d115_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/numba-0.60.0-py311h4bc866e_0.conda
-      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.1.8-pyhd8ed1ab_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.2.0-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/numcodecs-0.12.1-py311h4332511_1.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/numpy-1.26.4-py311h64a7726_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/openjpeg-2.5.2-h488ebb8_0.conda
@@ -4559,7 +4559,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/osx-64/nspr-4.35-hea0b92c_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-64/nss-3.103-he7eb89d_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-64/numba-0.60.0-py311h0e5bd6a_0.conda
-      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.1.8-pyhd8ed1ab_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.2.0-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-64/numcodecs-0.12.1-py311hbafa61a_1.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-64/numpy-1.26.4-py311hc43a94b_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-64/openjpeg-2.5.2-h7310d3a_0.conda
@@ -4877,7 +4877,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/nspr-4.35-hb7217d7_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/nss-3.103-hc42bcbf_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/numba-0.60.0-py311h9506ed5_0.conda
-      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.1.8-pyhd8ed1ab_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.2.0-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/numcodecs-0.12.1-py311hb9542d7_1.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/numpy-1.26.4-py311h7125741_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/openjpeg-2.5.2-h9f1df11_0.conda
@@ -5199,7 +5199,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/win-64/nh3-0.2.18-py311h9363f20_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/nodeenv-1.9.1-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/numba-0.60.0-py311h0673bce_0.conda
-      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.1.8-pyhd8ed1ab_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.2.0-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/numcodecs-0.12.1-py311hda3d55a_1.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/numpy-1.26.4-py311h0b4df5a_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/openjpeg-2.5.2-h3d672ee_0.conda
@@ -5559,7 +5559,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/linux-64/nspr-4.35-h27087fc_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/nss-3.103-h593d115_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/numba-0.60.0-py312h83e6fd3_0.conda
-      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.1.8-pyhd8ed1ab_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.2.0-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/numcodecs-0.12.1-py312h7070661_1.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/numpy-1.26.4-py312heda63a1_0.conda
       - conda: https://conda.anaconda.org/conda-forge/linux-64/openjpeg-2.5.2-h488ebb8_0.conda
@@ -5888,7 +5888,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/osx-64/nspr-4.35-hea0b92c_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-64/nss-3.103-he7eb89d_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-64/numba-0.60.0-py312hc3b515d_0.conda
-      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.1.8-pyhd8ed1ab_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.2.0-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-64/numcodecs-0.12.1-py312h28f332c_1.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-64/numpy-1.26.4-py312he3a82b2_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-64/openjpeg-2.5.2-h7310d3a_0.conda
@@ -6206,7 +6206,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/nspr-4.35-hb7217d7_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/nss-3.103-hc42bcbf_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/numba-0.60.0-py312h41cea2d_0.conda
-      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.1.8-pyhd8ed1ab_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.2.0-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/numcodecs-0.12.1-py312h5c2e7bc_1.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/numpy-1.26.4-py312h8442bc7_0.conda
       - conda: https://conda.anaconda.org/conda-forge/osx-arm64/openjpeg-2.5.2-h9f1df11_0.conda
@@ -6528,7 +6528,7 @@ environments:
       - conda: https://conda.anaconda.org/conda-forge/win-64/nh3-0.2.18-py312h68c23d6_0.conda
       - conda: https://conda.anaconda.org/conda-forge/noarch/nodeenv-1.9.1-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/numba-0.60.0-py312hcccf92d_0.conda
-      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.1.8-pyhd8ed1ab_0.conda
+      - conda: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.2.0-pyhd8ed1ab_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/numcodecs-0.12.1-py312h275cf98_1.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/numpy-1.26.4-py312h8753938_0.conda
       - conda: https://conda.anaconda.org/conda-forge/win-64/openjpeg-2.5.2-h3d672ee_0.conda
@@ -24663,13 +24663,13 @@ packages:
   timestamp: 1718888686982
 - kind: conda
   name: numba_celltree
-  version: 0.1.8
+  version: 0.2.0
   build: pyhd8ed1ab_0
   subdir: noarch
   noarch: python
-  url: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.1.8-pyhd8ed1ab_0.conda
-  sha256: 4b30d964bf034b41ce14842bf8fa3040b21878c0da1bdec15c5523ab785bc311
-  md5: 02b10d14b2e6693519804fe90d41c589
+  url: https://conda.anaconda.org/conda-forge/noarch/numba_celltree-0.2.0-pyhd8ed1ab_0.conda
+  sha256: 8ef116fc2af9d70afe8123cb4157f5efa55cfd042fd1ef36cad9aab65b36ca5a
+  md5: e2ed9d4ac5f28671045cd33b2269969a
   depends:
   - numba >=0.50
   - numpy
@@ -24678,8 +24678,8 @@ packages:
   license_family: MIT
   purls:
   - pkg:pypi/numba-celltree?source=conda-forge-mapping
-  size: 33524
-  timestamp: 1722763371483
+  size: 33566
+  timestamp: 1724401764094
 - kind: conda
   name: numcodecs
   version: 0.12.1
@@ -35614,9 +35614,9 @@ packages:
   timestamp: 1607292254607
 - kind: pypi
   name: xugrid
-  version: 0.11.0
+  version: 0.11.2
   path: .
-  sha256: 3ef164ee6ff38acd9f8c2bffc7a7bf5b1773cdd558a25de2fb8388e27203431f
+  sha256: fc26928e4cddec4bdbc1b7a4d23417aa2f6733526debb5626136185d9ed019ca
   requires_dist:
   - numba
   - numba-celltree
diff --git a/pyproject.toml b/pyproject.toml
index c2c68b1eb..4f1a75e90 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -81,7 +81,7 @@ geopandas = "*"
 mapbox_earcut = "*"
 matplotlib-base = "*"
 netcdf4 = "*"
-numba_celltree = ">=0.1.8"
+numba_celltree = ">=0.2.0"
 numpy = "<2.0"
 pip = "*"
 pooch = "*"

From 973068a7a1e2ea0e3f675a82ffeb01099e59cebc Mon Sep 17 00:00:00 2001
From: Huite Bootsma <huitebootsma@gmail.com>
Date: Fri, 23 Aug 2024 12:04:43 +0200
Subject: [PATCH 8/8] Another note in changelog

---
 docs/changelog.rst | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/docs/changelog.rst b/docs/changelog.rst
index f0ee42371..0f1bc18b9 100644
--- a/docs/changelog.rst
+++ b/docs/changelog.rst
@@ -18,6 +18,12 @@ Fixed
   beyond the centroids for structured to structured regridding, and it would
   give nearest results (equal to :class:`CentroidLocatorRegridder`) otherwise.
 
+Changed
+~~~~~~~
+
+- Selection operations such as :meth:`UgridDataArrayAccessor.sel_points` will
+  now also return points that are located on the edges of 2D topologies.
+
 [0.11.2] 2024-08-16
 -------------------