"""Graph and mesh topology utilities for finite element computations.
This module provides data structures and algorithms for handling:
- 2D simplex mesh connectivity (triangles to unique edges and orientations).
- Compressed index data representations for ragged/jagged arrays (CSR-like indexing).
- Upward mesh adjacencies (vertex-to-element reverse mappings).
- Directed and undirected graph connectivity checks and connected component coloring.
- Coordinate extraction for triangle elements (with optional Numba acceleration).
"""
from typing import Annotated, Optional, Tuple
import numpy as np
from numpy.typing import NDArray
# Tous les imports pyxfem regroupés ensemble
from pyxfem.core.config import (
FASTMATH,
USE_NUMBA,
FloatArray,
IndexArray,
IndexType,
npfloat,
)
# Numba n'est importé que si explicitement activé
if USE_NUMBA:
import numba
T2V = Annotated[IndexArray, "Shape(N, 3)"] # type for triangle to vertex connectivity table
E2V = Annotated[IndexArray, "Shape(N, 2)"] # type for edge to vertex connectivity table
T2E = T2V# type for triangle to edge connectivity table
DownwardAdjacencies = Annotated[IndexArray, "Shape(N, K)"]
[docs]
def create_edges_from_tris(t2v: T2V) -> Tuple[E2V, T2E, T2E]:
"""Extract unique edges and their local orientations from triangle connectivities.
Parameters
----------
t2v : NDArray of shape (nt, 3)
Triangle-to-vertex connectivity matrix, where each row contains
the indices of the three vertices defining a triangle.
Returns
-------
e2v : NDArray of shape (ne, 2)
Array of unique edges bounding the triangles. Each row `[vid0, vid1]`
contains vertex IDs ordered such that `vid0 < vid1`.
t2e : NDArray of shape (nt, 3)
Triangle-to-edge connectivity table. For each triangle, lists the global
IDs of its 3 unique bounding edges.
t2eflip : NDArray of shape (nt, 3)
Edge orientation flag matrix containing binary values (0 or 1).
A value of 1 indicates that the local triangle edge directed as
`v_i -> v_{i+1}` is reversed relative to the canonical global edge in `e2v`.
"""
t2ev = t2v[:, np.array([[0, 1], [1, 2], [2, 0]])]
e2v: E2V = t2ev.reshape(-1, 2)
asort = e2v.argsort(axis=1)
e2v, emap = np.unique(
np.take_along_axis(e2v, asort, 1), axis=0, return_inverse=True
)
t2e: T2E = emap.reshape((t2v.shape[0], 3))
t2eflip: T2E = np.where(t2ev[:, :, 0] > t2ev[:, :, 1], 1, 0)
return e2v, t2e, t2eflip
[docs]
class CompressedDataIndex:
"""Compressed Sparse Row (CSR) style storage for ragged 1D arrays.
Stores a list of variable-length sub-arrays within a single contiguous
data array, indexed by offset pointers.
Parameters
----------
index : NDArray
Array of starting offsets of length `N + 1`, where `index[k]`
to `index[k + 1]` defines the slice for entry `k`.
data : NDArray
Contiguous 1D buffer containing the concatenated elements.
check : bool, optional
If True, validates that all index bounds fall within `data` (default is False).
Raises
------
ValueError
If `check=True` and any index value is negative or exceeds the data length.
"""
def __init__(self, index: IndexArray, data: NDArray, check=False):
if check:
if not np.all(0 <= np.max(index) < data.shape[0]):
raise ValueError
self.index = index
self.data = data
[docs]
def ndata_perentry(self):
"""Compute the number of items stored in each entry.
Returns
-------
NDArray
1D array of length `N` containing the count of elements for each entry.
"""
return self.index[1:] - self.index[:-1]
def __getitem__(self, key):
"""Retrieve the slice of elements corresponding to a specific entry index.
Parameters
----------
key : int
Entry index.
Returns
-------
NDArray
View of the slice of `data` associated with entry `key`.
"""
return self.data[self.index[key] : self.index[key + 1]]
[docs]
def build_compressed_data_index_from_list_of_list(
lofl, check=False
) -> CompressedDataIndex:
"""Construct a CompressedDataIndex instance from a nested Python list.
Parameters
----------
lofl : list of list of int
List of variable-length integer sequences.
check : bool, optional
Whether to validate bounds on creation (default is False).
Returns
-------
CompressedDataIndex
Initialized compressed sparse index structure.
"""
nb_index = len(lofl)
index = np.zeros(nb_index+1, dtype = IndexType)
nbs = np.fromiter((len(i) for i in lofl), dtype = IndexType)
np.cumsum(nbs, dtype = IndexType, out = index[1:])
data = np.concatenate(lofl)
return CompressedDataIndex(index, data, check)
UpwardAdjacencies = CompressedDataIndex
[docs]
def upward_adjacencies(t: NDArray) -> UpwardAdjacencies:
"""Invert downward element-to-vertex adjacencies into vertex-to-element adjacencies.
Given an element connectivity array (e.g., triangle-to-vertex), constructs
the reverse mapping indicating which elements share each vertex.
Parameters
----------
t : NDArray of shape (num_elements, num_nodes_per_element)
Element connectivity table.
Returns
-------
UpwardAdjacencies
Compressed index structure where entry `i` yields the array of
elements sharing vertex `i`.
"""
flatt = t.flatten()
order = np.argsort(flatt)
fo = flatt[order]
fo = np.array(fo, dtype =IndexType)
bc = np.bincount(fo)
#elementstart = np.cumsum(np.hstack((np.array([0], dtype = 'int'),bc)))
element_index = np.array(np.cumsum(np.hstack(
(np.array([0], dtype = IndexType),bc))), dtype = IndexType)
element_data = np.array(order//t.shape[1], dtype = IndexType)
return UpwardAdjacencies(element_index, element_data)
[docs]
def which_in_upward(
upward: UpwardAdjacencies, downward: DownwardAdjacencies
) -> NDArray:
"""Locate the local node index of each vertex within its upward connected elements.
Parameters
----------
upward : UpwardAdjacencies
Upward vertex-to-element adjacencies.
downward : DownwardAdjacencies
Downward element-to-vertex connectivity table.
Returns
-------
NDArray of shape matching `upward.data`
Local column index (0, 1, ..., K-1) of vertex `vid` inside each element
referencing it in the upward adjacency structure.
"""
which_v = np.empty(upward.data.shape, dtype=IndexType)
for vid in range(len(upward.index) - 1):
s = upward.index[vid]
e = upward.index[vid+1]
which_v[s:e] = np.where(downward[upward.data[s:e]] ==vid)[1]
return which_v
[docs]
class Connections:
"""Directed graph connectivity representation using CSR-style indexing.
Parameters
----------
compressed_connections_data : CompressedDataIndex
Underlying compressed data representing the adjacency list.
check : bool, optional
If True, validates graph integrity on initialization (default is True).
Raises
------
ValueError
If `check=True` and the provided structure fails connectivity validation.
"""
[docs]
@staticmethod
def check(compressed_connections_data: CompressedDataIndex) -> bool:
"""Validate if a CompressedDataIndex conforms to a valid directed graph.
Checks dtype consistency, absence of negative vertex indices,
and ensures edge targets do not exceed the vertex count.
Parameters
----------
compressed_connections_data : CompressedDataIndex
Graph adjacency representation to check.
Returns
-------
bool
True if the data forms a valid directed graph, False otherwise.
"""
index = compressed_connections_data.index
if index.dtype != IndexType:
return False
data = compressed_connections_data.data
if np.max(index) > len(data) + 1:
return False
if np.any(data < 0):
return False
if np.any(data > (len(index) - 2)):
return False
return True
def __init__(self, compressed_connections_data:CompressedDataIndex, check = True):
if check :
if not Connections.check(compressed_connections_data):
raise ValueError
self.index = compressed_connections_data.index
self.data = compressed_connections_data.data.astype(IndexType)
[docs]
def isin(self, a: int) -> bool:
"""Check whether vertex index `a` exists in the graph.
Parameters
----------
a : int
Node index to check.
Returns
-------
bool
True if node `a` is within graph bounds, False otherwise.
"""
return 0 <= a < len(self.index - 1)
[docs]
def isconnected(self, a: int, b: int, check=False) -> bool:
"""Check if a directed edge exists from node `a` to node `b` (a -> b).
Parameters
----------
a : int
Source node index.
b : int
Target node index.
check : bool, optional
Whether to verify node existence bounds first (default is False).
Returns
-------
bool
True if a directed edge from `a` to `b` exists, False otherwise.
Raises
------
ValueError
If `check=True` and either node `a` or `b` is out of bounds.
"""
if check:
if not (self.isin(a) and self.isin(b)):
raise ValueError
index_sa = self.index[a]
index_ea = self.index[a+1]
return b in self.data[index_sa : index_ea]
[docs]
class UndirectedConnections:
"""Undirected graph connectivity representation using CSR-style indexing.
Enforces symmetry: if an edge exists from `a -> b`, the reverse edge
`b -> a` must also be present.
Parameters
----------
connect : Connections
Directed connection instance to interpret as undirected.
check : bool, optional
If True, validates reciprocity of edges (default is True).
Raises
------
ValueError
If `check=True` and reciprocity check fails.
"""
[docs]
@staticmethod
def check(connect: Connections) -> bool:
"""Verify that every edge in the graph is reciprocal (if a->b then b->a).
Parameters
----------
connect : Connections
Directed graph connection to inspect.
Returns
-------
bool
True if all connections are strictly undirected, False otherwise.
"""
index = connect.index
conn = connect.data
for i in range(len(index) - 1):
iconn = conn[index[i] : index[i + 1]]
for j in iconn:
jconn = conn[index[j]: index[j+1]]
if not np.isin(i, jconn):
return False
return True
def __init__(self, connect:Connections, check = True):
if check:
if not UndirectedConnections.check(connect):
raise ValueError
self.index = connect.index
self.data = connect.data
[docs]
def isin(self, a: int) -> bool:
"""Check whether vertex index `a` exists in the graph.
Parameters
----------
a : int
Node index.
Returns
-------
bool
True if node `a` is within graph bounds, False otherwise.
"""
return 0 <= a < len(self.index - 1)
[docs]
def isconnected(self, a: int, b: int, check=False) -> bool:
"""Check whether nodes `a` and `b` are connected.
Parameters
----------
a : int
First node index.
b : int
Second node index.
check : bool, optional
Whether to check bounds prior to lookup (default is False).
Returns
-------
bool
True if `a` and `b` share an edge, False otherwise.
Raises
------
ValueError
If `check=True` and either node is out of bounds.
"""
if check:
if not (self.isin(a) and self.isin(b)):
raise ValueError
index_sa = self.index[a]
index_ea = self.index[a+1]
return b in self.data[index_sa : index_ea]
[docs]
def build_undirected_connection(lofl, check=True) -> UndirectedConnections:
"""Build an UndirectedConnections instance directly from a list of adjacency lists.
Parameters
----------
lofl : list of list of int
Adjacency list representation where entry `i` contains neighbors of node `i`.
check : bool, optional
Whether to enforce connection validity and symmetry checks (default is True).
Returns
-------
UndirectedConnections
The resulting undirected graph connectivity structure.
Raises
------
ValueError
If `check=True` and the graph is invalid or asymmetric.
"""
cdata = build_compressed_data_index_from_list_of_list(lofl)
conn = Connections(cdata, check)
return UndirectedConnections(conn, check)
[docs]
def coloring_connection(conn: Connections) -> NDArray:
"""Label weakly connected components in a directed graph.
Performs graph traversal to assign a unique 1-based component identifier
to all reachable vertices.
Parameters
----------
conn : Connections
Directed graph connectivity structure.
Returns
-------
NDArray of int
Array of length `num_nodes` containing 1-based component color indices.
"""
color = 0
index = conn.index
data = conn.data
colors = np.zeros(len(index)-1, dtype =IndexType)
while np.any(colors == 0):
color += 1
queue = [ int(np.nonzero(colors == 0)[0][0])]
while len(queue) >0:
cid = queue.pop()
if colors[cid] < color :
colors[cid] = color
queue += list(data[ index[cid] : index[cid+1]].astype(int))
return colors
[docs]
def coloring_undirected_connection(conn: UndirectedConnections) -> NDArray:
"""Label connected components of an undirected graph with 0-based group IDs.
Parameters
----------
conn : UndirectedConnections
Undirected graph connectivity structure.
Returns
-------
NDArray of int
Array of component color IDs spanning `[0, num_components - 1]`.
"""
color = 0
index = conn.index
data = conn.data
colors = np.zeros(len(index)-1, dtype =IndexType)
while np.any(colors == 0):
color += 1
queue = [ int(np.nonzero(colors == 0)[0][0])]
while len(queue) >0:
cid = queue.pop()
if colors[cid] != color :
colors[cid] = color
queue += list(data[ index[cid] : index[cid+1]].astype(int))
# Passe finale : renumérotation contiguë de 0 à n_colors - 1
# Ce serait mieu, mais c'est a verifier que ca ne casse pas mesh_modifications
#_, contiguous_colors = np.unique(colors, return_inverse=True)
#return contiguous_colors.astype(IndexType)
return colors-1
if USE_NUMBA :