add contacts testcase that fails in hydrolib-core on macos 14 arm64

[veenstrajelmer]

What is the need for this task.
In Deltares/HYDROLIB-core#635 we found that HYDROLIB-core supports macos-12 and macos-13 on github. However, there is one failing testcase on macos-14. This is the only arm64 architecture, 12 and 13 are x86_64. This testcase uses meshkernel in the background. Also, the code was converted to meshkernelpy-only code and that also failed on the assertion. Therefore, it seems to point at a usecase that the meshkernelpy testbench does not cover yet.

What is the task?
Add this as a testcase to the meshkernelpy testbench. The meshkernel-only code is this:

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

def get_circle_gl(r, detail=100):
    t = np.r_[np.linspace(0, 2 * np.pi, detail), 0]
    polygon = GeometryList(np.cos(t) * r, np.sin(t) * r)
    return polygon


# Define line (spiral)
theta = np.arange(0.1, 20, 0.01)

y = np.sin(theta) * theta
x = np.cos(theta) * theta

dists = np.r_[0.0, np.cumsum(np.hypot(np.diff(x), np.diff(y)))]
dists = dists[np.arange(0, len(dists), 20)]

# Create Mesh1d
mki = MeshKernel()

# def mesh1d_add_branch(mki, branch)
mesh1d_node_x = np.array([  0.09950042,   0.28660095,   0.43879128,   0.53538953,
         0.55944897,   0.49895573,   0.34774848,   0.1061058 ,
        -0.21903564,  -0.61425018,  -1.06017682,  -1.53243485,
        -2.00285904,  -2.44099478,  -2.81577868,  -3.09731897,
        -3.25868324,  -3.27759841,  -3.13797012,  -2.83113599,
        -2.35677818,  -1.72343644,  -0.9485811 ,  -0.05822672,
         0.91391061,   1.92768649,   2.93818398,   3.89768376,
         4.75786287,   5.47212274,   5.99793747,   6.29910941,
         6.34781957,   6.12636709,   5.62850319,   4.86028133,
         3.84036588,   2.59976488,   1.18097874,  -0.36341679,
        -1.97270765,  -3.58042781,  -5.11710117,  -6.51322582,
        -7.70237336,  -8.62426658,  -9.22769553,  -9.47313548,
        -9.33493933,  -8.80299241,  -7.88373862,  -6.6005121 ,
        -4.99313774,  -3.11679531,  -1.04017448,   1.1570199 ,
         3.38712238,   5.55800473,   7.57687716,   9.35423652,
        10.80779395,  11.8662112 ,  12.47247849,  12.58677787,
        12.18869449,  11.27866268,   9.878564  ,   8.03142895,
         5.80023126,   3.26580248,   0.52393323,  -2.31823644,
        -5.14640187,  -7.84369048, -10.29548549, -12.39427773,
       -14.04434094, -15.16602867, -15.69950221, -15.60771894,
       -14.87853808, -13.52583451, -11.58955145,  -9.13466558,
        -6.24908366,  -3.04053512,   0.36743138,   3.84019936,
         7.23740149,  10.41859212,  13.24904611,  15.6054449 ,
        17.38121053,  18.49125831,  18.87595858,  18.50412697,
        17.3748994 ,  15.51839191,  12.99509423,   9.89399732])
