Coverage for pyxfem/core/graph.py: 78%

1"""Graph and mesh topology utilities for finite element computations. 

2 

3This module provides data structures and algorithms for handling: 

4- 2D simplex mesh connectivity (triangles to unique edges and orientations). 

5- Compressed index data representations for ragged/jagged arrays (CSR-like indexing). 

6- Upward mesh adjacencies (vertex-to-element reverse mappings). 

7- Directed and undirected graph connectivity checks and connected component coloring. 

8- Coordinate extraction for triangle elements (with optional Numba acceleration). 

9""" 

10 

11from typing import Annotated, Optional, Tuple 

12 

13import numpy as np 

14from numpy.typing import NDArray 

15 

16# Tous les imports pyxfem regroupés ensemble 

17from pyxfem.core.config import ( 

18 FASTMATH, 

19 USE_NUMBA, 

20 FloatArray, 

21 IndexArray, 

22 IndexType, 

23 npfloat, 

24) 

25 

26# Numba n'est importé que si explicitement activé 

27if USE_NUMBA: 

28 import numba 

29 

30T2V = Annotated[IndexArray, "Shape(N, 3)"] # type for triangle to vertex connectivity table 

31E2V = Annotated[IndexArray, "Shape(N, 2)"] # type for edge to vertex connectivity table 

32T2E = T2V# type for triangle to edge connectivity table 

33 

34DownwardAdjacencies = Annotated[IndexArray, "Shape(N, K)"] 

35 

36 

37def create_edges_from_tris(t2v: T2V) -> Tuple[E2V, T2E, T2E]: 

38 """Extract unique edges and their local orientations from triangle connectivities. 

39 

40 Parameters 

41 ---------- 

42 t2v : NDArray of shape (nt, 3) 

43 Triangle-to-vertex connectivity matrix, where each row contains 

44 the indices of the three vertices defining a triangle. 

45 

46 Returns 

47 ------- 

48 e2v : NDArray of shape (ne, 2) 

49 Array of unique edges bounding the triangles. Each row `[vid0, vid1]` 

50 contains vertex IDs ordered such that `vid0 < vid1`. 

51 t2e : NDArray of shape (nt, 3) 

52 Triangle-to-edge connectivity table. For each triangle, lists the global 

53 IDs of its 3 unique bounding edges. 

54 t2eflip : NDArray of shape (nt, 3) 

55 Edge orientation flag matrix containing binary values (0 or 1). 

56 A value of 1 indicates that the local triangle edge directed as 

57 `v_i -> v_{i+1}` is reversed relative to the canonical global edge in `e2v`. 

58 """ 

59 t2ev = t2v[:, np.array([[0, 1], [1, 2], [2, 0]])] 

60 e2v: E2V = t2ev.reshape(-1, 2) 

61 asort = e2v.argsort(axis=1) 

62 e2v, emap = np.unique( 

63 np.take_along_axis(e2v, asort, 1), axis=0, return_inverse=True 

64 ) 

65 t2e: T2E = emap.reshape((t2v.shape[0], 3)) 

66 t2eflip: T2E = np.where(t2ev[:, :, 0] > t2ev[:, :, 1], 1, 0) 

67 return e2v, t2e, t2eflip 

68 

69 

70class CompressedDataIndex: 

71 """Compressed Sparse Row (CSR) style storage for ragged 1D arrays. 

72 

73 Stores a list of variable-length sub-arrays within a single contiguous 

74 data array, indexed by offset pointers. 

75 

76 Parameters 

77 ---------- 

78 index : NDArray 

79 Array of starting offsets of length `N + 1`, where `index[k]` 

80 to `index[k + 1]` defines the slice for entry `k`. 

81 data : NDArray 

82 Contiguous 1D buffer containing the concatenated elements. 

83 check : bool, optional 

84 If True, validates that all index bounds fall within `data` (default is False). 

85 

86 Raises 

87 ------ 

88 ValueError 

89 If `check=True` and any index value is negative or exceeds the data length. 

90 """ 

91 

92 def __init__(self, index: IndexArray, data: NDArray, check=False): 

93 if check: 

94 if not np.all(0 <= np.max(index) < data.shape[0]): 

95 raise ValueError 

96 self.index = index 

97 self.data = data 

98 

99 def ndata_perentry(self): 

100 """Compute the number of items stored in each entry. 

101 

102 Returns 

103 ------- 

104 NDArray 

105 1D array of length `N` containing the count of elements for each entry. 

106 """ 

107 return self.index[1:] - self.index[:-1] 

108 

109 def __getitem__(self, key): 

110 """Retrieve the slice of elements corresponding to a specific entry index. 

111 

112 Parameters 

113 ---------- 

114 key : int 

115 Entry index. 

116 

117 Returns 

118 ------- 

119 NDArray 

120 View of the slice of `data` associated with entry `key`. 

121 """ 

