Give access to all neighbouring elements for any space domain element

[ThibHlln]

When the framework moves away from only supporting structured grids, an increment on X or Y will not work to find all neighbours of any given mesh element, so it needs to provide a functionality to query where the neighbours of any given mesh element are.

This is required when interpolating between two grids: if two cells are neighbours and it needs to be interpolated onto another grid, the interpolation scheme will need to know about this neighbouring relationship.

It would also be a useful facility to provide to model developers, e.g. when they need to route flow laterally (whether on the surface, or under the surface), they can find where the neighbours upstream and downstream are.

[ThibHlln]

For regular grids (LatLonGrid, RotatedLatLonGrid, BritishNationalGrid), the flow_direction attribute can provide some of this information.

For irregular grids, following the UGRID conventions https://ugrid-conventions.github.io, this information would be stored under the face_face_connectivity variable. However, I don't think the convention supports a directed connectivity. This may have to be inferred from topographic information on the faces?

[ThibHlln]

Related to the need for neighbouring relationships and directed connectivity, at the moment the flow_direction assumes that the connectivity is mono-directional, but there may be a need to consider multi-directional flow, e.g. if grid cell contains a ridge, there is a case for a fraction of the runoff towards two of more valleys.

E.g. for structured grids, groundwater models (MODFLOW) consider up to 8 neighbouring cells, and determine the ratio of flow going to each neighbour based on the hydraulic head gradient between a cell and its 8 neighbours. For surface runoff, some components may want to apply a similar approach, based on elevation gradients?

The flow_direction attribute may need to store a vector of 8 values for each grid cell containing weights to attribute fractions of runoff to each neighbour. A mono-directional case would still be supported by assigning all the weight to the given unique flow direction.