Add a refine_with_vertices method to Ugrid1d

[Huite]

After discussion with @JoerivanEngelen:

To translate data from networks, e.g. the water level output or bed level from 1D hydraulic simulations, to gridded (MODFLOW) input, we need some tooling that's capable of providing this data at the appropriate scale.

An appropriate method is to intersect the network topology with the grid topology, taking the midpoints of the resulting intersecting edges, and adding these as nodes to the original network. Then, we can apply some form of interpolation along the new network based on the connectivity to propagate features to arbitrary locations.

Some changes in numba_celltree will allow us to build a tree for 1D networks and search for edges efficiently.
Deltares/numba_celltree#11

Some implementation ideas. First we search for the new nodes, and remove far away vertices. Next, we need to create a new edge_node_connectivity. Since we're refining, we're just replacing original edges by refined smaller ones. The ones to refine can identified by their count.

edge_index = self.celltree.locate_points(vertices)
valid = edge_index != -1
edge_index = edge_index[valid]
new_vertices = vertices[valid]
n_edge = np.bincount(edge_index) + 1

We should probably also check for uniques of the old and new vertices.
Anyway, to set up the new edge node connectivity, we need to set the right order. E.g. 0 -> 1 may become 0 -> 3 -> 2 -> 1.
To find the appropriate order, we just need to compute the distance from the first vertex of the original edges.

distance = np.linalg.norm(vertices - edge_nodes[edge_index, 0])

Then we use a lexsort to first sort by edge_index, then distance. To make things a little easier, we add the original network vertices as well.

sorting_edge_index = np.concatenate((edge_index, self.edge_node_connectivity[:, 0], self.edge_node_connectivity[:, 1]))
distance = np.concatenate((distance, np.full(self.n_edge, 0.0), np.full(self.n_edge, FLOAT_MAX))
order = np.lexsort((distance, sorting_edge_index))

This puts everything in order. Now we just need to tie the nodes together to form the (sub-)edges.

node_index = np.concatenate((np.arange(self.n_node, self.n_node + len(edge_index), self.edge_node_connectivity[:, 0], self.edge_node_connectivity[:, 1]))[order]
new_edges = np.column_stack((node_index[:-1], node_index[1:]))

This potentially includes nodes that are tied "to themselves" (since they are the end of on edge, and the start of the next).
We can easily remove those and insert them into the original edge_node_connectivity to construct a new grid.

new_edges = new_edges[new_edges[:, 0] != new_edges[:, 1]]
n_edge = np.bincount(edge_index) + 1
edges = np.repeat(self.edge_node_connectivity, repeats=n_edge)
edges[n_edge > 1] = new_edges
new_grid = Ugrid1d(*new_vertices.T, self.fill_value, edges)