From 517dffc784f33e9a9d1695211fedcffb290f3347 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Mon, 10 Nov 2025 16:44:07 +0000 Subject: [PATCH 1/6] Use adjacencylists and pure numpy for arrays. no more mixture of lists and numpy arrays --- src/networks_fenicsx/assembly.py | 31 ++++++----- src/networks_fenicsx/mesh.py | 88 +++++++++++++++++++++++--------- tests/test_edge_info.py | 57 +++++++++++++++++++++ 3 files changed, 137 insertions(+), 39 deletions(-) create mode 100644 tests/test_edge_info.py diff --git a/src/networks_fenicsx/assembly.py b/src/networks_fenicsx/assembly.py index 13b84eb..5909a81 100644 --- a/src/networks_fenicsx/assembly.py +++ b/src/networks_fenicsx/assembly.py @@ -42,22 +42,26 @@ def compute_integration_data( """ # Pack all bifurcation in and out fluxes per colored edge in graph - influx_color_to_bifurcations: dict[int, list[int]] = { - int(color): [] for color in range(network_mesh._num_edge_colors) + influx_color_to_bifurcations: dict[int, npt.NDArray[np.int32]] = { + int(color): np.empty(0, dtype=np.int32) for color in range(network_mesh.num_edge_colors) } - outflux_color_to_bifurcations: dict[int, list[int]] = { - int(color): [] for color in range(network_mesh._num_edge_colors) + outflux_color_to_bifurcations: dict[int, npt.NDArray[np.int32]] = { + int(color): np.empty(0, dtype=np.int32) for color in range(network_mesh.num_edge_colors) } - for bifurcation in network_mesh.bifurcation_values: - for color in network_mesh.in_edges(bifurcation): - influx_color_to_bifurcations[color].append(int(bifurcation)) - for color in network_mesh.out_edges(bifurcation): - outflux_color_to_bifurcations[color].append(int(bifurcation)) + for i, bifurcation in enumerate(network_mesh.bifurcation_values): + for color in network_mesh.in_edges(i): + influx_color_to_bifurcations[color] = np.append( + influx_color_to_bifurcations[color], int(bifurcation) + ) + for color in network_mesh.out_edges(i): + outflux_color_to_bifurcations[color] = np.append( + outflux_color_to_bifurcations[color], int(bifurcation) + ) # Accumulate integration data for all in-edges on the same submesh. in_flux_entities: dict[int, npt.NDArray[np.int32]] = {} out_flux_entities: dict[int, npt.NDArray[np.int32]] = {} - for color in range(network_mesh._num_edge_colors): + for color in range(network_mesh.num_edge_colors): sm = network_mesh.submeshes[color] smfm = network_mesh.submesh_facet_markers[color] sm.topology.create_connectivity(sm.topology.dim - 1, sm.topology.dim) @@ -187,7 +191,7 @@ def compute_forms( f: source term p_bc: neumann bc for pressure """ - num_qs = self._network_mesh._num_edge_colors + num_flux_spaces = self._network_mesh.num_edge_colors test_functions = [ufl.TestFunction(fs) for fs in self.function_spaces] trial_functions = [ufl.TrialFunction(fs) for fs in self.function_spaces] @@ -214,7 +218,6 @@ def compute_forms( network_mesh = self._network_mesh # Assemble edge contributions to a and L - num_qs = len(self._network_mesh.submeshes) P1_e = fem.functionspace(network_mesh.mesh, ("Lagrange", 1)) p_bc = fem.Function(P1_e) p_bc.interpolate(p_bc_ex.eval) @@ -230,8 +233,8 @@ def compute_forms( ds_edge = ufl.Measure("ds", domain=submesh, subdomain_data=facet_marker) a[i][i] += qs[i] * vs[i] * dx_edge - a[num_qs][i] += phi * ufl.dot(ufl.grad(qs[i]), tangent) * dx_edge - a[i][num_qs] = -p * ufl.dot(ufl.grad(vs[i]), tangent) * dx_edge + a[num_flux_spaces][i] += phi * ufl.dot(ufl.grad(qs[i]), tangent) * dx_edge + a[i][num_flux_spaces] = -p * ufl.dot(ufl.grad(vs[i]), tangent) * dx_edge # Add all boundary contributions L[i] = p_bc * vs[i] * ds_edge(network_mesh.in_marker) - p_bc * vs[i] * ds_edge( diff --git a/src/networks_fenicsx/mesh.py b/src/networks_fenicsx/mesh.py index 97bb983..d4ecdb3 100644 --- a/src/networks_fenicsx/mesh.py +++ b/src/networks_fenicsx/mesh.py @@ -60,8 +60,9 @@ class NetworkMesh: # Graph properties _geom_dim: int _num_edge_colors: int - _bifurcation_in_color: dict[int, list[int]] - _bifurcation_out_color: dict[int, list[int]] + _bifurcation_values: npt.NDArray[np.int32] + _bifurcation_in_color: _graph.AdjacencyList + _bifurcation_out_color: _graph.AdjacencyList # Mesh properties _msh: mesh.Mesh | None @@ -71,7 +72,6 @@ class NetworkMesh: _edge_meshes: list[mesh.Mesh] _edge_entity_maps: list[mesh.EntityMap] _tangent: fem.Function - _bifurcation_values: npt.NDArray[np.int32] _boundary_values: npt.NDArray[np.int32] _lm_mesh: mesh.Mesh | None _lm_map: mesh.EntityMap | None @@ -158,8 +158,10 @@ def _build_mesh( bifurcation_values = None boundary_values = None number_of_nodes = None - bifurcation_in_color: dict[int, list[int]] | None = None - bifurcation_out_color: dict[int, list[int]] | None = None + bifurcation_in_color: npt.NDArray[np.int32] | None = None + bifurcation_in_offsets: npt.NDArray[np.int32] | None = None + bifurcation_out_color: npt.NDArray[np.int32] | None = None + bifurcation_out_offsets: npt.NDArray[np.int32] | None = None if comm.rank == graph_rank: assert isinstance(graph, nx.DiGraph), f"Directional graph not present of {graph_rank}" geom_dim = len(graph.nodes[1]["pos"]) @@ -175,29 +177,41 @@ def _build_mesh( boundary_values = np.flatnonzero(nodes_with_degree == 1) max_connections = np.max(nodes_with_degree) - bifurcation_in_color = {} - bifurcation_out_color = {} + bifurcation_in_color = np.empty(0, dtype=np.int32) + bifurcation_in_offsets = np.zeros(1, dtype=np.int32) + bifurcation_out_color = np.empty(0, dtype=np.int32) + bifurcation_out_offsets = np.zeros(1, dtype=np.int32) for bifurcation in bifurcation_values: in_edges = graph.in_edges(bifurcation) - bifurcation_in_color[int(bifurcation)] = [] - for edge in in_edges: - bifurcation_in_color[int(bifurcation)].append(edge_coloring[edge]) + bifurcation_in_offsets = np.append( + bifurcation_in_offsets, len(bifurcation_in_color) + len(in_edges) + ) + bifurcation_in_color = np.append( + bifurcation_in_color, [edge_coloring[edge] for edge in in_edges] + ) out_edges = graph.out_edges(bifurcation) - bifurcation_out_color[int(bifurcation)] = [] - for edge in out_edges: - bifurcation_out_color[int(bifurcation)].append(edge_coloring[edge]) + bifurcation_out_offsets = np.append( + bifurcation_out_offsets, len(bifurcation_out_color) + len(out_edges) + ) + bifurcation_out_color = np.append( + bifurcation_out_color, [edge_coloring[edge] for edge in out_edges] + ) # Map boundary_values to inlet and outlet data from graph - boundary_in_nodes = [] - boundary_out_nodes = [] + boundary_in_nodes = np.empty(0, dtype=np.int32) + boundary_out_nodes = np.empty(0, dtype=np.int32) for boundary in boundary_values: in_edges = graph.in_edges(boundary) out_edges = graph.out_edges(boundary) assert len(in_edges) + len(out_edges) == 1, "Boundary node with multiple edges" if len(in_edges) == 1: - boundary_in_nodes.append(int(boundary)) + boundary_in_nodes = ( + np.append(boundary_in_nodes, np.int32(boundary)).flatten().astype(np.int32) + ) else: - boundary_out_nodes.append(int(boundary)) + boundary_out_nodes = ( + np.append(boundary_out_nodes, np.int32(boundary)).flatten().astype(np.int32) + ) comm.bcast(num_edge_colors, root=graph_rank) num_edge_colors, number_of_nodes, max_connections, geom_dim = comm.bcast( @@ -208,8 +222,19 @@ def _build_mesh( (bifurcation_values, boundary_values), root=graph_rank ) comm.barrier() - bifurcation_in_color, bifurcation_out_color = comm.bcast( - (bifurcation_in_color, bifurcation_out_color), root=graph_rank + ( + bifurcation_in_color, + bifurcation_in_offsets, + bifurcation_out_color, + bifurcation_out_offsets, + ) = comm.bcast( + ( + bifurcation_in_color, + bifurcation_in_offsets, + bifurcation_out_color, + bifurcation_out_offsets, + ), + root=graph_rank, ) comm.barrier() @@ -219,8 +244,12 @@ def _build_mesh( self._boundary_values = boundary_values # Create lookup of in and out colors for each bifurcation - self._bifurcation_in_color = bifurcation_in_color - self._bifurcation_out_color = bifurcation_out_color + self._bifurcation_in_color = _graph.adjacencylist( + bifurcation_in_color, bifurcation_in_offsets + ) + self._bifurcation_out_color = _graph.adjacencylist( + bifurcation_out_color, bifurcation_out_offsets + ) # Generate mesh tangents: npt.NDArray[np.float64] = np.empty((0, self._geom_dim), dtype=np.float64) @@ -480,11 +509,20 @@ def bifurcation_values(self) -> npt.NDArray[np.int32]: def boundary_values(self) -> npt.NDArray[np.int32]: return self._boundary_values - def in_edges(self, bifurcation: int) -> list[int]: - return self._bifurcation_in_color[bifurcation] + def in_edges(self, bifurcation_idx: int) -> npt.NDArray[np.int32 | np.int64]: + """Return the list of in-edge colors for a given bifurcation node. + Index is is the index of the bifurcation in {py:meth}`self.bifurcation_values`.""" + return self._bifurcation_in_color.links(np.int32(bifurcation_idx)) - def out_edges(self, bifurcation: int) -> list[int]: - return self._bifurcation_out_color[bifurcation] + def out_edges(self, bifurcation_idx: int) -> npt.NDArray[np.int32 | np.int64]: + """Return the list of out-edge colors for a given bifurcation node. + Index is is the index of the bifurcation in {py:meth}`self.bifurcation_values`.""" + return self._bifurcation_out_color.links(np.int32(bifurcation_idx)) + + @property + def num_edge_colors(self) -> int: + """Return the number of edge colors in the network mesh.""" + return self._num_edge_colors def _compute_tangent( diff --git a/tests/test_edge_info.py b/tests/test_edge_info.py new file mode 100644 index 0000000..7401803 --- /dev/null +++ b/tests/test_edge_info.py @@ -0,0 +1,57 @@ +import networkx as nx +import numpy as np +import pytest + +from networks_fenicsx import Config, NetworkMesh + + +@pytest.mark.parametrize("lcar", [0.1, 0.05]) +def test_edge_info(lcar: float): + # Create manual graph where we have five bifurcations. + # One inlet (0) -> (1) -> (7). + G = nx.DiGraph() + G.add_node(0, pos=np.zeros(3)) + G.add_node(1, pos=np.array([0.0, 0.0, 1.0])) + G.add_node(2, pos=np.array([0.2, 0.2, 2.0])) + G.add_node(3, pos=np.array([-0.2, 0.3, 2.0])) + G.add_node(4, pos=np.array([0.0, 0.1, 2.1])) + G.add_node(5, pos=np.array([0.1, -0.1, 3.0])) + G.add_node(6, pos=np.array([-0.3, 0.4, 4.0])) + G.add_node(7, pos=1.1 * G.nodes[1]["pos"]) + G.add_edge(0, 1) + G.add_edge(1, 7) + # Split into three bifurcations (7)->(2), (7)->(3), (7)->(4) + G.add_edge(7, 2) + # First branch goes right to gathering + G.add_edge(2, 5) + # Second branch goes through an extra point(3), before meeting at point (4) + G.add_edge(7, 3) + G.add_edge(3, 4) + G.add_edge(4, 5) + # Last branch goes directly to gathering + G.add_edge(7, 4) + G.add_edge(5, 6) + + config = Config() + config.lcar = lcar + network_mesh = NetworkMesh(G, config) + assert len(network_mesh.bifurcation_values) == 6 + # Bifurcation values are sorted in increasing order + np.testing.assert_allclose([1, 2, 3, 4, 5, 7], network_mesh.bifurcation_values) + assert len(network_mesh.in_edges(0)) == 1 + assert len(network_mesh.out_edges(0)) == 1 + + assert len(network_mesh.in_edges(1)) == 1 + assert len(network_mesh.out_edges(1)) == 1 + + assert len(network_mesh.in_edges(2)) == 1 + assert len(network_mesh.out_edges(2)) == 1 + + assert len(network_mesh.in_edges(3)) == 2 + assert len(network_mesh.out_edges(3)) == 1 + + assert len(network_mesh.in_edges(4)) == 2 + assert len(network_mesh.out_edges(4)) == 1 + + assert len(network_mesh.in_edges(5)) == 1 + assert len(network_mesh.out_edges(5)) == 3 From 431d0b489534b2154cbbc7bbf0983febfe88f815 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= Date: Mon, 10 Nov 2025 16:46:22 +0000 Subject: [PATCH 2/6] Add sanity check --- src/networks_fenicsx/mesh.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/networks_fenicsx/mesh.py b/src/networks_fenicsx/mesh.py index d4ecdb3..98a367d 100644 --- a/src/networks_fenicsx/mesh.py +++ b/src/networks_fenicsx/mesh.py @@ -512,11 +512,13 @@ def boundary_values(self) -> npt.NDArray[np.int32]: def in_edges(self, bifurcation_idx: int) -> npt.NDArray[np.int32 | np.int64]: """Return the list of in-edge colors for a given bifurcation node. Index is is the index of the bifurcation in {py:meth}`self.bifurcation_values`.""" + assert bifurcation_idx < len(self.bifurcation_values) return self._bifurcation_in_color.links(np.int32(bifurcation_idx)) def out_edges(self, bifurcation_idx: int) -> npt.NDArray[np.int32 | np.int64]: """Return the list of out-edge colors for a given bifurcation node. Index is is the index of the bifurcation in {py:meth}`self.bifurcation_values`.""" + assert bifurcation_idx < len(self.bifurcation_values) return self._bifurcation_out_color.links(np.int32(bifurcation_idx)) @property From 9c7765d9a661ead3bef0bc2e8d84c822869080b6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=B8rgen=20Schartum=20Dokken?= Date: Mon, 10 Nov 2025 20:17:39 +0100 Subject: [PATCH 3/6] Apply suggestions from code review MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Paul T. Kühner <56360279+schnellerhase@users.noreply.github.com> --- src/networks_fenicsx/assembly.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/networks_fenicsx/assembly.py b/src/networks_fenicsx/assembly.py index 5909a81..92ac5f4 100644 --- a/src/networks_fenicsx/assembly.py +++ b/src/networks_fenicsx/assembly.py @@ -51,11 +51,11 @@ def compute_integration_data( for i, bifurcation in enumerate(network_mesh.bifurcation_values): for color in network_mesh.in_edges(i): influx_color_to_bifurcations[color] = np.append( - influx_color_to_bifurcations[color], int(bifurcation) + influx_color_to_bifurcations[color], bifurcation ) for color in network_mesh.out_edges(i): outflux_color_to_bifurcations[color] = np.append( - outflux_color_to_bifurcations[color], int(bifurcation) + outflux_color_to_bifurcations[color], bifurcation ) # Accumulate integration data for all in-edges on the same submesh. in_flux_entities: dict[int, npt.NDArray[np.int32]] = {} From 83d41b6bee69fa9b6931ee1cd20c9cd9ec2c9f32 Mon Sep 17 00:00:00 2001 From: jorgensd Date: Tue, 11 Nov 2025 19:54:22 +0100 Subject: [PATCH 4/6] Remove duplicate definition --- src/networks_fenicsx/mesh.py | 6 ------ 1 file changed, 6 deletions(-) diff --git a/src/networks_fenicsx/mesh.py b/src/networks_fenicsx/mesh.py index 99edc3d..51e1449 100644 --- a/src/networks_fenicsx/mesh.py +++ b/src/networks_fenicsx/mesh.py @@ -442,12 +442,6 @@ def _build_network_submeshes(self): mesh.meshtags(edge_mesh, 0, marked_vertices, marked_values) ) - def in_edges(self, bifurcation: int) -> list[int]: - return self._bifurcation_in_color[bifurcation] - - def out_edges(self, bifurcation: int) -> list[int]: - return self._bifurcation_out_color[bifurcation] - @property def submesh_facet_markers(self) -> list[mesh.MeshTags]: if self._submesh_facet_markers is None: From 4da86aa82634946021f64efde2d3c3a94d72a43c Mon Sep 17 00:00:00 2001 From: jorgensd Date: Tue, 11 Nov 2025 19:56:44 +0100 Subject: [PATCH 5/6] Remove config --- tests/test_edge_info.py | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/tests/test_edge_info.py b/tests/test_edge_info.py index 7401803..7f8d684 100644 --- a/tests/test_edge_info.py +++ b/tests/test_edge_info.py @@ -2,11 +2,11 @@ import numpy as np import pytest -from networks_fenicsx import Config, NetworkMesh +from networks_fenicsx import NetworkMesh -@pytest.mark.parametrize("lcar", [0.1, 0.05]) -def test_edge_info(lcar: float): +@pytest.mark.parametrize("N", [10, 50]) +def test_edge_info(N: int): # Create manual graph where we have five bifurcations. # One inlet (0) -> (1) -> (7). G = nx.DiGraph() @@ -32,9 +32,7 @@ def test_edge_info(lcar: float): G.add_edge(7, 4) G.add_edge(5, 6) - config = Config() - config.lcar = lcar - network_mesh = NetworkMesh(G, config) + network_mesh = NetworkMesh(G, N=N) assert len(network_mesh.bifurcation_values) == 6 # Bifurcation values are sorted in increasing order np.testing.assert_allclose([1, 2, 3, 4, 5, 7], network_mesh.bifurcation_values) From 691b2b0bc64f90b58705669e14b824615ada3bfa Mon Sep 17 00:00:00 2001 From: schnellerhase <56360279+schnellerhase@users.noreply.github.com> Date: Mon, 24 Nov 2025 20:11:48 +0100 Subject: [PATCH 6/6] Deactivate output --- demos/demo_perf.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/demos/demo_perf.py b/demos/demo_perf.py index c697527..45a16d6 100644 --- a/demos/demo_perf.py +++ b/demos/demo_perf.py @@ -35,7 +35,7 @@ def eval(self, x): G = None network_mesh = NetworkMesh(G, N=1, color_strategy="smallest_last") del G - network_mesh.export_tangent() + # network_mesh.export() export_submeshes(network_mesh, outdir / f"n{n}") assembler = HydraulicNetworkAssembler(network_mesh, flux_degree=1, pressure_degree=0) # Compute forms