use face_coordinates from mesh_topology variable

[veenstrajelmer]

face_coordinates are not read from mesh_topology variable, but computed instead. There is a discrepancy between the center of mass and circumcenter of triangles here. Solve by using face_coordinates if available.

face_coords_manual = np.c_[uds_test.mesh2d_face_x.to_numpy(),uds_test.mesh2d_face_y.to_numpy()]
face_coords_uds = uds.grid.face_coordinates
print((face_coords_manual == uds_test.grid.face_coordinates).all())
Out:
False

[veenstrajelmer]

from https://content.oss.deltares.nl/delft3d/D-Flow_FM_User_Manual.pdf

"D-Flow FM utilizes the circumcenter as the basis of the definition of the elementary flow variables ’water level’ and ’flow velocity’. The water level is defined at the circumcenter, whereas the face normal flow velocity is defined at the orthogonal projection of the circumcenter onto the cell face, i.e., the midpoint of the cell face."

xugrid computes the centroids, whereas FM writes the circumcenters in the face_coordinates variables.

[Huite]

Some thoughts:

• The centroid is a more "robust" coordinate, because it is guaranteed to be inside of a convex polygon. The circumcenter may be located outside.
• However, it's somewhat surprising that face_coordinates are overwritten silently during a roundtrip, if the coordinates are circumcenters.
• Best to clearly separate internally between the centroids and the circumcenters. Most methods and functions should exclusively and explicit call for centroids.
• When reading a file, we cannot know whether the face coordinates are one or the other. So we should recompute the centroids as needed, the first time.
• Centroids are stored in _centroids. They are written only if face_coordinates are not already present.

[Huite]

Reading this again: the problem is that a circumcenter is defined for triangles, but not for arbitrary n-gons.
This means that in any internal method, centroids should be used, because it is guaranteed to be available. All the regridding logic already does this, as far as I can tell. Centroids should be kept fully separate from face_coordinates.

The Ugrid2d class should probably have the face coordinates as an additional parameter, which is then maintained. When nothing exists, centroids are returned as the default option for face coordinates. The question then arises how to set specify the face coordinates as circumcenters. In principle, this could be made to work:

grid.face_coordinates = grid.circumcenters

But that does not provide easy access via an accessor for DataArrays or Datasets; that's probably best achieved via the assign_face_coords:

    def assign_face_coords(self, centroid: bool=True) -> UgridDataset:

Using circumcenters if false. This would then also set the face_coordinates on the underlying grid.