122 return self.data[self.index[key] : self.index[key + 1]] 

123 

124 

125def build_compressed_data_index_from_list_of_list( 

126 lofl, check=False 

127) -> CompressedDataIndex: 

128 """Construct a CompressedDataIndex instance from a nested Python list. 

129 

130 Parameters 

131 ---------- 

132 lofl : list of list of int 

133 List of variable-length integer sequences. 

134 check : bool, optional 

135 Whether to validate bounds on creation (default is False). 

136 

137 Returns 

138 ------- 

139 CompressedDataIndex 

140 Initialized compressed sparse index structure. 

141 """ 

142 nb_index = len(lofl) 

143 index = np.zeros(nb_index+1, dtype = IndexType) 

144 nbs = np.fromiter((len(i) for i in lofl), dtype = IndexType) 

145 np.cumsum(nbs, dtype = IndexType, out = index[1:]) 

146 data = np.concatenate(lofl) 

147 return CompressedDataIndex(index, data, check) 

148 

149UpwardAdjacencies = CompressedDataIndex 

150 

151 

152def upward_adjacencies(t: NDArray) -> UpwardAdjacencies: 

153 """Invert downward element-to-vertex adjacencies into vertex-to-element adjacencies. 

154 

155 Given an element connectivity array (e.g., triangle-to-vertex), constructs 

156 the reverse mapping indicating which elements share each vertex. 

157 

158 Parameters 

159 ---------- 

160 t : NDArray of shape (num_elements, num_nodes_per_element) 

161 Element connectivity table. 

162 

163 Returns 

164 ------- 

165 UpwardAdjacencies 

166 Compressed index structure where entry `i` yields the array of 

167 elements sharing vertex `i`. 

168 """ 

169 flatt = t.flatten() 

170 order = np.argsort(flatt) 

171 fo = flatt[order] 

172 fo = np.array(fo, dtype =IndexType) 

173 bc = np.bincount(fo) 

174 #elementstart = np.cumsum(np.hstack((np.array([0], dtype = 'int'),bc))) 

175 element_index = np.array(np.cumsum(np.hstack( 

176 (np.array([0], dtype = IndexType),bc))), dtype = IndexType) 

