diff --git a/docs/api.rst b/docs/api.rst index c6a7924a..05598896 100644 --- a/docs/api.rst +++ b/docs/api.rst @@ -218,6 +218,8 @@ UGRID1D Topology Ugrid1d.node_edge_connectivity Ugrid1d.node_node_connectivity + Ugrid1d.directed_node_node_connectivity + Ugrid1d.directed_edge_edge_connectivity Ugrid1d.get_connectivity_matrix Ugrid1d.copy @@ -300,6 +302,8 @@ UGRID2D Topology Ugrid2d.edge_node_connectivity Ugrid2d.face_edge_connectivity Ugrid2d.face_face_connectivity + Ugrid1d.directed_node_node_connectivity + Ugrid1d.directed_edge_edge_connectivity Ugrid2d.get_connectivity_matrix Ugrid2d.validate_edge_node_connectivity diff --git a/docs/changelog.rst b/docs/changelog.rst index 43aeb216..ba7b3f69 100644 --- a/docs/changelog.rst +++ b/docs/changelog.rst @@ -30,7 +30,8 @@ Added - This method is used in :meth:`xugrid.UgridDataArray.from_structured2d` and :meth:`xugrid.UgridDataset.from_structured2d` when the optional arguments ``x_bounds`` and ``y_bounds`` are provided. - +- Added :attr:`xugrid.Ugrid1d.directed_edge_edge_connectivity` and + :attr:`xugrid.Ugrid2d.directed_edge_edge_connectivity`. [0.12.2] 2025-01-31 ------------------- diff --git a/tests/test_connectivity.py b/tests/test_connectivity.py index 9be9bf63..ec772db3 100644 --- a/tests/test_connectivity.py +++ b/tests/test_connectivity.py @@ -492,6 +492,45 @@ def test_directed_node_node_connectivity(): assert np.array_equal(expected_edges, coo.data) +def test_directed_edge_edge_connectivity(): + r""" + Derive connectivity for this network: + + 5 + 4 --- 6 + / + / 3 + 0 1 2 / + 0 --- 1 --- 2 --- 3 + \ + \ 4 + \ + 5 + """ + edge_node_connectivity = np.array( + [ + [0, 1], # edge 0 + [1, 2], # edge 1 + [2, 3], # edge 2 + [3, 4], # edge 3 + [3, 5], # edge 4 + [4, 6], # edge 5 + ] + ) + node_edge_connectivity = connectivity.invert_dense_to_sparse(edge_node_connectivity) + csr = connectivity.directed_edge_edge_connectivity( + edge_node_connectivity, node_edge_connectivity + ) + assert isinstance(csr, sparse.csr_matrix) + + coo = csr.tocoo() + actual = np.column_stack([coo.row, coo.col]) + expected = np.array([[0, 1], [1, 2], [2, 3], [2, 4], [3, 5]]) + assert np.array_equal(actual, expected) + # Test through which node the connection is formed. + assert np.array_equal(coo.data, [1, 2, 3, 3, 4]) + + def test_face_face_connectivity(): edge_faces = np.array( [ diff --git a/tests/test_ugrid1d.py b/tests/test_ugrid1d.py index dd121589..2926343e 100644 --- a/tests/test_ugrid1d.py +++ b/tests/test_ugrid1d.py @@ -95,6 +95,8 @@ def test_ugrid1d_properties(): assert grid.bounds == (0.0, 0.0, 2.0, 2.0) assert isinstance(grid.node_edge_connectivity, sparse.csr_matrix) assert isinstance(grid.node_node_connectivity, sparse.csr_matrix) + assert isinstance(grid.directed_node_node_connectivity, sparse.csr_matrix) + assert isinstance(grid.directed_edge_edge_connectivity, sparse.csr_matrix) expected_coords = [ [[0.0, 0.0], [1.0, 1.0]], diff --git a/tests/test_ugrid2d.py b/tests/test_ugrid2d.py index 1b30f668..9ac4b2ea 100644 --- a/tests/test_ugrid2d.py +++ b/tests/test_ugrid2d.py @@ -151,6 +151,8 @@ def test_ugrid2d_properties(): assert grid.bounds == (0.0, 0.0, 2.0, 2.0) assert isinstance(grid.node_node_connectivity, sparse.csr_matrix) assert isinstance(grid.node_edge_connectivity, sparse.csr_matrix) + assert isinstance(grid.directed_node_node_connectivity, sparse.csr_matrix) + assert isinstance(grid.directed_edge_edge_connectivity, sparse.csr_matrix) edge_node_coords = grid.edge_node_coordinates face_node_coords = grid.face_node_coordinates assert edge_node_coords.shape == (10, 2, 2) diff --git a/xugrid/ugrid/connectivity.py b/xugrid/ugrid/connectivity.py index 5df555ed..0a17cbdf 100644 --- a/xugrid/ugrid/connectivity.py +++ b/xugrid/ugrid/connectivity.py @@ -526,6 +526,28 @@ def node_node_connectivity(edge_node_connectivity: IntArray) -> sparse.csr_matri return coo_matrix.tocsr() +def directed_edge_edge_connectivity( + edge_node_connectivity: IntArray, + node_edge_connectivity: sparse.csr_matrix, +) -> sparse.csr_matrix: + # - Get the second node of each edge. + # - Find which edges are connected to this second node. + # - This will include the edge we started with. + # - Hence, make room for one more. + # - Then filter self -> self away. + n_edge = len(edge_node_connectivity) + second_node = edge_node_connectivity[:, 1] + n_downstream = node_edge_connectivity.getnnz(axis=1)[second_node] + upstream_edges = np.repeat(np.arange(n_edge), n_downstream) + downstream_edges = node_edge_connectivity[second_node].indices + node_index = np.repeat(second_node, n_downstream) + valid = downstream_edges != upstream_edges + return sparse.csr_matrix( + (node_index[valid], (upstream_edges[valid], downstream_edges[valid])), + shape=(n_edge, n_edge), + ) + + def structured_connectivity(active: IntArray) -> AdjacencyMatrix: nrow, ncol = active.shape nodes = np.arange(nrow * ncol).reshape(nrow, ncol) diff --git a/xugrid/ugrid/ugridbase.py b/xugrid/ugrid/ugridbase.py index f94f164b..065bcc9d 100644 --- a/xugrid/ugrid/ugridbase.py +++ b/xugrid/ugrid/ugridbase.py @@ -720,6 +720,25 @@ def directed_node_node_connectivity(self) -> csr_matrix: """ return connectivity.directed_node_node_connectivity(self.edge_node_connectivity) + @property + def directed_edge_edge_connectivity(self) -> csr_matrix: + """ + Directed edge to edge connectivity. + + The connectivity is represented as an adjacency matrix in CSR format, + with the row and column indices as a (0-based) edge index. The data of + the matrix contains the node index of the common node through which + the connection is formed. + + Returns + ------- + connectivity: csr_matrix + """ + return connectivity.directed_edge_edge_connectivity( + self.edge_node_connectivity, + self.node_edge_connectivity, + ) + @staticmethod def _connectivity_weights(connectivity: csr_matrix, coordinates: FloatArray): xy = coordinates