added testcase and replaced Link1d2d._process() with properties

hydrolib/core/dflowfm/net/models.py

@@ -589,14 +589,29 @@ class Link1d2d(BaseModel):
     """
 
     meshkernel: mk.MeshKernel = Field(default_factory=mk.MeshKernel)
-
-    link1d2d_id: np.ndarray = Field(default_factory=lambda: np.empty(0, object))
-    link1d2d_long_name: np.ndarray = Field(default_factory=lambda: np.empty(0, object))
-    link1d2d_contact_type: np.ndarray = Field(
-        default_factory=lambda: np.empty(0, np.int32)
-    )
-    link1d2d: np.ndarray = Field(default_factory=lambda: np.empty((0, 2), np.int32))
-
+
+    @property
+    def link1d2d(self) -> np.ndarray[np.int32]:
+        contacts = self.meshkernel.contacts_get()
+        link1d2d_arr = np.stack([contacts.mesh1d_indices, contacts.mesh2d_indices], axis=1)
+        return link1d2d_arr
+
+    @property
+    def link1d2d_contact_type(self) -> np.ndarray[np.int32]:
+        contacts = self.meshkernel.contacts_get()
+        link1d2d_contact_type_arr = np.full(contacts.mesh1d_indices.size, 3)
+        return link1d2d_contact_type_arr
+
+    @property
+    def link1d2d_id(self) -> np.ndarray[object]:
+        link1d2d_id_arr = np.array([f"{n1d:d}_{f2d:d}" for n1d, f2d in self.link1d2d])
+        return link1d2d_id_arr
+
+    @property
+    def link1d2d_long_name(self) -> np.ndarray[object]:
+        link1d2d_id_arr = np.array([f"{n1d:d}_{f2d:d}" for n1d, f2d in self.link1d2d])
+        return link1d2d_id_arr
+
     def is_empty(self) -> bool:
         """Whether this Link1d2d is currently empty.
 
@@ -626,23 +641,6 @@ def clear(self) -> None:
         self.meshkernel._allocate_state(self.meshkernel.get_projection())
         self.meshkernel.contacts_get()
 
-    def _process(self) -> None:
-        """
-        Get links from meshkernel and add to the array with link administration
-        """
-        contacts = self.meshkernel.contacts_get()
-
-        self.link1d2d = np.append(
-            self.link1d2d,
-            np.stack([contacts.mesh1d_indices, contacts.mesh2d_indices], axis=1),
-            axis=0,
-        )
-        self.link1d2d_contact_type = np.append(
-            self.link1d2d_contact_type, np.full(contacts.mesh1d_indices.size, 3)
-        )
-        self.link1d2d_id = np.array([f"{n1d:d}_{f2d:d}" for n1d, f2d in self.link1d2d])
-        self.link1d2d_long_name = np.array([f"{n1d:d}_{f2d:d}" for n1d, f2d in self.link1d2d])
-
     def _link_from_1d_to_2d(
         self, node_mask: np.ndarray, polygon: mk.GeometryList = None
     ):