177 element_data = np.array(order//t.shape[1], dtype = IndexType) 

178 return UpwardAdjacencies(element_index, element_data) 

179 

180def which_in_upward( 

181 upward: UpwardAdjacencies, downward: DownwardAdjacencies 

182) -> NDArray: 

183 """Locate the local node index of each vertex within its upward connected elements. 

184 

185 Parameters 

186 ---------- 

187 upward : UpwardAdjacencies 

188 Upward vertex-to-element adjacencies. 

189 downward : DownwardAdjacencies 

190 Downward element-to-vertex connectivity table. 

191 

192 Returns 

193 ------- 

194 NDArray of shape matching `upward.data` 

195 Local column index (0, 1, ..., K-1) of vertex `vid` inside each element 

196 referencing it in the upward adjacency structure. 

197 """ 

198 which_v = np.empty(upward.data.shape, dtype=IndexType) 

199 for vid in range(len(upward.index) - 1): 

200 s = upward.index[vid] 

201 e = upward.index[vid+1] 

202 which_v[s:e] = np.where(downward[upward.data[s:e]] ==vid)[1] 

203 return which_v 

204 

205 

206class Connections: 

207 """Directed graph connectivity representation using CSR-style indexing. 

208 

209 Parameters 

210 ---------- 

211 compressed_connections_data : CompressedDataIndex 

212 Underlying compressed data representing the adjacency list. 

213 check : bool, optional 

214 If True, validates graph integrity on initialization (default is True). 

215 

216 Raises 

217 ------ 

218 ValueError 

219 If `check=True` and the provided structure fails connectivity validation. 

220 """ 

221 

222 @staticmethod 

223 def check(compressed_connections_data: CompressedDataIndex) -> bool: 

224 """Validate if a CompressedDataIndex conforms to a valid directed graph. 

225 

226 Checks dtype consistency, absence of negative vertex indices, 

227 and ensures edge targets do not exceed the vertex count. 

228 

229 Parameters 

230 ---------- 

231 compressed_connections_data : CompressedDataIndex 

232 Graph adjacency representation to check. 

233 

234 Returns 

235 ------- 

236 bool 

237 True if the data forms a valid directed graph, False otherwise. 

238 """ 

239 index = compressed_connections_data.index 

240 if index.dtype != IndexType: 

241 return False 

242 data = compressed_connections_data.data 

243 if np.max(index) > len(data) + 1: 

244 return False 

245 if np.any(data < 0): 

246 return False 

247 if np.any(data > (len(index) - 2)): 

248 return False 

249 return True 

250 

251 def __init__(self, compressed_connections_data:CompressedDataIndex, check = True): 

252 if check : 

253 if not Connections.check(compressed_connections_data): 

254 raise ValueError 

255 self.index = compressed_connections_data.index 

256 self.data = compressed_connections_data.data.astype(IndexType) 

257 

258 def isin(self, a: int) -> bool: 

259 """Check whether vertex index `a` exists in the graph. 

260 

261 Parameters 

262 ---------- 

263 a : int 

264 Node index to check. 

265 

266 Returns 

267 ------- 

268 bool 

269 True if node `a` is within graph bounds, False otherwise. 

270 """ 

271 return 0 <= a < len(self.index - 1) 

272 

273 def isconnected(self, a: int, b: int, check=False) -> bool: 

274 """Check if a directed edge exists from node `a` to node `b` (a -> b). 

275 

276 Parameters 

277 ---------- 

278 a : int 

279 Source node index. 

280 b : int 

281 Target node index. 

282 check : bool, optional 

283 Whether to verify node existence bounds first (default is False). 

284 

285 Returns 

286 ------- 

287 bool 

288 True if a directed edge from `a` to `b` exists, False otherwise. 

289 

290 Raises 

291 ------ 

292 ValueError 

293 If `check=True` and either node `a` or `b` is out of bounds. 

294 """ 

295 if check: 

296 if not (self.isin(a) and self.isin(b)): 

297 raise ValueError 

298 index_sa = self.index[a] 

299 index_ea = self.index[a+1] 

300 return b in self.data[index_sa : index_ea] 

301 

302class UndirectedConnections: 

303 """Undirected graph connectivity representation using CSR-style indexing. 

304 

305 Enforces symmetry: if an edge exists from `a -> b`, the reverse edge 

306 `b -> a` must also be present. 

307 

308 Parameters 

309 ---------- 

310 connect : Connections 

311 Directed connection instance to interpret as undirected. 

312 check : bool, optional 

313 If True, validates reciprocity of edges (default is True). 

314 

315 Raises 

316 ------ 

317 ValueError 

318 If `check=True` and reciprocity check fails. 

319 """ 

320 

321 @staticmethod 

322 def check(connect: Connections) -> bool: 

323 """Verify that every edge in the graph is reciprocal (if a->b then b->a). 

324 

325 Parameters 

326 ---------- 

327 connect : Connections 

328 Directed graph connection to inspect. 

329 

330 Returns 

331 ------- 

332 bool 

333 True if all connections are strictly undirected, False otherwise. 

334 """ 

335 index = connect.index 

336 conn = connect.data 

337 for i in range(len(index) - 1): 

338 iconn = conn[index[i] : index[i + 1]] 

339 for j in iconn: 

340 jconn = conn[index[j]: index[j+1]] 

341 if not np.isin(i, jconn): 

342 return False 

343 return True 

344 

345 def __init__(self, connect:Connections, check = True): 

346 if check: 

347 if not UndirectedConnections.check(connect): 

348 raise ValueError 

349 self.index = connect.index 

350 self.data = connect.data 

351 

352 def isin(self, a: int) -> bool: 

353 """Check whether vertex index `a` exists in the graph. 

354 

355 Parameters 

356 ---------- 

357 a : int 

358 Node index. 

359 

360 Returns 

361 ------- 

362 bool 

363 True if node `a` is within graph bounds, False otherwise. 

364 """ 

365 return 0 <= a < len(self.index - 1) 

366 

367 def isconnected(self, a: int, b: int, check=False) -> bool: 

368 """Check whether nodes `a` and `b` are connected. 

369 

370 Parameters 

371 ---------- 

372 a : int 

373 First node index. 

374 b : int 

375 Second node index. 

376 check : bool, optional 

377 Whether to check bounds prior to lookup (default is False). 

378 

379 Returns 

380 ------- 

381 bool 

382 True if `a` and `b` share an edge, False otherwise. 

383 

384 Raises 

385 ------ 

386 ValueError 

387 If `check=True` and either node is out of bounds. 

388 """ 

389 if check: 

390 if not (self.isin(a) and self.isin(b)): 

391 raise ValueError 

392 index_sa = self.index[a] 

393 index_ea = self.index[a+1] 

394 return b in self.data[index_sa : index_ea] 

395 

396def build_undirected_connection(lofl, check=True) -> UndirectedConnections: 

397 """Build an UndirectedConnections instance directly from a list of adjacency lists. 

398 

399 Parameters 

400 ---------- 

401 lofl : list of list of int 

402 Adjacency list representation where entry `i` contains neighbors of node `i`. 

403 check : bool, optional 

404 Whether to enforce connection validity and symmetry checks (default is True). 

405 

406 Returns 

407 ------- 

408 UndirectedConnections 

409 The resulting undirected graph connectivity structure. 

410 

411 Raises 

412 ------ 

413 ValueError 

414 If `check=True` and the graph is invalid or asymmetric. 

415 """ 

416 cdata = build_compressed_data_index_from_list_of_list(lofl) 

417 conn = Connections(cdata, check) 

418 return UndirectedConnections(conn, check) 

419 

420 

421def coloring_connection(conn: Connections) -> NDArray: 

422 """Label weakly connected components in a directed graph. 

423 

424 Performs graph traversal to assign a unique 1-based component identifier 

425 to all reachable vertices. 

426 

427 Parameters 

428 ---------- 

429 conn : Connections 

430 Directed graph connectivity structure. 

431 

432 Returns 

433 ------- 

434 NDArray of int 

435 Array of length `num_nodes` containing 1-based component color indices. 

436 """ 

437 color = 0 

438 index = conn.index 

439 data = conn.data 

440 colors = np.zeros(len(index)-1, dtype =IndexType) 

441 while np.any(colors == 0): 

442 color += 1 

443 queue = [ int(np.nonzero(colors == 0)[0][0])] 

444 while len(queue) >0: 

445 cid = queue.pop() 

446 if colors[cid] < color : 

447 colors[cid] = color 

448 queue += list(data[ index[cid] : index[cid+1]].astype(int)) 

449 return colors 

450 

451 

452def coloring_undirected_connection(conn: UndirectedConnections) -> NDArray: 

453 """Label connected components of an undirected graph with 0-based group IDs. 

454 

455 Parameters 

456 ---------- 

457 conn : UndirectedConnections 

458 Undirected graph connectivity structure. 

459 

460 Returns 

461 ------- 

462 NDArray of int 

463 Array of component color IDs spanning `[0, num_components - 1]`. 

464 """ 

465 color = 0 

466 index = conn.index 

467 data = conn.data 

468 colors = np.zeros(len(index)-1, dtype =IndexType) 

469 while np.any(colors == 0): 

470 color += 1 

471 queue = [ int(np.nonzero(colors == 0)[0][0])] 

472 while len(queue) >0: 

473 cid = queue.pop() 

474 if colors[cid] != color : 

475 colors[cid] = color 

476 queue += list(data[ index[cid] : index[cid+1]].astype(int)) 

477 

478 # Passe finale : renumérotation contiguë de 0 à n_colors - 1 

479 # Ce serait mieu, mais c'est a verifier que ca ne casse pas mesh_modifications 

480 #_, contiguous_colors = np.unique(colors, return_inverse=True) 

481 #return contiguous_colors.astype(IndexType) 

482 

483 return colors-1 

484 

485 

486def extract_txy( 

487 xy: FloatArray, t2v: T2V, out: Optional[FloatArray] = None 

488) -> FloatArray: 

489 """Extract spatial coordinates for all nodes in each triangle element. 

490 

491 Parameters 

492 ---------- 

493 xy : NDArray of shape (num_nodes, dim) 

494 Global vertex coordinate table. 

495 t2v : NDArray of shape (num_elements, 3) 

496 Triangle-to-vertex connectivity table. 

497 out : NDArray of shape (num_elements, 3, dim), optional 

498 Preallocated output buffer. If None, an array is allocated. 

499 

500 Returns 

501 ------- 

502 NDArray of shape (num_elements, 3, dim) 

503 Coordinates of the 3 vertices for each triangle element. 

504 """ 

505 if out is None: 

506 shape = (t2v.shape) + xy.shape[-1:] 

507 out = np.empty(shape, npfloat) 

508 if USE_NUMBA : 

509 # print('bla') 

510 extract_txy_numba(xy, t2v, out) 

511 else : 

512 np.take(xy, t2v, axis =0, out =out) 

513 return out 

514 

515if USE_NUMBA : 

516 @numba.njit(parallel=True, cache=True, fastmath= FASTMATH) 

517 def extract_txy_numba(xy, t2v, out_txy): 

518 """Accelerated parallel kernel for triangle vertex coordinate extraction. 

519 

520 Parameters 

521 ---------- 

522 xy : NDArray of shape (N, dim) 

523 Node coordinate table. 

524 t2v : NDArray of shape (num_elements, 3) 

525 Triangle-to-vertex connectivity table. 

526 out_txy : NDArray of shape (num_elements, 3, dim) 

527 Preallocated array written in place. 

528 """ 

529 N = t2v.shape[0] 

530 for ie in numba.prange(N): #pylint: disable=not-an-iterable 

531 v0, v1, v2 = t2v[ie] 

532 out_txy[ie, 0, :] = xy[v0] 

533 out_txy[ie, 1, :] = xy[v1] 

534 out_txy[ie, 2, :] = xy[v2]