Skip to content

Commit c8146fc

Browse files
zasextonzasexton
committed
dicts for element mapping
1 parent adc25d1 commit c8146fc

File tree

1 file changed

+72
-56
lines changed

1 file changed

+72
-56
lines changed

‎svv/simulation/utils/extract_faces.py‎

Lines changed: 72 additions & 56 deletions
Original file line numberDiff line numberDiff line change
@@ -396,16 +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-
# Build each face as a sorted triple of node IDs for exact matching
405-
face_nodes = numpy.delete(global_elements[i], j)
406-
face_nodes = numpy.sort(face_nodes)
407-
tet_faces.append(face_nodes)
408-
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
409411
#for i, cap in enumerate(iscap):
410412
# if not cap == 1:
411413
# walls.append(faces[i])
@@ -418,27 +420,33 @@ def compute_circularity(loop_polydata):
418420
wall_cells = surface.extract_cells(walls[i])
419421
wall_surface = wall_cells.extract_surface()
420422
if not isinstance(mesh, type(None)):
421-
wall_surface.point_data["GlobalNodeID"] = numpy.zeros(wall_surface.n_points, dtype=int)
422-
wall_surface.cell_data['GlobalElementID'] = numpy.zeros(wall_surface.n_cells, dtype=int)
423-
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)
424426
_, indices = global_node_tree.query(wall_surface.points)
425427
wall_surface.point_data["GlobalNodeID"] = indices.astype(int)
426428
# Assign Global Element IDs
427429
wall_faces = wall_surface.point_data["GlobalNodeID"][wall_surface.faces]
428-
wall_faces = wall_faces.reshape(-1, 4)[:, 1:]
429-
wall_faces = numpy.sort(wall_faces, axis=1)
430-
dists, indices = tet_face_tree.query(wall_faces)
431-
if not numpy.all(numpy.isclose(dists, 0.0)):
432-
# Identify a small sample of mismatches for debugging
433-
bad_idx = numpy.where(~numpy.isclose(dists, 0.0))[0]
434-
sample = bad_idx[:5]
435-
examples = wall_faces[sample]
436-
raise ValueError(
437-
f"Failed to map all wall surface faces to volume mesh faces: {bad_idx.size} mismatches. "
438-
f"Example face node triples (GlobalNodeID) that failed exact match: {examples.tolist()}"
439-
)
440-
wall_surface.cell_data["GlobalElementID"] = indices // 4
441-
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+
442450
boundaries = wall_surface.extract_feature_edges(boundary_edges=True, manifold_edges=False,
443451
feature_edges=False, non_manifold_edges=False)
444452
boundaries = boundaries.split_bodies()
@@ -455,27 +463,31 @@ def compute_circularity(loop_polydata):
455463
cap_cells = surface.extract_cells(face_cap)
456464
cap_surface = cap_cells.extract_surface()
457465
if not isinstance(mesh, type(None)):
458-
cap_surface.point_data["GlobalNodeID"] = numpy.zeros(cap_surface.n_points, dtype=int)
459-
cap_surface.cell_data["GlobalElementID"] = numpy.zeros(cap_surface.n_cells, dtype=int)
460-
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)
461469
# Assign Global Node IDs
462470
_, indices = global_node_tree.query(cap_surface.points)
463471
cap_surface.point_data["GlobalNodeID"] = indices.astype(int)
464472
# Assign Global Element IDs
465473
cap_faces = cap_surface.point_data["GlobalNodeID"][cap_surface.faces]
466-
cap_faces = cap_faces.reshape(-1, 4)[:, 1:]
467-
cap_faces = numpy.sort(cap_faces, axis=1)
468-
dists, indices = tet_face_tree.query(cap_faces)
469-
if not numpy.all(numpy.isclose(dists, 0.0)):
470-
bad_idx = numpy.where(~numpy.isclose(dists, 0.0))[0]
471-
sample = bad_idx[:5]
472-
examples = cap_faces[sample]
473-
raise ValueError(
474-
f"Failed to map all cap surface faces to volume mesh faces: {bad_idx.size} mismatches. "
475-
f"Example face node triples (GlobalNodeID) that failed exact match: {examples.tolist()}"
476-
)
477-
cap_surface.cell_data["GlobalElementID"] = indices // 4
478-
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)
479491
boundaries = cap_surface.extract_feature_edges(boundary_edges=True, manifold_edges=False,
480492
feature_edges=False, non_manifold_edges=False)
481493
boundaries = boundaries.split_bodies()
@@ -501,19 +513,23 @@ def compute_circularity(loop_polydata):
501513
lumen_surface.point_data["GlobalNodeID"] = indices.astype(int)
502514
# Assign Global Element IDs
503515
lumen_faces = lumen_surface.point_data["GlobalNodeID"][lumen_surface.faces]
504-
lumen_faces = lumen_faces.reshape(-1, 4)[:, 1:]
505-
lumen_faces = numpy.sort(lumen_faces, axis=1)
506-
dists, indices = tet_face_tree.query(lumen_faces)
507-
if not numpy.all(numpy.isclose(dists, 0.0)):
508-
bad_idx = numpy.where(~numpy.isclose(dists, 0.0))[0]
509-
sample = bad_idx[:5]
510-
examples = lumen_faces[sample]
511-
raise ValueError(
512-
f"Failed to map all lumen surface faces to volume mesh faces: {bad_idx.size} mismatches. "
513-
f"Example face node triples (GlobalNodeID) that failed exact match: {examples.tolist()}"
514-
)
515-
lumen_surface.cell_data["GlobalElementID"] = indices // 4
516-
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)
517533
boundaries = lumen_surface.extract_feature_edges(boundary_edges=True, manifold_edges=False,
518534
feature_edges=False, non_manifold_edges=False)
519535
boundaries = boundaries.split_bodies()

0 commit comments

Comments
 (0)