@@ -662,8 +660,7 @@ def _link_from_1d_to_2d(
         self.meshkernel.contacts_compute_single(
             node_mask=node_mask, polygons=polygon, projection_factor=1.0
         )
-        self._process()
-
+
         # Note that the function "contacts_compute_multiple" also computes the connections, but does not take into account
         # a bounding polygon or the end points of the 1d mesh.
 
@@ -672,8 +669,7 @@ def _link_from_2d_to_1d_embedded(
     ):
         """"""
         self.meshkernel.contacts_compute_with_points(node_mask=node_mask, points=points)
-        self._process()
-
+
     def _link_from_2d_to_1d_lateral(
         self,
         node_mask: np.ndarray,
@@ -687,7 +683,6 @@ def _link_from_2d_to_1d_lateral(
         self.meshkernel.contacts_compute_boundary(
             node_mask=node_mask, polygons=polygon, search_radius=search_radius
         )
-        self._process()
 
 
 class Mesh1d(BaseModel):

hydrolib/core/dflowfm/net/reader.py

@@ -121,12 +121,13 @@ def read_link1d2d(self, link1d2d: Link1d2d) -> None:
         ds = nc.Dataset(self._ncfile_path)  # type: ignore[import]
 
         # Read mesh1d
-        for meshkey, nckey in self._explorer.link1d2d_var_name_mapping.items():
-            setattr(link1d2d, meshkey, self._read_nc_attribute(ds[nckey]))
-
+        link1d2d_contact_type_arr = self._read_nc_attribute(ds['link1d2d_contact_type'])
+        assert (link1d2d_contact_type_arr == 3).all()
+
+        link1d2d_arr = self._read_nc_attribute(ds['link1d2d'])
         # set contacts on meshkernel, use .copy() to avoid strided arrays
-        mesh1d_indices = link1d2d.link1d2d[:, 0].copy()
-        mesh2d_indices = link1d2d.link1d2d[:, 1].copy()
+        mesh1d_indices = link1d2d_arr[:, 0].copy()
+        mesh2d_indices = link1d2d_arr[:, 1].copy()
         contacts = Contacts(
             mesh1d_indices=mesh1d_indices, mesh2d_indices=mesh2d_indices
         )

hydrolib/core/dflowfm/net/writer.py

@@ -330,7 +330,6 @@ def _set_1dmesh(self, ncfile: nc.Dataset, mesh1d: Mesh1d) -> None:  # type: igno
         mesh1d_node_offset.units = "m"
         mesh1d_node_offset[:] = mesh1d.mesh1d_node_branch_offset
 
-
         mesh_edge_x = ncfile.createVariable("mesh1d_edge_x", np.float64, "mesh1d_nEdges")
         mesh_edge_x.standard_name = "projection_x_coordinate"
         mesh_edge_x.long_name = "Characteristic x-coordinate of the mesh edge (e.g. midpoint)"

tests/dflowfm/test_net.py

@@ -100,7 +100,6 @@ def get_circle_gl(r, detail=100):
 
 @pytest.mark.plots
 def test_create_1d_2d_1d2d():
-
     # Define line (spiral)
     theta = np.arange(0.1, 20, 0.01)
 
@@ -171,8 +170,53 @@ def test_create_1d_2d_1d2d():
     network2.plot(ax=ax2)
 
 
-def test_create_2d():
+def test_create_1d_2d_1d2d_call_link_generation_twice():
+    """
+    a second generation of links should overwrite the existing links,
+    this testcase checks whether that is the case.
+    Related issue: https://github.com/Deltares/HYDROLIB-core/issues/546
+    """
+    # 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 branch
+    branch = Branch(geometry=np.stack([x, y], axis=1), branch_offsets=dists)
+
+    # Create Mesh1d
+    network = Network()
+    network.mesh1d_add_branch(branch, name="branch1")
+
+    branch = Branch(geometry=np.array([[-25.0, 0.0], [x[0], y[0]]]))
+    branch.generate_nodes(mesh1d_edge_length=2.5)
+    network.mesh1d_add_branch(branch, name="branch2")
+
+    # Add Mesh2d
+    network.mesh2d_create_rectilinear_within_extent(
+        extent=(-22, -22, 22, 22), dx=2, dy=2
+    )
+    network.mesh2d_clip_mesh(geometrylist=get_circle_gl(22))
 
+    network.mesh2d_refine_mesh(polygon=get_circle_gl(11), level=1)
+    network.mesh2d_refine_mesh(polygon=get_circle_gl(3), level=1)
+
+    # Add links, do this twice to check if contacts are overwritten and not appended
+    network.link1d2d_from_1d_to_2d(branchids=["branch1"], polygon=get_circle_gl(19))
+    network.link1d2d_from_1d_to_2d(branchids=["branch1"], polygon=get_circle_gl(19))
+
+    network1_link1d2d = network._link1d2d.link1d2d
+    assert network1_link1d2d.shape == (21, 2)
+    assert len(network._link1d2d.link1d2d_contact_type) == 21
+    assert len(network._link1d2d.link1d2d_id) == 21
+    assert len(network._link1d2d.link1d2d_long_name) == 21
+
+
+def test_create_2d():
     # Define polygon
     bbox = (1.0, -2.0, 3.0, 4.0)