mesh1d_node_y = np.array([ 9.98334166e-03,  8.86560620e-02,  2.39712769e-01,  4.50952381e-01,
        7.04994219e-01,  9.80328096e-01,  1.25262564e+00,  1.49624248e+00,
        1.68583018e+00,  1.79797017e+00,  1.81273967e+00,  1.71512199e+00,
        1.49618036e+00,  1.15392568e+00,  6.93823055e-01,  1.28900054e-01,
       -5.20560791e-01, -1.22774130e+00, -1.96039372e+00, -2.68228802e+00,
       -3.35493616e+00, -3.93951353e+00, -4.39888553e+00, -4.69963931e+00,
       -4.81401780e+00, -4.72165488e+00, -4.41101744e+00, -3.88047179e+00,
       -3.13890759e+00, -2.20587232e+00, -1.11119128e+00,  1.05927573e-01,
        1.39827992e+00,  2.71249447e+00,  3.99123437e+00,  5.17568018e+00,
        6.20818733e+00,  7.03499983e+00,  7.60889540e+00,  7.89163660e+00,
        7.85610747e+00,  7.48802622e+00,  6.78714046e+00,  5.76783230e+00,
        4.45908562e+00,  2.90379510e+00,  1.15742614e+00, -7.13935644e-01,
       -2.63607808e+00, -4.52960535e+00, -6.31321355e+00, -7.90716384e+00,
       -9.23680548e+00, -1.02359947e+01, -1.08502552e+01, -1.10395337e+01,
       -1.07804175e+01, -1.00677000e+01, -8.91520793e+00, -7.35583164e+00,
       -5.44073432e+00, -3.23775103e+00, -8.29023717e-01,  1.69204693e+00,
        4.22442026e+00,  6.66346518e+00,  8.90527784e+00,  1.08510898e+01,
        1.24115800e+01,  1.35109043e+01,  1.40902624e+01,  1.41108391e+01,
        1.35559783e+01,  1.24324784e+01,  1.07709321e+01,  8.62507273e+00,
        6.07013077e+00,  3.20024597e+00,  1.25021985e-01, -3.03465144e+00,
       -6.15134982e+00, -9.09625202e+00, -1.17444581e+01, -1.39802677e+01,
       -1.57021958e+01, -1.68275116e+01, -1.72960977e+01, -1.70734551e+01,
       -1.61527094e+01, -1.45555123e+01, -1.23317792e+01, -9.55824719e+00,
       -6.33589144e+00, -2.78628178e+00,  9.52988800e-01,  4.73363337e+00,
        8.40255146e+00,  1.18080275e+01,  1.48059963e+01,  1.72661176e+01])
nlinks = len(mesh1d_node_y)-1
new_edge_nodes = np.stack([np.arange(nlinks), np.arange(nlinks) + 1], axis=1)
m1d = Mesh1d(node_x=mesh1d_node_x.astype(np.float64),
                        node_y=mesh1d_node_y.astype(np.float64),
                        edge_nodes=new_edge_nodes.ravel().astype(np.int32),
                        )
mki.mesh1d_set(m1d)

# Add Mesh2d
xmin, ymin, xmax, ymax = (-22, -22, 22, 22)
dx=2
dy=2
rows = int((ymax - ymin) / dy)
columns = int((xmax - xmin) / dx)
params = MakeGridParameters(
    num_columns=columns,
    num_rows=rows,
    origin_x=xmin,
    origin_y=ymin,
    block_size_x=dx,
    block_size_y=dy,
)
mki.curvilinear_compute_rectangular_grid(params)
mki.curvilinear_convert_to_mesh2d()  # convert to ugrid/mesh2d

# clip mesh inner circle
mki.mesh2d_delete(
    geometry_list=get_circle_gl(22),
    delete_option=DeleteMeshOption.INSIDE_NOT_INTERSECTED,
    invert_deletion=False,
)

# Add links
mki.contacts_compute_single(
    node_mask=np.full(nlinks+1, True, dtype=bool),
    polygons=get_circle_gl(19), 
    projection_factor=1.0
)

fig, ax = plt.subplots()
mesh2d_output = mki.mesh2d_get()
mesh1d_output = mki.mesh1d_get()
links_output = mki.contacts_get()
mesh2d_output.plot_edges(ax=ax, color="r")
mesh1d_output.plot_edges(ax=ax, color="g")
links_output.plot_edges(
    ax=ax, mesh1d=mesh1d_output, mesh2d=mesh2d_output, color="k"
)

contacts = mki.contacts_get()
network1_link1d2d = np.stack([contacts.mesh1d_indices, contacts.mesh2d_indices], axis=1)
assert network1_link1d2d.shape == (21, 2)

This tests also only fails on macos-14 (on Github) on the line assert network1_link1d2d.shape == (21, 2) and with a comparable error message: assert (19, 2) == (21, 2) (19 instead of 20 as in the original issue). This 19 should be 21 just like on other platforms/architecture.

The code also produces this figure (on windows, so this is the reference/correct figure):

[backeb]

Running the above example on MacOS 14 arm64 produces the following output

---------------------------------------------------------------------------
AssertionError                            Traceback (most recent call last)
Cell In[1], [line 130](vscode-notebook-cell:?execution_count=1&line=130)
    [128](vscode-notebook-cell:?execution_count=1&line=128) contacts = mki.contacts_get()
    [129](vscode-notebook-cell:?execution_count=1&line=129) network1_link1d2d = np.stack([contacts.mesh1d_indices, contacts.mesh2d_indices], axis=1)
--> [130](vscode-notebook-cell:?execution_count=1&line=130) assert network1_link1d2d.shape == (21, 2)

AssertionError:

sw_vers output

ProductName:            macOS
ProductVersion:         14.4.1
BuildVersion:           23E224