Internal boundaries polygon is actually deletion polygons before mesh2d_set()

[veenstrajelmer]

Describe the bug
When deleting parts of the grid, holes are now represented in the face-node-connectivity. In order to preserve these, an array of polygons around deleted cells is maintained internally. When calling mk.mesh2d_get_mesh_inner_boundaries_as_polygons(), these polygons are exposed for backwards-compatibility reasons. However, the inner bounds should actually contain the inner bounds, not these internal polygons around individual deleted cells. When retrieving all bounds (which also include the inner bounds), the inner boundaries are correct.

Furthermore, if we read+write the grid to a netcdf file (or use mesh2d_set), the inner boundaries will be correct because mesh2d set retrieves the inner boundary from the actual face-node-connectivity.

To Reproduce

import numpy as np
from meshkernel import (ProjectionType, MakeGridParameters, MeshKernel, 
                        MeshRefinementParameters, GriddedSamples, RefinementType,
                        GeometryList, DeleteMeshOption,
                        )
import matplotlib.pyplot as plt
plt.close("all")

def get_illegal_xy(mk, all_bounds=False):
    # derive illegalcells as geodataframe
    if all_bounds:
        illegalcells_geom = mk.mesh2d_get_mesh_boundaries_as_polygons()
    else:
        illegalcells_geom = mk.mesh2d_get_mesh_inner_boundaries_as_polygons()
    illegal_xx = illegalcells_geom.x_coordinates
    illegal_yy = illegalcells_geom.y_coordinates
    illegal_xx[illegal_xx==geometry_separator] = np.nan
    illegal_yy[illegal_xx==geometry_separator] = np.nan
    return illegal_xx, illegal_yy


lon_min, lon_max, lat_min, lat_max = 147.75, 147.9, -40.4, -40.25
dxy = 0.05 #0.05

# grid generation and refinement with GEBCO bathymetry
# create base grid
projection = ProjectionType.SPHERICAL
make_grid_parameters = MakeGridParameters(angle=0,
                                          origin_x=lon_min,
                                          origin_y=lat_min,
                                          upper_right_x=lon_max,
                                          upper_right_y=lat_max,
                                          block_size_x=dxy,
                                          block_size_y=dxy)

mk = MeshKernel(projection=projection)
mk.curvilinear_compute_rectangular_grid_on_extension(make_grid_parameters)
mk.curvilinear_convert_to_mesh2d() #convert to ugrid/mesh2d

# refine
min_edge_size = 500 #in meters
gebco_elev = np.array([[-37, -37, -36, -36],
        [-41, -39, -38, -34],
        [-42, -32, -27,  -6],
        [-38,  -6,  -2,   2],
        [-43, -42, -31, -20]],dtype=np.int16)
lat_np = np.array([-40.397917, -40.364583, -40.33125 , -40.297917, -40.264583])
lon_np = np.array([147.76875 , 147.802083, 147.835417, 147.86875 ])
values_np = gebco_elev.flatten()
gridded_samples = GriddedSamples(x_coordinates=lon_np,y_coordinates=lat_np,values=values_np)

#refinement
mesh_refinement_parameters = MeshRefinementParameters(min_edge_size=min_edge_size, #always in meters
                                                      refinement_type=RefinementType(1), #Wavecourant/1,
                                                      connect_hanging_nodes=True, #set to False to do multiple refinement steps (e.g. for multiple regions)
                                                      smoothing_iterations=2,
                                                      max_courant_time=120)

mk.mesh2d_refine_based_on_gridded_samples(gridded_samples=gridded_samples,
                                          mesh_refinement_params=mesh_refinement_parameters,
                                          )    
# cutcells
geometry_separator = -999
xx1, yy1 = np.array([[147.77326613, -40.30182344],
       [147.79012097, -40.30091608],
       [147.78967742, -40.32541485],
       [147.77237903, -40.32586853],
       [147.77326613, -40.30182344]]).T

xx2, yy2 = np.array([[147.89701613, -40.34038632],
       [147.95068548, -40.33675687],
       [147.95245968, -40.40208693],
       [147.89790323, -40.40299429],
       [147.89701613, -40.34038632]]).T

xx3 = np.array([147.76125 , 147.764556, 147.780833, 147.802528, 147.829139,
       147.836222, 147.810861, 147.774111, 147.78125 , 147.76125 ])
yy3 = np.array([-40.253028, -40.249667, -40.241778, -40.25425 , -40.239222,
       -40.249639, -40.274278, -40.272611, -40.257222, -40.253028])

xx4 = np.array([147.83625 , 147.839556, 147.855833, 147.877528, 147.904139,
        147.911222, 147.885861, 147.849111, 147.85625 , 147.83625 ])
yy4 = np.array([-40.305028, -40.301667, -40.293778, -40.30625 , -40.291222,
        -40.301639, -40.326278-0.05, -40.324611-0.05, -40.309222, -40.305028]) + 0.052
# update yy to still generate two illegalcells after changing 
xx_list = [xx1, xx2, xx3, xx4]
yy_list = [yy1, yy2, yy3, yy4]
for xx, yy in zip(xx_list, yy_list):
    delete_pol_geom = GeometryList(x_coordinates=xx, y_coordinates=yy, geometry_separator=geometry_separator)
    mk.mesh2d_delete(geometry_list=delete_pol_geom, 
                     delete_option=DeleteMeshOption.INSIDE_NOT_INTERSECTED,
                     invert_deletion=False)

illegal1_xx, illegal1_yy = get_illegal_xy(mk)
bounds1_xx, bounds1_yy = get_illegal_xy(mk, all_bounds=True)

fig, ax = plt.subplots()
mk.mesh2d_get().plot_edges(ax=ax)
ax.plot(illegal1_xx, illegal1_yy, "r-")
ax.plot(bounds1_xx, bounds1_yy, "k--")
ax.plot(xx1, yy1, "c-")
ax.plot(xx2, yy2, "c-")
ax.plot(xx3, yy3, "c-")
ax.plot(xx4, yy4, "c-")

# create a fresh mk instance with the same grid
mk2 = MeshKernel(projection=mk.get_projection())
mk2.mesh2d_set(mk.mesh2d_get())

illegal2_xx, illegal2_yy = get_illegal_xy(mk2)
bounds2_xx, bounds2_yy = get_illegal_xy(mk2, all_bounds=True)

fig, ax = plt.subplots()
mk2.mesh2d_get().plot_edges(ax=ax)
ax.plot(illegal2_xx, illegal2_yy, "r-")
ax.plot(bounds2_xx, bounds2_yy, "k--")

Results in the following figure (mesh in blue, all bounds in black, innerbounds in red, deletion polygon in cyan):

Versus the same grid after applying mesh2d_set():

Expected behavior
The figures should be identical: the inner boundaries between the meshes are not changed and should not consist of a collection of polygons around deleted cells.

Version info:

• OS: Windows
• Version 8.2.1

Additional information
To discuss:

• Originally in interacter, illegalcells only contained polygons that would result in high ortho (so max 6 edges).
• What does this issue mean for preservation of deleted cells that are manually being re-added?
• The current implementation of the internal polygons maintains polygons for all the deleted cells (even cells outside of the outer boundary of the grid), that are looped over in deletion of faces after each re-administration of the topology. Can this not be very troublesome for performance when having deleted a lot of cells (say ~1million)?
• if the polygons around deleted cells are being used to delete cells from the readministrated grid, the new centroids might be exactly on the original cell edges, causing them not to be deleted (with some other arm or float behaviour)