inconsistency between number of edges and number of orthogonality values

[veenstrajelmer]

Describe the bug
The number of values in the orthogonality array is not the same as the number of edges, while this should be the case.

To Reproduce

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

# input params
lon_min, lon_max, lat_min, lat_max = 147.75, 147.9, -40.4, -40.25
dxy = 0.05 #0.05
crs = 'EPSG:4326'

# 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
print('Refining basegrid')
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,
                                           use_nodal_refinement=True) #TODO: what does this do?

print("max ortho after refining", mk.mesh2d_get_orthogonality().values.max()) #TODO: couple back to uds, currently ordering mismatch: https://github.com/Deltares/MeshKernelPy/issues/72

#cutcells
print('Cutting cells')
xx = np.array([147.83625 , 147.839556, 147.855833, 147.877528, 147.904139,
       147.911222, 147.885861, 147.849111, 147.85625 , 147.83625 ])
yy = np.array([-40.305028, -40.301667, -40.293778, -40.30625 , -40.291222,
       -40.301639, -40.326278, -40.324611, -40.309222, -40.305028])
delete_pol_geom = GeometryList(x_coordinates=xx, y_coordinates=yy) #TODO: .copy()/to_numpy() makes the array contiguous in memory, which is necessary for meshkernel.mesh2d_delete()
mk.mesh2d_delete(geometry_list=delete_pol_geom, 
                 delete_option=DeleteMeshOption.INSIDE_NOT_INTERSECTED,
                 invert_deletion=False)

num_edges = len(mk.mesh2d_get().edge_x)
num_orthovalues = len(mk.mesh2d_get_orthogonality().values)
print("num_edges:", num_edges)
print("num_orthovalues:", num_orthovalues)
assert num_edges == num_orthovalues

Raises AssertionError:

num_edges:  383
num_orthovalues:  384
Traceback (most recent call last):

  File ~\Anaconda3\envs\dfm_tools_env\Lib\site-packages\spyder_kernels\py3compat.py:356 in compat_exec
    exec(code, globals, locals)

  File c:\data\dfm_tools\dfm_tools\untitled1.py:79
    assert num_edges == num_orthovalues

AssertionError

Expected behavior
The number of values in the orthogonality array should be equal to the number of edges.

Version info (please complete the following information):

• OS: Windows
• Version meshkernelpy wheel from 17-05-2024