curvilinear_make_uniform_on_extension does not accept a geometry_list

[veenstrajelmer]

Describe the bug

• curvilinear_make_uniform accepts a geometry_list, but gives unexpected result (figure below)
• curvilinear_make_uniform_on_extension does not accept a geometry_list, it raises "TypeError: curvilinear_make_uniform_on_extension() takes 2 positional arguments but 3 were given"
• when using is_geographic=True with curvilinear_make_uniform and a geometry_list argument (creating spherical grid in polygon) it raises: "MeshKernelError: bad allocation"

To Reproduce

import meshkernel
import matplotlib.pyplot as plt
plt.close('all')
import numpy as np

on_extension = False
grid_in_pol = False #grid_in_pol=True with on_extension=False gives non-covering grid, grid_in_pol=True/False with on_extension=True raises "TypeError: curvilinear_make_uniform_on_extension() takes 2 positional arguments but 3 were given"
is_geographic = False #False/True/True: creating spherical grid in pol raises: "MeshKernelError: bad allocation"

#general settings
lon_min,lon_max = -6,2
lat_min,lat_max = 48.5,51.2
lon_res,lat_res = 0.2,0.2

# Create an instance of MakeGridParameters and set the values
make_grid_parameters = meshkernel.MakeGridParameters()
make_grid_parameters.angle = 0.0
make_grid_parameters.origin_x = lon_min
make_grid_parameters.origin_y = lat_min
if on_extension:
    make_grid_parameters.upper_right_x = lon_max
    make_grid_parameters.upper_right_y = lat_max
else:
    num_x = int(np.ceil((lon_max-lon_min)/lon_res))
    num_y = int(np.ceil((lat_max-lat_min)/lat_res))
    if is_geographic:
        num_y = num_y*2 #TODO: remove *2, necessary to get correct lat grid extent with is_geographic=True
    make_grid_parameters.num_columns = num_x
    make_grid_parameters.num_rows = num_y
make_grid_parameters.block_size_x = lon_res
make_grid_parameters.block_size_y = lat_res


# If a polygon is provided it will be used in the generation of the curvilinear grid. The polygon must be closed
if grid_in_pol: #can be used instead of origin_x/origin_y and num_x/num_y
    pol_x = np.array([-6,-4,0,-6], dtype=np.double)
    pol_y = np.array([48,51,49.5,48], dtype=np.double)
    geometry_list = meshkernel.GeometryList(pol_x, pol_y)
else:
    geometry_list = None

mk2 = meshkernel.MeshKernel(is_geographic=is_geographic)
if on_extension:
    mk2.curvilinear_make_uniform_on_extension(make_grid_parameters,geometry_list)
else:
    mk2.curvilinear_make_uniform(make_grid_parameters,geometry_list)
mk2.curvilinear_convert_to_mesh2d() #convert to ugrid/mesh2d

mesh2d = mk2.mesh2d_get()
fig, ax = plt.subplots()
if grid_in_pol:
    ax.plot(pol_x,pol_y,'r-')
mesh2d.plot_edges(ax)
print(mk2.mesh2d_get_orthogonality().values.max())

Expected behavior
...

Version info (please complete the following information):

• OS: [e.g. Windows]
• Version branch 2.1.0 wheel of 4-7-2023

[veenstrajelmer]

With the meshkernel 3.0.0 pre-release the following code results in the desired behaviour in all cases:

import meshkernel
import matplotlib.pyplot as plt
plt.close('all')
import numpy as np

on_extension = True
grid_in_pol = False
projection = False

if on_extension and grid_in_pol:
    raise ValueError('on_extension and grid_in_pol cannot be True at the same time')

#general settings
lon_min,lon_max = -6,2
lat_min,lat_max = 48.5,51.2
lon_res,lat_res = 0.2,0.2

# Create an instance of MakeGridParameters and set the values
make_grid_parameters = meshkernel.MakeGridParameters()
make_grid_parameters.angle = 0.0
make_grid_parameters.origin_x = lon_min
make_grid_parameters.origin_y = lat_min
if on_extension:
    make_grid_parameters.upper_right_x = lon_max
    make_grid_parameters.upper_right_y = lat_max
else:
    num_x = int(np.ceil((lon_max-lon_min)/lon_res))
    num_y = int(np.ceil((lat_max-lat_min)/lat_res))
    if projection:
        num_y = num_y*2 #TODO: remove *2, necessary to get correct lat grid extent with is_geographic=True
    make_grid_parameters.num_columns = num_x
    make_grid_parameters.num_rows = num_y
make_grid_parameters.block_size_x = lon_res
make_grid_parameters.block_size_y = lat_res


# If a polygon is provided it will be used in the generation of the curvilinear grid. The polygon must be closed
if grid_in_pol: #can be used instead of origin_x/origin_y and num_x/num_y
    pol_x = np.array([-6,-4,0,-6], dtype=np.double)
    pol_y = np.array([48,51,49.5,48], dtype=np.double)
    geometry_list = meshkernel.GeometryList(pol_x, pol_y)
else:
    geometry_list = None

mk2 = meshkernel.MeshKernel(projection=projection)
if on_extension:
    mk2.curvilinear_compute_rectangular_grid_on_extension(make_grid_parameters)
elif geometry_list is None:
    mk2.curvilinear_compute_rectangular_grid(make_grid_parameters)
else:
    mk2.curvilinear_compute_rectangular_grid_from_polygon(make_grid_parameters,geometry_list)
mk2.curvilinear_convert_to_mesh2d() #convert to ugrid/mesh2d

mesh2d = mk2.mesh2d_get()
fig, ax = plt.subplots()
if grid_in_pol:
    ax.plot(pol_x,pol_y,'r-')
mesh2d.plot_edges(ax)
print(mk2.mesh2d_get_orthogonality().values.max())