@@ -396,14 +396,18 @@ def compute_circularity(loop_polydata):
396396 if not isinstance (mesh , type (None )):
397397 global_nodes = mesh .points
398398 global_node_tree = cKDTree (global_nodes )
399- global_elements = mesh .cell_connectivity .reshape (- 1 , 4 )
400- global_elements = numpy .sort (global_elements , axis = 1 )
401- tet_faces = []
402- for i in tqdm .trange (global_elements .shape [0 ], desc = "Building tetrahedral faces" , leave = False ):
403- for j in range (4 ):
404- idx = set (list (range (4 ))) - set ([j ])
405- tet_faces .append (global_elements [i , list (idx )])
406- tet_face_tree = cKDTree (tet_faces )
399+ # Build a robust face->cell index using canonical keys (sorted tuples)
400+ face_to_cell = {}
401+ n_cells = mesh .n_cells
402+ for i in tqdm .trange (n_cells , desc = "Indexing cell faces" , leave = False ):
403+ cell = mesh .GetCell (i )
404+ nfaces = cell .GetNumberOfFaces ()
405+ for j in range (nfaces ):
406+ face = cell .GetFace (j )
407+ npts = face .GetNumberOfPoints ()
408+ # Canonicalize face nodes to a sorted tuple so orientation/order doesn't matter
409+ key = tuple (sorted (face .GetPointId (k ) for k in range (npts )))
410+ face_to_cell [key ] = i
407411 #for i, cap in enumerate(iscap):
408412 # if not cap == 1:
409413 # walls.append(faces[i])
@@ -416,18 +420,33 @@ def compute_circularity(loop_polydata):
416420 wall_cells = surface .extract_cells (walls [i ])
417421 wall_surface = wall_cells .extract_surface ()
418422 if not isinstance (mesh , type (None )):
419- wall_surface .point_data ["GlobalNodeID" ] = numpy .zeros (wall_surface .n_points , dtype = int )
420- wall_surface .cell_data ['GlobalElementID' ] = numpy .zeros (wall_surface .n_cells , dtype = int )
421- wall_surface .cell_data ['ModelFaceID' ] = numpy .ones (wall_surface .n_cells , dtype = int )
423+ wall_surface .point_data ["GlobalNodeID" ] = numpy .zeros (wall_surface .n_points , dtype = numpy . int32 )
424+ wall_surface .cell_data ['GlobalElementID' ] = numpy .zeros (wall_surface .n_cells , dtype = numpy . int32 )
425+ wall_surface .cell_data ['ModelFaceID' ] = numpy .ones (wall_surface .n_cells , dtype = numpy . int32 )
422426 _ , indices = global_node_tree .query (wall_surface .points )
423427 wall_surface .point_data ["GlobalNodeID" ] = indices .astype (int )
424428 # Assign Global Element IDs
425429 wall_faces = wall_surface .point_data ["GlobalNodeID" ][wall_surface .faces ]
426- wall_faces = wall_faces .reshape (- 1 , 4 )[:, 1 :]
427- wall_faces = numpy .sort (wall_faces , axis = 1 )
428- _ , indices = tet_face_tree .query (wall_faces )
429- wall_surface .cell_data ["GlobalElementID" ] = indices // 4
430- wall_surface .cell_data ["GlobalElementID" ] = wall_surface .cell_data ["GlobalElementID" ].astype (numpy .int32 )
430+ wall_faces = wall_faces .reshape (- 1 , 4 )[:, 1 :].tolist ()
431+ elem_ids = []
432+ for face in wall_faces :
433+ # Use the same canonical key used to build the map
434+ elem_ids .append (face_to_cell [tuple (sorted (face ))])
435+ wall_surface .cell_data ["GlobalElementID" ] = numpy .array (elem_ids , dtype = numpy .int32 )
436+ #wall_faces = numpy.sort(wall_faces, axis=1)
437+ #dists, indices = tet_face_tree.query(wall_faces)
438+ #if not numpy.all(numpy.isclose(dists, 0.0)):
439+ # # Identify a small sample of mismatches for debugging
440+ # bad_idx = numpy.where(~numpy.isclose(dists, 0.0))[0]
441+ # sample = bad_idx[:5]
442+ # examples = wall_faces[sample]
443+ # raise ValueError(
444+ # f"Failed to map all wall surface faces to volume mesh faces: {bad_idx.size} mismatches. "
445+ # f"Example face node triples (GlobalNodeID) that failed exact match: {examples.tolist()}"
446+ # )
447+ #wall_surface.cell_data["GlobalElementID"] = indices // 4
448+ #wall_surface.cell_data["GlobalElementID"] = wall_surface.cell_data["GlobalElementID"].astype(numpy.int32)
449+
431450 boundaries = wall_surface .extract_feature_edges (boundary_edges = True , manifold_edges = False ,
432451 feature_edges = False , non_manifold_edges = False )
433452 boundaries = boundaries .split_bodies ()
@@ -444,19 +463,31 @@ def compute_circularity(loop_polydata):
444463 cap_cells = surface .extract_cells (face_cap )
445464 cap_surface = cap_cells .extract_surface ()
446465 if not isinstance (mesh , type (None )):
447- cap_surface .point_data ["GlobalNodeID" ] = numpy .zeros (cap_surface .n_points , dtype = int )
448- cap_surface .cell_data ["GlobalElementID" ] = numpy .zeros (cap_surface .n_cells , dtype = int )
449- cap_surface .cell_data ["ModelFaceID" ] = numpy .ones (cap_surface .n_cells , dtype = int ) * ( i + 2 )
466+ cap_surface .point_data ["GlobalNodeID" ] = numpy .zeros (cap_surface .n_points , dtype = numpy . int32 )
467+ cap_surface .cell_data ["GlobalElementID" ] = numpy .zeros (cap_surface .n_cells , dtype = numpy . int32 )
468+ cap_surface .cell_data ["ModelFaceID" ] = numpy .ones (cap_surface .n_cells , dtype = numpy . int32 )
450469 # Assign Global Node IDs
451470 _ , indices = global_node_tree .query (cap_surface .points )
452471 cap_surface .point_data ["GlobalNodeID" ] = indices .astype (int )
453472 # Assign Global Element IDs
454473 cap_faces = cap_surface .point_data ["GlobalNodeID" ][cap_surface .faces ]
455- cap_faces = cap_faces .reshape (- 1 , 4 )[:, 1 :]
456- cap_faces = numpy .sort (cap_faces , axis = 1 )
457- _ , indices = tet_face_tree .query (cap_faces )
458- cap_surface .cell_data ["GlobalElementID" ] = indices // 4
459- cap_surface .cell_data ["GlobalElementID" ] = cap_surface .cell_data ["GlobalElementID" ].astype (numpy .int32 )
474+ cap_faces = cap_faces .reshape (- 1 , 4 )[:, 1 :].tolist ()
475+ elem_ids = []
476+ for face in cap_faces :
477+ elem_ids .append (face_to_cell [tuple (sorted (face ))])
478+ cap_surface .cell_data ["GlobalElementID" ] = numpy .array (elem_ids , dtype = numpy .int32 )
479+ #cap_faces = numpy.sort(cap_faces, axis=1)
480+ #dists, indices = tet_face_tree.query(cap_faces)
481+ #if not numpy.all(numpy.isclose(dists, 0.0)):
482+ # bad_idx = numpy.where(~numpy.isclose(dists, 0.0))[0]
483+ # sample = bad_idx[:5]
484+ # examples = cap_faces[sample]
485+ # raise ValueError(
486+ # f"Failed to map all cap surface faces to volume mesh faces: {bad_idx.size} mismatches. "
487+ # f"Example face node triples (GlobalNodeID) that failed exact match: {examples.tolist()}"
488+ # )
489+ #cap_surface.cell_data["GlobalElementID"] = indices // 4
490+ #cap_surface.cell_data["GlobalElementID"] = cap_surface.cell_data["GlobalElementID"].astype(numpy.int32)
460491 boundaries = cap_surface .extract_feature_edges (boundary_edges = True , manifold_edges = False ,
461492 feature_edges = False , non_manifold_edges = False )
462493 boundaries = boundaries .split_bodies ()
@@ -482,11 +513,23 @@ def compute_circularity(loop_polydata):
482513 lumen_surface .point_data ["GlobalNodeID" ] = indices .astype (int )
483514 # Assign Global Element IDs
484515 lumen_faces = lumen_surface .point_data ["GlobalNodeID" ][lumen_surface .faces ]
485- lumen_faces = lumen_faces .reshape (- 1 , 4 )[:, 1 :]
486- lumen_faces = numpy .sort (lumen_faces , axis = 1 )
487- _ , indices = tet_face_tree .query (lumen_faces )
488- lumen_surface .cell_data ["GlobalElementID" ] = indices // 4
489- lumen_surface .cell_data ["GlobalElementID" ] = lumen_surface .cell_data ["GlobalElementID" ].astype (numpy .int32 )
516+ lumen_faces = lumen_faces .reshape (- 1 , 4 )[:, 1 :].tolist ()
517+ elem_ids = []
518+ for face in lumen_faces :
519+ elem_ids .append (face_to_cell [tuple (sorted (face ))])
520+ lumen_surface .cell_data ["GlobalElementID" ] = numpy .array (elem_ids , dtype = numpy .int32 )
521+ #lumen_faces = numpy.sort(lumen_faces, axis=1)
522+ #dists, indices = tet_face_tree.query(lumen_faces)
523+ #if not numpy.all(numpy.isclose(dists, 0.0)):
524+ # bad_idx = numpy.where(~numpy.isclose(dists, 0.0))[0]
525+ # sample = bad_idx[:5]
526+ # examples = lumen_faces[sample]
527+ # raise ValueError(
528+ # f"Failed to map all lumen surface faces to volume mesh faces: {bad_idx.size} mismatches. "
529+ # f"Example face node triples (GlobalNodeID) that failed exact match: {examples.tolist()}"
530+ # )
531+ #lumen_surface.cell_data["GlobalElementID"] = indices // 4
532+ #lumen_surface.cell_data["GlobalElementID"] = lumen_surface.cell_data["GlobalElementID"].astype(numpy.int32)
490533 boundaries = lumen_surface .extract_feature_edges (boundary_edges = True , manifold_edges = False ,
491534 feature_edges = False , non_manifold_edges = False )
492535 boundaries = boundaries .split_bodies ()
0 commit comments