From c789ab2d255ce9426e647ddaaa98740cc7f57943 Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Fri, 7 Nov 2025 14:50:45 +0100
Subject: [PATCH 01/17] Tangent as UFL expression

---
 demos/demo_perf.py               |  1 -
 src/networks_fenicsx/assembly.py |  6 ++-
 src/networks_fenicsx/mesh.py     | 72 +-------------------------------
 3 files changed, 7 insertions(+), 72 deletions(-)

diff --git a/demos/demo_perf.py b/demos/demo_perf.py
index 2d11c0a..47401fd 100644
--- a/demos/demo_perf.py
+++ b/demos/demo_perf.py
@@ -49,7 +49,6 @@ def eval(self, x):
         G = None
     network_mesh = NetworkMesh(G, cfg)
     del G
-    network_mesh.export_tangent()
     export_submeshes(network_mesh, cfg.outdir / f"n{n}")
     assembler = HydraulicNetworkAssembler(cfg, network_mesh)
     # Compute forms
diff --git a/src/networks_fenicsx/assembly.py b/src/networks_fenicsx/assembly.py
index 13b84eb..c77cf5a 100644
--- a/src/networks_fenicsx/assembly.py
+++ b/src/networks_fenicsx/assembly.py
@@ -218,7 +218,11 @@ def compute_forms(
         P1_e = fem.functionspace(network_mesh.mesh, ("Lagrange", 1))
         p_bc = fem.Function(P1_e)
         p_bc.interpolate(p_bc_ex.eval)
-        tangent = self._network_mesh.tangent
+
+        tangent = ufl.Jacobian(self._network_mesh.mesh)
+        tangent = tangent[:, 0]
+        tangent /= ufl.sqrt(ufl.inner(tangent, tangent))
+
         for i, (submesh, entity_map, facet_marker) in enumerate(
             zip(
                 network_mesh.submeshes,
diff --git a/src/networks_fenicsx/mesh.py b/src/networks_fenicsx/mesh.py
index 97bb983..bbcae39 100644
--- a/src/networks_fenicsx/mesh.py
+++ b/src/networks_fenicsx/mesh.py
@@ -8,7 +8,7 @@
 https://github.com/IngeborgGjerde/fenics-networks by Ingeborg Gjerde.
 """
 
-from typing import Any, Callable, Iterable
+from typing import Callable, Iterable
 
 from mpi4py import MPI
 
@@ -18,9 +18,9 @@
 
 import basix.ufl
 import ufl
-from dolfinx import fem, mesh
 from dolfinx import graph as _graph
 from dolfinx import io as _io
+from dolfinx import mesh
 from networks_fenicsx import config
 
 __all__ = ["NetworkMesh"]
@@ -70,7 +70,6 @@ class NetworkMesh:
     _submesh_facet_markers: list[mesh.MeshTags]
     _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
@@ -223,7 +222,6 @@ def _build_mesh(
         self._bifurcation_out_color = bifurcation_out_color
 
         # Generate mesh
-        tangents: npt.NDArray[np.float64] = np.empty((0, self._geom_dim), dtype=np.float64)
         if comm.rank == graph_rank:
             assert isinstance(graph, nx.DiGraph), (
                 f"No directional graph present on rank {comm.rank}"
@@ -235,14 +233,10 @@ def _build_mesh(
             mesh_nodes = vertex_coords.copy()
             cells = []
             cell_markers = []
-            local_tangents = []
             if len(line_weights) == 0:
                 for segment in cells_array:
                     cells.append(np.array([segment[0], segment[1]], dtype=np.int64))
                     cell_markers.append(edge_coloring[(segment[0], segment[1])])
-                    start = vertex_coords[segment[0]]
-                    end = vertex_coords[segment[1]]
-                    local_tangents.append(_compute_tangent(start, end))
             else:
                 for segment in cells_array:
                     start_coord_pos = mesh_nodes.shape[0]
@@ -250,9 +244,6 @@ def _build_mesh(
                     end = vertex_coords[segment[1]]
                     internal_line_coords = start * (1 - line_weights) + end * line_weights
 
-                    tangent = _compute_tangent(start, end)
-                    local_tangents.append(np.tile(tangent, (internal_line_coords.shape[0] + 1, 1)))
-
                     mesh_nodes = np.vstack((mesh_nodes, internal_line_coords))
                     cells.append(np.array([segment[0], start_coord_pos], dtype=np.int64))
                     segment_connectivity = (
@@ -280,12 +271,10 @@ def _build_mesh(
                     )
             cells_ = np.vstack(cells).astype(np.int64)
             cell_markers_ = np.array(cell_markers, dtype=np.int32)
-            tangents = np.vstack(local_tangents).astype(np.float64)
         else:
             cells_ = np.empty((0, 2), dtype=np.int64)
             mesh_nodes = np.empty((0, self._geom_dim), dtype=np.float64)
             cell_markers_ = np.empty((0,), dtype=np.int32)
-        tangents = tangents[:, : self._geom_dim].copy()
         partitioner = mesh.create_cell_partitioner(mesh.GhostMode.shared_facet)
         graph_mesh = mesh.create_mesh(
             comm,
@@ -304,10 +293,6 @@ def _build_mesh(
             cell_markers_,
         )
 
-        # Distribute tangents
-        tangent_space = fem.functionspace(graph_mesh, ("DG", 0, (self._geom_dim,)))
-        self._tangent = fem.Function(tangent_space)
-
         # Which cells from the input mesh (rank 0) correspond to the local cells on this rank
         original_cell_indices = self.mesh.topology.original_cell_index[
             : self.mesh.topology.index_map(self.mesh.topology.dim).size_local
@@ -320,26 +305,6 @@ def _build_mesh(
             assert isinstance(recv_counts, Iterable)
             assert sum(recv_counts) == cells_.shape[0]
 
-        # Send to rank 0 which tangents that we need
-        recv_buffer: None | list[npt.NDArray[np.int64] | list[Any] | MPI.Datatype] = None
-        if comm.rank == graph_rank:
-            assert recv_counts is not None
-            recv_buffer = [np.empty(cells_.shape[0], dtype=np.int64), recv_counts, MPI.INT64_T]
-
-        comm.Gatherv(sendbuf=original_cell_indices, recvbuf=recv_buffer, root=graph_rank)
-        send_t_buffer = None
-        if comm.rank == graph_rank:
-            assert recv_buffer is not None
-            recv_data = recv_buffer[0]
-            assert isinstance(recv_data, np.ndarray)
-            send_t_buffer = [
-                tangents[recv_data].flatten(),
-                np.array(recv_counts) * self._geom_dim,
-            ]
-
-        comm.Scatterv(sendbuf=send_t_buffer, recvbuf=self.tangent.x.array, root=graph_rank)
-        self.tangent.x.scatter_forward()
-
         self._subdomains = mesh.meshtags_from_entities(
             self.mesh,
             self.mesh.topology.dim,
@@ -459,19 +424,6 @@ def entity_maps(self):
             )
         return self._edge_entity_maps
 
-    @property
-    def tangent(self):
-        """Return DG-0 field containing the tangent vector of the graph."""
-        return self._tangent
-
-    def export_tangent(self):
-        if self.cfg.export:
-            with _io.XDMFFile(self.comm, self.cfg.outdir / "mesh/tangent.xdmf", "w") as file:
-                file.write_mesh(self.mesh)
-                file.write_function(self.tangent)
-        else:
-            print("Export of tangent skipped as cfg.export is set to False.")
-
     @property
     def bifurcation_values(self) -> npt.NDArray[np.int32]:
         return self._bifurcation_values
@@ -485,23 +437,3 @@ def in_edges(self, bifurcation: int) -> list[int]:
 
     def out_edges(self, bifurcation: int) -> list[int]:
         return self._bifurcation_out_color[bifurcation]
-
-
-def _compute_tangent(
-    start: npt.NDArray[np.float64], end: npt.NDArray[np.float64]
-) -> npt.NDArray[np.float64]:
-    """Compute the tangent vector from start to end points.
-
-    Tangent is oriented according to positive y-axis.
-    If perpendicular to y-axis, align with x-axis.
-
-    Args:
-        start: Starting point coordinates.
-        end: Ending point coordinates.
-
-    Returns:
-        Normalized tangent vector from start to end.
-    """
-    tangent = end - start
-    tangent /= np.linalg.norm(tangent, axis=-1)
-    return tangent

From d8a8488b84cddc2a6f1e3609e3d6f99ef060f769 Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Sat, 8 Nov 2025 12:12:48 +0100
Subject: [PATCH 02/17] Introduce orientation

---
 .github/workflows/test_package.yml |  8 ++-
 demos/demo_perf.py                 |  1 +
 src/networks_fenicsx/assembly.py   |  7 +--
 src/networks_fenicsx/mesh.py       | 83 ++++++++++++++++++++++++------
 tests/test_orientation.py          | 64 +++++++++++++++++++++++
 5 files changed, 143 insertions(+), 20 deletions(-)
 create mode 100644 tests/test_orientation.py

diff --git a/.github/workflows/test_package.yml b/.github/workflows/test_package.yml
index ef83679..6d9577e 100644
--- a/.github/workflows/test_package.yml
+++ b/.github/workflows/test_package.yml
@@ -37,7 +37,13 @@ jobs:
       - name: Run tests
         run: python3 -m pytest -vxs
 
-      - name: Run tests in parallel
+      - name: Run tests in parallel (np=2)
+        run: mpiexec -n 3 python3 -m pytest -v  -m "not serial"
+
+      - name: Run tests in parallel (np=3)
+        run: mpiexec -n 3 python3 -m pytest -v  -m "not serial"
+
+      - name: Run tests in parallel (np=4)
         run: mpiexec -n 4 python3 -m pytest -v  -m "not serial"
 
       - name: Run demos (sequential)
diff --git a/demos/demo_perf.py b/demos/demo_perf.py
index 47401fd..8a2dcfa 100644
--- a/demos/demo_perf.py
+++ b/demos/demo_perf.py
@@ -49,6 +49,7 @@ def eval(self, x):
         G = None
     network_mesh = NetworkMesh(G, cfg)
     del G
+    network_mesh.export_orientation()
     export_submeshes(network_mesh, cfg.outdir / f"n{n}")
     assembler = HydraulicNetworkAssembler(cfg, network_mesh)
     # Compute forms
diff --git a/src/networks_fenicsx/assembly.py b/src/networks_fenicsx/assembly.py
index c77cf5a..e3f70ce 100644
--- a/src/networks_fenicsx/assembly.py
+++ b/src/networks_fenicsx/assembly.py
@@ -219,10 +219,11 @@ def compute_forms(
         p_bc = fem.Function(P1_e)
         p_bc.interpolate(p_bc_ex.eval)
 
-        tangent = ufl.Jacobian(self._network_mesh.mesh)
-        tangent = tangent[:, 0]
-        tangent /= ufl.sqrt(ufl.inner(tangent, tangent))
+        J = ufl.Jacobian(network_mesh.mesh)
+        t = J[:, 0]
+        t /= ufl.sqrt(ufl.inner(t, t))
 
+        tangent = self._network_mesh.orientation * t
         for i, (submesh, entity_map, facet_marker) in enumerate(
             zip(
                 network_mesh.submeshes,
diff --git a/src/networks_fenicsx/mesh.py b/src/networks_fenicsx/mesh.py
index bbcae39..55524c7 100644
--- a/src/networks_fenicsx/mesh.py
+++ b/src/networks_fenicsx/mesh.py
@@ -18,9 +18,9 @@
 
 import basix.ufl
 import ufl
+from dolfinx import fem, mesh
 from dolfinx import graph as _graph
 from dolfinx import io as _io
-from dolfinx import mesh
 from networks_fenicsx import config
 
 __all__ = ["NetworkMesh"]
@@ -70,6 +70,7 @@ class NetworkMesh:
     _submesh_facet_markers: list[mesh.MeshTags]
     _edge_meshes: list[mesh.Mesh]
     _edge_entity_maps: list[mesh.EntityMap]
+    _orientation: fem.Function
     _bifurcation_values: npt.NDArray[np.int32]
     _boundary_values: npt.NDArray[np.int32]
     _lm_mesh: mesh.Mesh | None
@@ -233,10 +234,15 @@ def _build_mesh(
             mesh_nodes = vertex_coords.copy()
             cells = []
             cell_markers = []
+            orientations = np.array([], dtype=np.float64)
             if len(line_weights) == 0:
                 for segment in cells_array:
                     cells.append(np.array([segment[0], segment[1]], dtype=np.int64))
                     cell_markers.append(edge_coloring[(segment[0], segment[1])])
+                    start = vertex_coords[segment[0]]
+                    end = vertex_coords[segment[1]]
+                    in_order = segment[0] < segment[1]
+                    orientations = np.append(orientations, [1 if in_order else -1])
             else:
                 for segment in cells_array:
                     start_coord_pos = mesh_nodes.shape[0]
@@ -244,6 +250,9 @@ def _build_mesh(
                     end = vertex_coords[segment[1]]
                     internal_line_coords = start * (1 - line_weights) + end * line_weights
 
+                    in_order = segment[0] < segment[1]
+                    orientations = np.append(orientations, [1 if in_order else -1])
+
                     mesh_nodes = np.vstack((mesh_nodes, internal_line_coords))
                     cells.append(np.array([segment[0], start_coord_pos], dtype=np.int64))
                     segment_connectivity = (
@@ -275,6 +284,8 @@ def _build_mesh(
             cells_ = np.empty((0, 2), dtype=np.int64)
             mesh_nodes = np.empty((0, self._geom_dim), dtype=np.float64)
             cell_markers_ = np.empty((0,), dtype=np.int32)
+            orientations = np.empty(0, dtype=np.float64)
+
         partitioner = mesh.create_cell_partitioner(mesh.GhostMode.shared_facet)
         graph_mesh = mesh.create_mesh(
             comm,
@@ -286,24 +297,51 @@ def _build_mesh(
         )
         self._msh = graph_mesh
 
-        local_entities, local_values = _io.distribute_entity_data(
-            self.mesh,
-            self.mesh.topology.dim,
-            cells_,
-            cell_markers_,
+        tdim = self.mesh.topology.dim
+
+        # Distribute orientations as meshtag
+        local_cells, local_orientations = _io.distribute_entity_data(
+            self.mesh, tdim, cells_, orientations
         )
 
-        # Which cells from the input mesh (rank 0) correspond to the local cells on this rank
-        original_cell_indices = self.mesh.topology.original_cell_index[
-            : self.mesh.topology.index_map(self.mesh.topology.dim).size_local
-        ]
+        meshtag_orientation = mesh.meshtags_from_entities(
+            self.mesh, tdim, _graph.adjacencylist(local_cells), local_orientations
+        )
 
-        # First compute how much data we are getting from each rank
-        send_count = original_cell_indices.size
-        recv_counts = comm.gather(send_count, root=graph_rank)
-        if comm.rank == graph_rank:
-            assert isinstance(recv_counts, Iterable)
-            assert sum(recv_counts) == cells_.shape[0]
+        # Assemble orientations function
+        orientation_space = fem.functionspace(graph_mesh, ("DG", 0))
+        self._orientation = fem.Function(orientation_space)
+        self._orientation.x.array[orientation_space.dofmap.list.flatten()] = (
+            meshtag_orientation.values
+        )
+
+        # Correct orientations for possible reorder
+        e_to_f = self.mesh.topology.connectivity(tdim, tdim - 1).array.reshape(-1, 2)
+
+        e_idx = np.arange(self.mesh.topology.index_map(tdim).size_local)
+        e = e_to_f[e_idx]
+        global_input = self.mesh.geometry.input_global_indices
+        in_order = global_input[e[:, 0]] < global_input[e[:, 1]]
+
+        # Four cases that might arise (per edge), with naming i := input_in_order, n := in_order:
+        # 1) i & n:
+        #   orientation = 1 (input=1, no action needed)
+        # 2) i & !n:
+        #   orientation = -1 (input=1, need sign flip)
+        # 3) !i & n:
+        #   orientation = -1 (input=-1, no action needed)
+        # 4) !i & !n
+        #   orientation = 1 (input=-1, need sign flip)
+        #
+        # so we need to sign flip when
+        #   (i & !n) | (!i & !n) = !n
+        self.orientation.x.array[: e_idx.size][~in_order] *= -1
+
+        self.orientation.x.scatter_forward()
+
+        local_entities, local_values = _io.distribute_entity_data(
+            self.mesh, tdim, cells_, cell_markers_
+        )
 
         self._subdomains = mesh.meshtags_from_entities(
             self.mesh,
@@ -424,6 +462,19 @@ def entity_maps(self):
             )
         return self._edge_entity_maps
 
+    @property
+    def orientation(self):
+        """Return DG-0 field containing the tangent vector of the graph."""
+        return self._orientation
+
+    def export_orientation(self):
+        if self.cfg.export:
+            with _io.XDMFFile(self.comm, self.cfg.outdir / "mesh/orientation.xdmf", "w") as file:
+                file.write_mesh(self.mesh)
+                file.write_function(self.orientation)
+        else:
+            print("Export of tangent skipped as cfg.export is set to False.")
+
     @property
     def bifurcation_values(self) -> npt.NDArray[np.int32]:
         return self._bifurcation_values
diff --git a/tests/test_orientation.py b/tests/test_orientation.py
new file mode 100644
index 0000000..7acba38
--- /dev/null
+++ b/tests/test_orientation.py
@@ -0,0 +1,64 @@
+import networkx as nx
+import numpy as np
+import pytest
+
+import dolfinx
+import ufl
+from networks_fenicsx.config import Config
+from networks_fenicsx.mesh import NetworkMesh
+
+
+def linear_graph(n: int, dim: int = 2, ordered=lambda _: True) -> nx.DiGraph:
+    G = nx.DiGraph()
+    G.add_nodes_from(range(n))
+
+    for i in range(n - 1):
+        if ordered(i):
+            G.add_edge(i, i + 1)
+        else:
+            G.add_edge(i + 1, i)
+
+    for i in range(n):
+        pos = np.zeros(dim)
+        pos[0] = i / (n - 1)
+        G.nodes[i]["pos"] = pos
+
+    return G
+
+
+@pytest.mark.parametrize("n", [30])
+@pytest.mark.parametrize("order", ["in", "reverse", "alternating"])
+def test_orientation(n: int, order: str) -> None:
+    if order == "in":
+        ordered = lambda _: True
+    elif order == "reverse":
+        ordered = lambda _: False
+    elif order == "alternating":
+        ordered = lambda k: k % 2
+    else:
+        raise RuntimeError()
+
+    G = linear_graph(n, ordered=ordered)
+    # import matplotlib.pyplot as plt
+    # nx.draw(G)
+    # plt.show()
+
+    cfg = Config()
+    cfg.export = False
+
+    network_mesh = NetworkMesh(G, cfg)
+
+    J = ufl.Jacobian(network_mesh.mesh)
+    t = J[:, 0]
+    t /= ufl.sqrt(ufl.inner(t, t))
+    f = dolfinx.fem.form(ufl.inner(ufl.as_vector((1, 0)), t) * network_mesh.orientation * ufl.dx)
+
+    val = network_mesh.comm.allreduce(dolfinx.fem.assemble_scalar(f))
+
+    if order == "in":
+        assert np.isclose(val, 1.0)
+    elif order == "reverse":
+        assert np.isclose(val, -1.0)
+    else:
+        edge_count = n - 1
+        assert np.isclose(val, edge_count % 2 * -1 / edge_count)

From e047ad3b202c45160bdd2bfae9a413d3ef27a0cd Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Paul=20T=2E=20K=C3=BChner?=
 <56360279+schnellerhase@users.noreply.github.com>
Date: Mon, 10 Nov 2025 17:55:08 +0100
Subject: [PATCH 03/17] Update src/networks_fenicsx/mesh.py
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Co-authored-by: Jørgen Schartum Dokken <dokken92@gmail.com>
---
 src/networks_fenicsx/mesh.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/networks_fenicsx/mesh.py b/src/networks_fenicsx/mesh.py
index 55524c7..574bafb 100644
--- a/src/networks_fenicsx/mesh.py
+++ b/src/networks_fenicsx/mesh.py
@@ -318,7 +318,7 @@ def _build_mesh(
         # Correct orientations for possible reorder
         e_to_f = self.mesh.topology.connectivity(tdim, tdim - 1).array.reshape(-1, 2)
 
-        e_idx = np.arange(self.mesh.topology.index_map(tdim).size_local)
+        e_idx = np.arange(self.mesh.topology.index_map(tdim).size_local, dtype=np.int32)
         e = e_to_f[e_idx]
         global_input = self.mesh.geometry.input_global_indices
         in_order = global_input[e[:, 0]] < global_input[e[:, 1]]

From ae1353a655f5365fc59be2bb17dd586a915e472e Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Mon, 10 Nov 2025 18:10:41 +0100
Subject: [PATCH 04/17] Reorder and meshtag indices

---
 src/networks_fenicsx/mesh.py | 29 ++++++++++++++---------------
 1 file changed, 14 insertions(+), 15 deletions(-)

diff --git a/src/networks_fenicsx/mesh.py b/src/networks_fenicsx/mesh.py
index 574bafb..9f0b6a3 100644
--- a/src/networks_fenicsx/mesh.py
+++ b/src/networks_fenicsx/mesh.py
@@ -299,7 +299,19 @@ def _build_mesh(
 
         tdim = self.mesh.topology.dim
 
-        # Distribute orientations as meshtag
+        # Compute subdomain tags.
+        local_entities, local_values = _io.distribute_entity_data(
+            self.mesh, tdim, cells_, cell_markers_
+        )
+
+        self._subdomains = mesh.meshtags_from_entities(
+            self.mesh,
+            self.mesh.topology.dim,
+            _graph.adjacencylist(local_entities),
+            local_values,
+        )
+
+        # Distribute orientations as meshtag.
         local_cells, local_orientations = _io.distribute_entity_data(
             self.mesh, tdim, cells_, orientations
         )
@@ -311,9 +323,7 @@ def _build_mesh(
         # Assemble orientations function
         orientation_space = fem.functionspace(graph_mesh, ("DG", 0))
         self._orientation = fem.Function(orientation_space)
-        self._orientation.x.array[orientation_space.dofmap.list.flatten()] = (
-            meshtag_orientation.values
-        )
+        self._orientation.x.array[meshtag_orientation.indices] = meshtag_orientation.values
 
         # Correct orientations for possible reorder
         e_to_f = self.mesh.topology.connectivity(tdim, tdim - 1).array.reshape(-1, 2)
@@ -339,17 +349,6 @@ def _build_mesh(
 
         self.orientation.x.scatter_forward()
 
-        local_entities, local_values = _io.distribute_entity_data(
-            self.mesh, tdim, cells_, cell_markers_
-        )
-
-        self._subdomains = mesh.meshtags_from_entities(
-            self.mesh,
-            self.mesh.topology.dim,
-            _graph.adjacencylist(local_entities),
-            local_values,
-        )
-
         self._in_marker = 3 * number_of_nodes
         self._out_marker = 5 * number_of_nodes
         if comm.rank == graph_rank:

From a7a65b689e6b6ee8b0c8b4e0a6d5f242b6953055 Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Mon, 10 Nov 2025 18:19:04 +0100
Subject: [PATCH 05/17] Fix silent assumption on topology to geometry mapping

---
 src/networks_fenicsx/mesh.py | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/src/networks_fenicsx/mesh.py b/src/networks_fenicsx/mesh.py
index 9f0b6a3..3f11e5a 100644
--- a/src/networks_fenicsx/mesh.py
+++ b/src/networks_fenicsx/mesh.py
@@ -326,12 +326,11 @@ def _build_mesh(
         self._orientation.x.array[meshtag_orientation.indices] = meshtag_orientation.values
 
         # Correct orientations for possible reorder
-        e_to_f = self.mesh.topology.connectivity(tdim, tdim - 1).array.reshape(-1, 2)
-
         e_idx = np.arange(self.mesh.topology.index_map(tdim).size_local, dtype=np.int32)
-        e = e_to_f[e_idx]
+        e_geo = mesh.entities_to_geometry(self.mesh, tdim, e_idx)
+
         global_input = self.mesh.geometry.input_global_indices
-        in_order = global_input[e[:, 0]] < global_input[e[:, 1]]
+        in_order = global_input[e_geo[:, 0]] < global_input[e_geo[:, 1]]
 
         # Four cases that might arise (per edge), with naming i := input_in_order, n := in_order:
         # 1) i & n:

From b22b7e8bfff108fe3a94f18a5c3c5f453b5a42bc Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Mon, 10 Nov 2025 18:49:11 +0100
Subject: [PATCH 06/17] Extend to non unit lcar

---
 src/networks_fenicsx/mesh.py | 17 ++++++++++++++++-
 tests/test_orientation.py    |  4 +++-
 2 files changed, 19 insertions(+), 2 deletions(-)

diff --git a/src/networks_fenicsx/mesh.py b/src/networks_fenicsx/mesh.py
index 3f11e5a..a89f4dc 100644
--- a/src/networks_fenicsx/mesh.py
+++ b/src/networks_fenicsx/mesh.py
@@ -251,7 +251,6 @@ def _build_mesh(
                     internal_line_coords = start * (1 - line_weights) + end * line_weights
 
                     in_order = segment[0] < segment[1]
-                    orientations = np.append(orientations, [1 if in_order else -1])
 
                     mesh_nodes = np.vstack((mesh_nodes, internal_line_coords))
                     cells.append(np.array([segment[0], start_coord_pos], dtype=np.int64))
@@ -278,8 +277,24 @@ def _build_mesh(
                             dtype=np.int32,
                         )
                     )
+                    orientations = np.append(
+                        orientations, [1 if segment[0] < start_coord_pos else -1]
+                    )
+                    orientations = np.append(
+                        orientations, np.full(internal_line_coords.shape[0] - 1, 1)
+                    )
+                    orientations = np.append(
+                        orientations,
+                        [
+                            1
+                            if start_coord_pos + internal_line_coords.shape[0] - 1 < segment[1]
+                            else -1
+                        ],
+                    )
             cells_ = np.vstack(cells).astype(np.int64)
             cell_markers_ = np.array(cell_markers, dtype=np.int32)
+
+            assert cell_markers_.shape == orientations.shape
         else:
             cells_ = np.empty((0, 2), dtype=np.int64)
             mesh_nodes = np.empty((0, self._geom_dim), dtype=np.float64)
diff --git a/tests/test_orientation.py b/tests/test_orientation.py
index 7acba38..629b9c6 100644
--- a/tests/test_orientation.py
+++ b/tests/test_orientation.py
@@ -28,7 +28,8 @@ def linear_graph(n: int, dim: int = 2, ordered=lambda _: True) -> nx.DiGraph:
 
 @pytest.mark.parametrize("n", [30])
 @pytest.mark.parametrize("order", ["in", "reverse", "alternating"])
-def test_orientation(n: int, order: str) -> None:
+@pytest.mark.parametrize("lcar", [1, 1 / 4, 1 / 8])
+def test_orientation(n: int, order: str, lcar: int) -> None:
     if order == "in":
         ordered = lambda _: True
     elif order == "reverse":
@@ -45,6 +46,7 @@ def test_orientation(n: int, order: str) -> None:
 
     cfg = Config()
     cfg.export = False
+    cfg.lcar = lcar
 
     network_mesh = NetworkMesh(G, cfg)
 

From e44040dc047a5717615b19fdcc0ae17fa420959b Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Mon, 10 Nov 2025 19:04:46 +0100
Subject: [PATCH 07/17] Compute from cells

---
 src/networks_fenicsx/mesh.py | 20 ++++----------------
 1 file changed, 4 insertions(+), 16 deletions(-)

diff --git a/src/networks_fenicsx/mesh.py b/src/networks_fenicsx/mesh.py
index a89f4dc..91e97fc 100644
--- a/src/networks_fenicsx/mesh.py
+++ b/src/networks_fenicsx/mesh.py
@@ -234,7 +234,6 @@ def _build_mesh(
             mesh_nodes = vertex_coords.copy()
             cells = []
             cell_markers = []
-            orientations = np.array([], dtype=np.float64)
             if len(line_weights) == 0:
                 for segment in cells_array:
                     cells.append(np.array([segment[0], segment[1]], dtype=np.int64))
@@ -242,7 +241,6 @@ def _build_mesh(
                     start = vertex_coords[segment[0]]
                     end = vertex_coords[segment[1]]
                     in_order = segment[0] < segment[1]
-                    orientations = np.append(orientations, [1 if in_order else -1])
             else:
                 for segment in cells_array:
                     start_coord_pos = mesh_nodes.shape[0]
@@ -277,23 +275,13 @@ def _build_mesh(
                             dtype=np.int32,
                         )
                     )
-                    orientations = np.append(
-                        orientations, [1 if segment[0] < start_coord_pos else -1]
-                    )
-                    orientations = np.append(
-                        orientations, np.full(internal_line_coords.shape[0] - 1, 1)
-                    )
-                    orientations = np.append(
-                        orientations,
-                        [
-                            1
-                            if start_coord_pos + internal_line_coords.shape[0] - 1 < segment[1]
-                            else -1
-                        ],
-                    )
+
             cells_ = np.vstack(cells).astype(np.int64)
             cell_markers_ = np.array(cell_markers, dtype=np.int32)
 
+            orientations= np.full_like(cell_markers_, 1.0, dtype=np.float64)
+            orientations[cells_[:,0] > cells_[:,1]] = -1.0
+
             assert cell_markers_.shape == orientations.shape
         else:
             cells_ = np.empty((0, 2), dtype=np.int64)

From a9998d1a3cbcffd5c26cf52f19b045c1409b68b6 Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Mon, 10 Nov 2025 19:25:59 +0100
Subject: [PATCH 08/17] ruff

---
 src/networks_fenicsx/mesh.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/networks_fenicsx/mesh.py b/src/networks_fenicsx/mesh.py
index 91e97fc..110459b 100644
--- a/src/networks_fenicsx/mesh.py
+++ b/src/networks_fenicsx/mesh.py
@@ -279,8 +279,8 @@ def _build_mesh(
             cells_ = np.vstack(cells).astype(np.int64)
             cell_markers_ = np.array(cell_markers, dtype=np.int32)
 
-            orientations= np.full_like(cell_markers_, 1.0, dtype=np.float64)
-            orientations[cells_[:,0] > cells_[:,1]] = -1.0
+            orientations = np.full_like(cell_markers_, 1.0, dtype=np.float64)
+            orientations[cells_[:, 0] > cells_[:, 1]] = -1.0
 
             assert cell_markers_.shape == orientations.shape
         else:

From 70f81028b891dd50b2e3e6ca183674b04e85d544 Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Tue, 11 Nov 2025 17:54:34 +0100
Subject: [PATCH 09/17] Addapt up to edge coloring

---
 demos/demo_perf.py           |  6 ++++--
 src/networks_fenicsx/mesh.py |  8 --------
 tests/test_orientation.py    | 11 +++--------
 3 files changed, 7 insertions(+), 18 deletions(-)

diff --git a/demos/demo_perf.py b/demos/demo_perf.py
index 0e7863e..0ac70dc 100644
--- a/demos/demo_perf.py
+++ b/demos/demo_perf.py
@@ -3,7 +3,7 @@
 
 from mpi4py import MPI
 
-from dolfinx.io import VTXWriter
+from dolfinx.io import VTXWriter, XDMFFile
 from networks_fenicsx import (
     HydraulicNetworkAssembler,
     NetworkMesh,
@@ -35,7 +35,9 @@ def eval(self, x):
         G = None
     network_mesh = NetworkMesh(G, N=1, color_strategy="smallest_last")
     del G
-    network_mesh.export_orientation()
+    with XDMFFile(network_mesh.comm, outdir / "mesh" / "orientation.xdmf", "w") as file:
+        file.write_mesh(network_mesh.mesh)
+        file.write_function(network_mesh.orientation)
     export_submeshes(network_mesh, outdir / f"n{n}")
     assembler = HydraulicNetworkAssembler(network_mesh, flux_degree=1, pressure_degree=0)
     # Compute forms
diff --git a/src/networks_fenicsx/mesh.py b/src/networks_fenicsx/mesh.py
index 7e157db..4389b62 100644
--- a/src/networks_fenicsx/mesh.py
+++ b/src/networks_fenicsx/mesh.py
@@ -470,14 +470,6 @@ def orientation(self):
         """Return DG-0 field containing the tangent vector of the graph."""
         return self._orientation
 
-    def export_orientation(self):
-        if self.cfg.export:
-            with _io.XDMFFile(self.comm, self.cfg.outdir / "mesh/orientation.xdmf", "w") as file:
-                file.write_mesh(self.mesh)
-                file.write_function(self.orientation)
-        else:
-            print("Export of tangent skipped as cfg.export is set to False.")
-
     @property
     def bifurcation_values(self) -> npt.NDArray[np.int32]:
         return self._bifurcation_values
diff --git a/tests/test_orientation.py b/tests/test_orientation.py
index 629b9c6..a43ce73 100644
--- a/tests/test_orientation.py
+++ b/tests/test_orientation.py
@@ -4,7 +4,6 @@
 
 import dolfinx
 import ufl
-from networks_fenicsx.config import Config
 from networks_fenicsx.mesh import NetworkMesh
 
 
@@ -28,8 +27,8 @@ def linear_graph(n: int, dim: int = 2, ordered=lambda _: True) -> nx.DiGraph:
 
 @pytest.mark.parametrize("n", [30])
 @pytest.mark.parametrize("order", ["in", "reverse", "alternating"])
-@pytest.mark.parametrize("lcar", [1, 1 / 4, 1 / 8])
-def test_orientation(n: int, order: str, lcar: int) -> None:
+@pytest.mark.parametrize("N", [1, 4, 8])
+def test_orientation(n: int, order: str, N: int) -> None:
     if order == "in":
         ordered = lambda _: True
     elif order == "reverse":
@@ -44,11 +43,7 @@ def test_orientation(n: int, order: str, lcar: int) -> None:
     # nx.draw(G)
     # plt.show()
 
-    cfg = Config()
-    cfg.export = False
-    cfg.lcar = lcar
-
-    network_mesh = NetworkMesh(G, cfg)
+    network_mesh = NetworkMesh(G, N=N)
 
     J = ufl.Jacobian(network_mesh.mesh)
     t = J[:, 0]

From 53aaf7594a73899533b7d909ec2cc90ab7730222 Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Thu, 13 Nov 2025 16:13:36 +0100
Subject: [PATCH 10/17] Fix rebase

---
 src/networks_fenicsx/assembly.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/src/networks_fenicsx/assembly.py b/src/networks_fenicsx/assembly.py
index 44b0a8a..630102f 100644
--- a/src/networks_fenicsx/assembly.py
+++ b/src/networks_fenicsx/assembly.py
@@ -215,7 +215,7 @@ def compute_forms(
         # 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)
         if isinstance(p_bc_ex, ufl.core.expr.Expr):
             try:
                 expr = fem.Expression(p_bc_ex, P1_e.element.interpolation_points())  # type: ignore[operator]
@@ -230,7 +230,6 @@ def compute_forms(
         t /= ufl.sqrt(ufl.inner(t, t))
 
         tangent = self._network_mesh.orientation * t
-        tangent = self._network_mesh.tangent
         for i, (submesh, entity_map, facet_marker) in enumerate(
             zip(
                 network_mesh.submeshes,

From 0bc96d4354620a0bf416eef30a8f270c53ed8fe3 Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Tue, 17 Feb 2026 19:33:08 +0100
Subject: [PATCH 11/17] Fixup merge

---
 src/networks_fenicsx/assembly.py | 2 +-
 src/networks_fenicsx/mesh.py     | 3 +--
 2 files changed, 2 insertions(+), 3 deletions(-)

diff --git a/src/networks_fenicsx/assembly.py b/src/networks_fenicsx/assembly.py
index 4cb4666..cd0c843 100644
--- a/src/networks_fenicsx/assembly.py
+++ b/src/networks_fenicsx/assembly.py
@@ -229,7 +229,7 @@ def compute_forms(
             p_bc.interpolate(expr)
         else:
             p_bc.interpolate(p_bc_ex)
-        
+
         dx_global = ufl.Measure("dx", domain=network_mesh.mesh)
 
         J = ufl.Jacobian(network_mesh.mesh)
diff --git a/src/networks_fenicsx/mesh.py b/src/networks_fenicsx/mesh.py
index fc966c3..0167576 100644
--- a/src/networks_fenicsx/mesh.py
+++ b/src/networks_fenicsx/mesh.py
@@ -9,7 +9,7 @@
 """
 
 import inspect
-from typing import Any, Callable, Iterable
+from typing import Callable, Iterable
 
 from mpi4py import MPI
 
@@ -296,7 +296,6 @@ def _build_mesh(
             cell_markers_ = np.empty((0,), dtype=np.int32)
             orientations = np.empty(0, dtype=np.float64)
 
-        tangents = tangents[:, : self._geom_dim].copy()
         sig = inspect.signature(mesh.create_cell_partitioner)
         part_kwargs = {}
         if "max_facet_to_cell_links" in list(sig.parameters.keys()):

From 5362e8a1e1d1e8fa32d0d2a34ff0127cd939897c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= <dokken@simula.no>
Date: Wed, 18 Feb 2026 16:14:31 +0100
Subject: [PATCH 12/17] Fixes for dtype

---
 src/networks_fenicsx/mesh.py | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/src/networks_fenicsx/mesh.py b/src/networks_fenicsx/mesh.py
index 8cd1317..8db2971 100644
--- a/src/networks_fenicsx/mesh.py
+++ b/src/networks_fenicsx/mesh.py
@@ -177,7 +177,6 @@ def _build_mesh(
             assert isinstance(graph, nx.DiGraph), f"Directional graph not present of {graph_rank}"
             geom_dim = len(graph.nodes[1]["pos"])
             edge_coloring = color_graph(graph, color_strategy)
-
             num_edge_colors = len(set(edge_coloring.values()))
             cells_array = np.asarray([[u, v] for u, v in graph.edges()])
             number_of_nodes = graph.number_of_nodes()
@@ -198,14 +197,16 @@ def _build_mesh(
                     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]
+                    bifurcation_in_color,
+                    np.asarray([edge_coloring[edge] for edge in in_edges], dtype=np.int32),
                 )
                 out_edges = graph.out_edges(bifurcation)
                 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]
+                    bifurcation_out_color,
+                    np.asarray([edge_coloring[edge] for edge in out_edges], dtype=np.int32),
                 )
 
             # Map boundary_values to inlet and outlet data from graph

From ced5bb0e1186b218b6801b1b077776ef73423843 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=B8rgen=20S=2E=20Dokken?= <dokken@simula.no>
Date: Wed, 18 Feb 2026 16:27:04 +0100
Subject: [PATCH 13/17] Remove duplicate value

---
 src/networks_fenicsx/mesh.py | 1 -
 1 file changed, 1 deletion(-)

diff --git a/src/networks_fenicsx/mesh.py b/src/networks_fenicsx/mesh.py
index 8db2971..611b5d4 100644
--- a/src/networks_fenicsx/mesh.py
+++ b/src/networks_fenicsx/mesh.py
@@ -76,7 +76,6 @@ class NetworkMesh:
     _edge_meshes: list[mesh.Mesh]
     _edge_entity_maps: list[mesh.EntityMap]
     _orientation: fem.Function
-    _bifurcation_values: npt.NDArray[np.int32]
 
     _boundary_values: npt.NDArray[np.int32]
     _lm_mesh: mesh.Mesh | None

From 2aae09c578f36af57fe062f5bff81bff02f6414f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Paul=20T=2E=20K=C3=BChner?=
 <56360279+schnellerhase@users.noreply.github.com>
Date: Wed, 18 Feb 2026 16:36:20 +0100
Subject: [PATCH 14/17] Update tests/test_orientation.py
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Co-authored-by: Jørgen Schartum Dokken <dokken92@gmail.com>
---
 tests/test_orientation.py | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/tests/test_orientation.py b/tests/test_orientation.py
index a43ce73..599926a 100644
--- a/tests/test_orientation.py
+++ b/tests/test_orientation.py
@@ -39,9 +39,6 @@ def test_orientation(n: int, order: str, N: int) -> None:
         raise RuntimeError()
 
     G = linear_graph(n, ordered=ordered)
-    # import matplotlib.pyplot as plt
-    # nx.draw(G)
-    # plt.show()
 
     network_mesh = NetworkMesh(G, N=N)
 

From 4f743d97c5c071d1f747008fe192989e15395328 Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Wed, 18 Feb 2026 17:30:07 +0100
Subject: [PATCH 15/17] Pull out order functionality

---
 src/networks_fenicsx/mesh.py | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/src/networks_fenicsx/mesh.py b/src/networks_fenicsx/mesh.py
index 611b5d4..438b42c 100644
--- a/src/networks_fenicsx/mesh.py
+++ b/src/networks_fenicsx/mesh.py
@@ -262,6 +262,10 @@ def _build_mesh(
             bifurcation_out_color, bifurcation_out_offsets
         )
 
+        def _original_order(a, b):
+            """Order on edge entities to detect order reversals during mesh generation."""
+            return a < b
+
         # Generate mesh
         if comm.rank == graph_rank:
             assert isinstance(graph, nx.DiGraph), (
@@ -278,7 +282,7 @@ def _build_mesh(
                     cell_markers.append(edge_coloring[(segment[0], segment[1])])
                     start = vertex_coords[segment[0]]
                     end = vertex_coords[segment[1]]
-                    in_order = segment[0] < segment[1]
+                    in_order = _original_order(segment[0], segment[1])
             else:
                 for segment in cells_array:
                     start_coord_pos = mesh_nodes.shape[0]
@@ -286,7 +290,7 @@ def _build_mesh(
                     end = vertex_coords[segment[1]]
                     internal_line_coords = start * (1 - line_weights) + end * line_weights
 
-                    in_order = segment[0] < segment[1]
+                    in_order = _original_order(segment[0], segment[1])
 
                     mesh_nodes = np.vstack((mesh_nodes, internal_line_coords))
                     cells.append(np.array([segment[0], start_coord_pos], dtype=np.int64))
@@ -318,7 +322,7 @@ def _build_mesh(
             cell_markers_ = np.array(cell_markers, dtype=np.int32)
 
             orientations = np.full_like(cell_markers_, 1.0, dtype=np.float64)
-            orientations[cells_[:, 0] > cells_[:, 1]] = -1.0
+            orientations[~_original_order(cells_[:, 0], cells_[:, 1])] = -1.0
 
             assert cell_markers_.shape == orientations.shape
         else:

From 944e3e320203b0e5b2de4fc79b6857ca2c1a63ef Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Wed, 18 Feb 2026 22:56:34 +0100
Subject: [PATCH 16/17] Missed one

---
 src/networks_fenicsx/mesh.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/networks_fenicsx/mesh.py b/src/networks_fenicsx/mesh.py
index 438b42c..c1786fa 100644
--- a/src/networks_fenicsx/mesh.py
+++ b/src/networks_fenicsx/mesh.py
@@ -384,7 +384,7 @@ def _original_order(a, b):
         e_geo = mesh.entities_to_geometry(self.mesh, tdim, e_idx)
 
         global_input = self.mesh.geometry.input_global_indices
-        in_order = global_input[e_geo[:, 0]] < global_input[e_geo[:, 1]]
+        in_order = _original_order(global_input[e_geo[:, 0]], global_input[e_geo[:, 1]])
 
         # Four cases that might arise (per edge), with naming i := input_in_order, n := in_order:
         # 1) i & n:

From 106168caf49f36148af89591e2f73876991728ad Mon Sep 17 00:00:00 2001
From: schnellerhase <56360279+schnellerhase@users.noreply.github.com>
Date: Thu, 19 Feb 2026 12:58:07 +0100
Subject: [PATCH 17/17] Simplify (remove in_order checks)

---
 src/networks_fenicsx/mesh.py | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/src/networks_fenicsx/mesh.py b/src/networks_fenicsx/mesh.py
index c1786fa..f77683c 100644
--- a/src/networks_fenicsx/mesh.py
+++ b/src/networks_fenicsx/mesh.py
@@ -282,7 +282,6 @@ def _original_order(a, b):
                     cell_markers.append(edge_coloring[(segment[0], segment[1])])
                     start = vertex_coords[segment[0]]
                     end = vertex_coords[segment[1]]
-                    in_order = _original_order(segment[0], segment[1])
             else:
                 for segment in cells_array:
                     start_coord_pos = mesh_nodes.shape[0]
@@ -290,8 +289,6 @@ def _original_order(a, b):
                     end = vertex_coords[segment[1]]
                     internal_line_coords = start * (1 - line_weights) + end * line_weights
 
-                    in_order = _original_order(segment[0], segment[1])
-
                     mesh_nodes = np.vstack((mesh_nodes, internal_line_coords))
                     cells.append(np.array([segment[0], start_coord_pos], dtype=np.int64))
                     segment_connectivity = (