add contacts testcase

tests/test_mesh1d_basics.py

@@ -1,10 +1,18 @@
+import matplotlib.pyplot as plt
 import numpy as np
 import pytest
 from mesh2d_factory import Mesh2dFactory
 from numpy import ndarray
 from numpy.testing import assert_array_equal
 
-from meshkernel import Contacts, DeleteMeshOption, GeometryList, Mesh1d, MeshKernel
+from meshkernel import (
+    Contacts,
+    DeleteMeshOption,
+    GeometryList,
+    MakeGridParameters,
+    Mesh1d,
+    MeshKernel,
+)
 
 
 def sort_contacts_by_mesh2d_indices(contacts):
@@ -574,3 +582,288 @@ def test_contacts_compute_with_points_after_deletion(
 
     assert contacts.mesh2d_indices[0] == 0
     assert contacts.mesh2d_indices[1] == 7
+
+
+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
+
+
+def test_contacts_compute_single_circle():
+    # 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.25262564e00,
+            1.49624248e00,
+            1.68583018e00,
+            1.79797017e00,
+            1.81273967e00,
+            1.71512199e00,
+            1.49618036e00,
+            1.15392568e00,
+            6.93823055e-01,
+            1.28900054e-01,
+            -5.20560791e-01,
+            -1.22774130e00,
+            -1.96039372e00,
+            -2.68228802e00,
+            -3.35493616e00,
+            -3.93951353e00,
+            -4.39888553e00,
+            -4.69963931e00,
+            -4.81401780e00,
+            -4.72165488e00,
+            -4.41101744e00,
+            -3.88047179e00,
+            -3.13890759e00,
+            -2.20587232e00,
+            -1.11119128e00,
+            1.05927573e-01,
+            1.39827992e00,
+            2.71249447e00,
+            3.99123437e00,
+            5.17568018e00,
+            6.20818733e00,
+            7.03499983e00,
+            7.60889540e00,
+            7.89163660e00,
+            7.85610747e00,
+            7.48802622e00,
+            6.78714046e00,
+            5.76783230e00,
+            4.45908562e00,
+            2.90379510e00,
+            1.15742614e00,
+            -7.13935644e-01,
+            -2.63607808e00,
+            -4.52960535e00,
+            -6.31321355e00,
+            -7.90716384e00,
+            -9.23680548e00,
+            -1.02359947e01,
+            -1.08502552e01,
+            -1.10395337e01,
+            -1.07804175e01,
+            -1.00677000e01,
+            -8.91520793e00,
+            -7.35583164e00,
+            -5.44073432e00,
+            -3.23775103e00,
+            -8.29023717e-01,
+            1.69204693e00,
+            4.22442026e00,
+            6.66346518e00,
+            8.90527784e00,
+            1.08510898e01,
+            1.24115800e01,
+            1.35109043e01,
+            1.40902624e01,
+            1.41108391e01,
+            1.35559783e01,
+            1.24324784e01,
+            1.07709321e01,
+            8.62507273e00,
+            6.07013077e00,
+            3.20024597e00,
+            1.25021985e-01,
+            -3.03465144e00,
+            -6.15134982e00,
+            -9.09625202e00,
+            -1.17444581e01,
+            -1.39802677e01,
+            -1.57021958e01,
+            -1.68275116e01,
+            -1.72960977e01,
+            -1.70734551e01,
+            -1.61527094e01,
+            -1.45555123e01,
+            -1.23317792e01,
+            -9.55824719e00,
+            -6.33589144e00,
+            -2.78628178e00,
+            9.52988800e-01,
+            4.73363337e00,
+            8.40255146e00,
+            1.18080275e01,
+            1.48059963e01,
+            1.72661176e01,
+        ]
+    )
+    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)