Skip to content

Commit 7cef22e

Browse files
pjuangphpjuangph
committed
added verify periodicity
1 parent a8f96be commit 7cef22e

File tree

6 files changed

+346
-6
lines changed

6 files changed

+346
-6
lines changed

‎python/plot3d/__init__.py‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -5,13 +5,13 @@
55
from .block import Block
66
from .blockfunctions import rotate_block, get_outer_bounds, block_connection_matrix,split_blocks, plot_blocks, reduce_blocks, find_matching_faces
77
from .block_merging_mixed_facepairs import combine_nxnxn_cubes_mixed_pairs
8-
from .connectivity import find_matching_blocks, get_face_intersection, connectivity_fast, face_matches_to_dict
8+
from .connectivity import find_matching_blocks, get_face_intersection, connectivity_fast, face_matches_to_dict, verify_connectivity
99
from .face import Face
1010
from .facefunctions import create_face_from_diagonals, get_outer_faces, find_bounding_faces,split_face,find_face_nearest_point,match_faces_dict_to_list,outer_face_dict_to_list,find_closest_block
1111
from .read import read_plot3D, read_ap_nasa
1212
from .write import write_plot3D
1313
from .differencing import find_edges, find_face_edges
14-
from .periodicity import periodicity, periodicity_fast, create_rotation_matrix, rotated_periodicity, translational_periodicity
14+
from .periodicity import periodicity, periodicity_fast, create_rotation_matrix, rotated_periodicity, translational_periodicity, verify_periodicity
1515
from .point_match import point_match
1616
from .split_block import split_blocks, Direction
1717
from .listfunctions import unique_pairs

‎python/plot3d/connectivity.py‎

Lines changed: 156 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -565,3 +565,159 @@ def face_matches_to_dict(face1:Face, face2:Face,block1:Block,block2:Block):
565565
match['block2']['JMAX'] = int(df.iloc[0]['J'])
566566
match['block2']['KMAX'] = int(df.iloc[0]['K'])
567567
return match
568+
569+
570+
def verify_connectivity(blocks: List[Block], face_matches: list, tol: float = 1E-6):
571+
"""Verifies that the diagonal corners of face_matches are spatially consistent.
572+
573+
For each face_match, checks that block1's lower corner coordinates match
574+
block2's lower corner coordinates (and similarly for upper corners) within
575+
the specified tolerance. If the stored diagonal doesn't match, tries all
576+
permutations of block2's face corners. If a valid permutation is found,
577+
the face_match is corrected and added to the verified list.
578+
579+
Uses GCD reduction (same as connectivity_fast) for efficient coordinate lookups.
580+
581+
Args:
582+
blocks (List[Block]): List of all blocks (original full-resolution)
583+
face_matches (list): List of face_match dicts from connectivity or periodicity
584+
tol (float, optional): Euclidean distance tolerance. Defaults to 1E-6.
585+
586+
Returns:
587+
(list): verified face_matches whose diagonals are confirmed or corrected
588+
(list): mismatched face_matches where no corner permutation matched
589+
"""
590+
# Compute GCD and reduce blocks (same pattern as connectivity_fast)
591+
gcd_array = list()
592+
for block_indx in range(len(blocks)):
593+
block = blocks[block_indx]
594+
gcd_array.append(math.gcd(block.IMAX-1, math.gcd(block.JMAX-1, block.KMAX-1)))
595+
gcd_to_use = min(gcd_array)
596+
reduced_blocks = reduce_blocks(deepcopy(blocks), gcd_to_use)
597+
598+
# Scale down face_matches indices by GCD
599+
scaled_matches = deepcopy(face_matches)
600+
for fm in scaled_matches:
601+
for side in ['block1', 'block2']:
602+
for key in ['IMIN', 'JMIN', 'KMIN', 'IMAX', 'JMAX', 'KMAX']:
603+
fm[side][key] = fm[side][key] // gcd_to_use
604+
605+
verified = list()
606+
mismatched = list()
607+
608+
for idx in range(len(scaled_matches)):
609+
fm = scaled_matches[idx]
610+
b1 = fm['block1']
611+
b2 = fm['block2']
612+
b1_idx = b1['block_index']
613+
b2_idx = b2['block_index']
614+
block1 = reduced_blocks[b1_idx]
615+
block2 = reduced_blocks[b2_idx]
616+
617+
# Block1 diagonal coordinates
618+
x1_l = block1.X[b1['IMIN'], b1['JMIN'], b1['KMIN']]
619+
y1_l = block1.Y[b1['IMIN'], b1['JMIN'], b1['KMIN']]
620+
z1_l = block1.Z[b1['IMIN'], b1['JMIN'], b1['KMIN']]
621+
622+
x1_u = block1.X[b1['IMAX'], b1['JMAX'], b1['KMAX']]
623+
y1_u = block1.Y[b1['IMAX'], b1['JMAX'], b1['KMAX']]
624+
z1_u = block1.Z[b1['IMAX'], b1['JMAX'], b1['KMAX']]
625+
626+
# Enumerate unique corners of block2's face
627+
I2 = [b2['IMIN'], b2['IMAX']]
628+
J2 = [b2['JMIN'], b2['JMAX']]
629+
K2 = [b2['KMIN'], b2['KMAX']]
630+
631+
unique_corners = list()
632+
seen = set()
633+
for i in I2:
634+
for j in J2:
635+
for k in K2:
636+
key = (i, j, k)
637+
if key not in seen:
638+
seen.add(key)
639+
unique_corners.append(key)
640+
641+
# Check stored diagonal first
642+
x2_l = block2.X[b2['IMIN'], b2['JMIN'], b2['KMIN']]
643+
y2_l = block2.Y[b2['IMIN'], b2['JMIN'], b2['KMIN']]
644+
z2_l = block2.Z[b2['IMIN'], b2['JMIN'], b2['KMIN']]
645+
646+
x2_u = block2.X[b2['IMAX'], b2['JMAX'], b2['KMAX']]
647+
y2_u = block2.Y[b2['IMAX'], b2['JMAX'], b2['KMAX']]
648+
z2_u = block2.Z[b2['IMAX'], b2['JMAX'], b2['KMAX']]
649+
650+
dx = x2_l - x1_l; dy = y2_l - y1_l; dz = z2_l - z1_l
651+
d_lower = math.sqrt(dx*dx + dy*dy + dz*dz)
652+
dx = x2_u - x1_u; dy = y2_u - y1_u; dz = z2_u - z1_u
653+
d_upper = math.sqrt(dx*dx + dy*dy + dz*dz)
654+
655+
if d_lower < tol and d_upper < tol:
656+
verified.append(face_matches[idx])
657+
continue
658+
659+
# Try all permutations of block2's corners
660+
found = False
661+
best_d_lower = d_lower
662+
best_d_upper = d_upper
663+
664+
for corner_lower in unique_corners:
665+
for corner_upper in unique_corners:
666+
if corner_lower == corner_upper:
667+
continue
668+
669+
il, jl, kl = corner_lower
670+
iu, ju, ku = corner_upper
671+
672+
x2_l = block2.X[il, jl, kl]
673+
y2_l = block2.Y[il, jl, kl]
674+
z2_l = block2.Z[il, jl, kl]
675+
676+
x2_u = block2.X[iu, ju, ku]
677+
y2_u = block2.Y[iu, ju, ku]
678+
z2_u = block2.Z[iu, ju, ku]
679+
680+
dx = x2_l - x1_l; dy = y2_l - y1_l; dz = z2_l - z1_l
681+
dl = math.sqrt(dx*dx + dy*dy + dz*dz)
682+
dx = x2_u - x1_u; dy = y2_u - y1_u; dz = z2_u - z1_u
683+
du = math.sqrt(dx*dx + dy*dy + dz*dz)
684+
685+
if dl < best_d_lower:
686+
best_d_lower = dl
687+
if du < best_d_upper:
688+
best_d_upper = du
689+
690+
if dl < tol and du < tol:
691+
corrected = deepcopy(face_matches[idx])
692+
corrected['block2']['IMIN'] = il * gcd_to_use
693+
corrected['block2']['JMIN'] = jl * gcd_to_use
694+
corrected['block2']['KMIN'] = kl * gcd_to_use
695+
corrected['block2']['IMAX'] = iu * gcd_to_use
696+
corrected['block2']['JMAX'] = ju * gcd_to_use
697+
corrected['block2']['KMAX'] = ku * gcd_to_use
698+
verified.append(corrected)
699+
if b1_idx == b2_idx:
700+
print("verify_connectivity: Self-match corrected for block index {0}".format(b1_idx))
701+
found = True
702+
break
703+
if found:
704+
break
705+
706+
if not found:
707+
orig = face_matches[idx]
708+
b1_orig = orig['block1']
709+
b2_orig = orig['block2']
710+
print(f"verify_connectivity: MISMATCH at face_match index {idx}")
711+
print(f" block1 (block_index={b1_orig['block_index']}): "
712+
f"lower=({b1_orig['IMIN']},{b1_orig['JMIN']},{b1_orig['KMIN']}) "
713+
f"upper=({b1_orig['IMAX']},{b1_orig['JMAX']},{b1_orig['KMAX']})")
714+
print(f" block2 (block_index={b2_orig['block_index']}): "
715+
f"lower=({b2_orig['IMIN']},{b2_orig['JMIN']},{b2_orig['KMIN']}) "
716+
f"upper=({b2_orig['IMAX']},{b2_orig['JMAX']},{b2_orig['KMAX']})")
717+
print(f" block1 lower xyz = ({x1_l:.6e}, {y1_l:.6e}, {z1_l:.6e})")
718+
print(f" block1 upper xyz = ({x1_u:.6e}, {y1_u:.6e}, {z1_u:.6e})")
719+
print(f" Closest block2 corner dist to block1 lower: {best_d_lower:.6e}")
720+
print(f" Closest block2 corner dist to block1 upper: {best_d_upper:.6e}")
721+
mismatched.append(face_matches[idx])
722+
723+
return verified, mismatched

‎python/plot3d/gridpro/import_functions.py‎

Lines changed: 12 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,8 +1,7 @@
1-
from typing import Dict, List, Tuple, Optional, Union, Set
1+
from typing import Dict, List, Tuple, Optional, Union, Set, Any
22
import numpy as np
33
import pandas as pd
44
from tqdm import tqdm
5-
from typing import List, Tuple, Optional, Any
65
from plot3d import Block
76
from collections import defaultdict
87

@@ -221,6 +220,8 @@ def parse_patch(tokens: List[str]) -> None:
221220
parse_patch(toks)
222221

223222
patches_df = pd.DataFrame(patch_rows)
223+
patches_to_connectivity = patches_df[(patches_df['sb1'] > 0) & (patches_df['sb2'] > 0) & (patches_df['pty'] > 1)]
224+
224225

225226
# ---------------------------- helpers ----------------------------
226227
def face_dict(sb: int, imin: int, jmin: int, kmin: int, imax: int, jmax: int, kmax: int, pty:int=-1) -> Dict[str, int]:
@@ -253,6 +254,14 @@ def pair_dict(sb1: int, r1: Tuple[int,int,int,int,int,int],
253254
)
254255
)
255256

257+
# These are patches that should be in connectivity but aren't
258+
connectivity_to_check: List[Dict[str, Dict[str, int]]] = []
259+
for p in patches_to_connectivity.itertuples(index=False):
260+
connectivity_to_check.append(pair_dict(
261+
p.sb1, (p.L1i, p.L1j, p.L1k, p.H1i, p.H1j, p.H1k),
262+
p.sb2, (p.L2i, p.L2j, p.L2k, p.H2i, p.H2j, p.H2k),
263+
))
264+
256265
# Boundary-condition groups (by pty)
257266
def bc_faces_for(name: str, pty_vals: List[int], bc_type: str) -> List[Dict[str, int]]:
258267
targets = set(pty_vals)
@@ -384,6 +393,7 @@ def add_bc_group(name: str, ids: List[int], bc_type: str) -> None:
384393
"periodic_faces": periodic_faces,
385394
"volume_zones": volume_zones,
386395
"blocksizes": superblock_sizes,
396+
"connectivity_to_check": connectivity_to_check,
387397
"patches": patches_df,
388398
}
389399

‎python/plot3d/periodicity.py‎

Lines changed: 174 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -904,3 +904,177 @@ def __periodicity_check__(face1:Face, face2:Face,block1:Block,block2:Block,tol:f
904904
return df,periodic_faces,split_faces
905905

906906

907+
def verify_periodicity(blocks: List[Block], face_matches: list, theta: float, rotation_axis: str = 'x', tol: float = 1E-6):
908+
"""Verifies that the diagonal corners of periodic face_matches are spatially
909+
consistent after rotating block1 by ±theta about the given axis.
910+
911+
For each face_match, rotates block1 by +theta (and if needed -theta) and checks
912+
that the rotated lower/upper corner coordinates match block2's lower/upper
913+
corners within tolerance. If the stored diagonal doesn't match, tries all
914+
permutations of block2's face corners. If a valid permutation is found,
915+
the face_match is corrected and added to the verified list.
916+
917+
Uses GCD reduction (same as connectivity_fast / periodicity_fast) for
918+
efficient coordinate lookups.
919+
920+
Args:
921+
blocks (List[Block]): List of all blocks (original full-resolution)
922+
face_matches (list): List of face_match dicts from periodicity or rotated_periodicity
923+
theta (float): Rotation angle in degrees
924+
rotation_axis (str, optional): Axis of rotation 'x', 'y', or 'z'. Defaults to 'x'.
925+
tol (float, optional): Euclidean distance tolerance. Defaults to 1E-6.
926+
927+
Returns:
928+
(list): verified face_matches whose diagonals are confirmed or corrected
929+
(list): mismatched face_matches where no corner permutation matched
930+
"""
931+
# Compute GCD and reduce blocks (same pattern as connectivity_fast)
932+
gcd_array = list()
933+
for block_indx in range(len(blocks)):
934+
block = blocks[block_indx]
935+
gcd_array.append(math.gcd(block.IMAX-1, math.gcd(block.JMAX-1, block.KMAX-1)))
936+
gcd_to_use = min(gcd_array)
937+
reduced_blocks = reduce_blocks(deepcopy(blocks), gcd_to_use)
938+
939+
# Build rotation matrices for +theta and -theta
940+
rotation_matrix_pos = create_rotation_matrix(radians(theta), rotation_axis)
941+
rotation_matrix_neg = create_rotation_matrix(radians(-theta), rotation_axis)
942+
943+
# Pre-rotate all reduced blocks in both directions
944+
rotated_blocks_pos = [rotate_block(b, rotation_matrix_pos) for b in reduced_blocks]
945+
rotated_blocks_neg = [rotate_block(b, rotation_matrix_neg) for b in reduced_blocks]
946+
947+
# Scale down face_matches indices by GCD
948+
scaled_matches = deepcopy(face_matches)
949+
for fm in scaled_matches:
950+
for side in ['block1', 'block2']:
951+
for key in ['IMIN', 'JMIN', 'KMIN', 'IMAX', 'JMAX', 'KMAX']:
952+
fm[side][key] = fm[side][key] // gcd_to_use
953+
954+
verified = list()
955+
mismatched = list()
956+
957+
for idx in range(len(scaled_matches)):
958+
fm = scaled_matches[idx]
959+
b1 = fm['block1']
960+
b2 = fm['block2']
961+
b1_idx = b1['block_index']
962+
b2_idx = b2['block_index']
963+
block2 = reduced_blocks[b2_idx]
964+
965+
# Enumerate unique corners of block2's face
966+
I2 = [b2['IMIN'], b2['IMAX']]
967+
J2 = [b2['JMIN'], b2['JMAX']]
968+
K2 = [b2['KMIN'], b2['KMAX']]
969+
970+
unique_corners = list()
971+
seen = set()
972+
for i in I2:
973+
for j in J2:
974+
for k in K2:
975+
key = (i, j, k)
976+
if key not in seen:
977+
seen.add(key)
978+
unique_corners.append(key)
979+
980+
found = False
981+
best_d_lower = float('inf')
982+
best_d_upper = float('inf')
983+
984+
# Try +theta rotation first, then -theta
985+
for rotated_blocks in [rotated_blocks_pos, rotated_blocks_neg]:
986+
if found:
987+
break
988+
989+
block1_rotated = rotated_blocks[b1_idx]
990+
991+
# Block1 rotated diagonal coordinates
992+
x1_l = block1_rotated.X[b1['IMIN'], b1['JMIN'], b1['KMIN']]
993+
y1_l = block1_rotated.Y[b1['IMIN'], b1['JMIN'], b1['KMIN']]
994+
z1_l = block1_rotated.Z[b1['IMIN'], b1['JMIN'], b1['KMIN']]
995+
996+
x1_u = block1_rotated.X[b1['IMAX'], b1['JMAX'], b1['KMAX']]
997+
y1_u = block1_rotated.Y[b1['IMAX'], b1['JMAX'], b1['KMAX']]
998+
z1_u = block1_rotated.Z[b1['IMAX'], b1['JMAX'], b1['KMAX']]
999+
1000+
# Check stored diagonal first
1001+
x2_l = block2.X[b2['IMIN'], b2['JMIN'], b2['KMIN']]
1002+
y2_l = block2.Y[b2['IMIN'], b2['JMIN'], b2['KMIN']]
1003+
z2_l = block2.Z[b2['IMIN'], b2['JMIN'], b2['KMIN']]
1004+
1005+
x2_u = block2.X[b2['IMAX'], b2['JMAX'], b2['KMAX']]
1006+
y2_u = block2.Y[b2['IMAX'], b2['JMAX'], b2['KMAX']]
1007+
z2_u = block2.Z[b2['IMAX'], b2['JMAX'], b2['KMAX']]
1008+
1009+
dx = x2_l - x1_l; dy = y2_l - y1_l; dz = z2_l - z1_l
1010+
d_lower = math.sqrt(dx*dx + dy*dy + dz*dz)
1011+
dx = x2_u - x1_u; dy = y2_u - y1_u; dz = z2_u - z1_u
1012+
d_upper = math.sqrt(dx*dx + dy*dy + dz*dz)
1013+
1014+
if d_lower < best_d_lower:
1015+
best_d_lower = d_lower
1016+
if d_upper < best_d_upper:
1017+
best_d_upper = d_upper
1018+
1019+
if d_lower < tol and d_upper < tol:
1020+
verified.append(face_matches[idx])
1021+
found = True
1022+
break
1023+
1024+
# Try all permutations of block2's corners
1025+
for corner_lower in unique_corners:
1026+
for corner_upper in unique_corners:
1027+
if corner_lower == corner_upper:
1028+
continue
1029+
1030+
il, jl, kl = corner_lower
1031+
iu, ju, ku = corner_upper
1032+
1033+
x2_l = block2.X[il, jl, kl]
1034+
y2_l = block2.Y[il, jl, kl]
1035+
z2_l = block2.Z[il, jl, kl]
1036+
1037+
x2_u = block2.X[iu, ju, ku]
1038+
y2_u = block2.Y[iu, ju, ku]
1039+
z2_u = block2.Z[iu, ju, ku]
1040+
1041+
dx = x2_l - x1_l; dy = y2_l - y1_l; dz = z2_l - z1_l
1042+
dl = math.sqrt(dx*dx + dy*dy + dz*dz)
1043+
dx = x2_u - x1_u; dy = y2_u - y1_u; dz = z2_u - z1_u
1044+
du = math.sqrt(dx*dx + dy*dy + dz*dz)
1045+
1046+
if dl < best_d_lower:
1047+
best_d_lower = dl
1048+
if du < best_d_upper:
1049+
best_d_upper = du
1050+
1051+
if dl < tol and du < tol:
1052+
corrected = deepcopy(face_matches[idx])
1053+
corrected['block2']['IMIN'] = il * gcd_to_use
1054+
corrected['block2']['JMIN'] = jl * gcd_to_use
1055+
corrected['block2']['KMIN'] = kl * gcd_to_use
1056+
corrected['block2']['IMAX'] = iu * gcd_to_use
1057+
corrected['block2']['JMAX'] = ju * gcd_to_use
1058+
corrected['block2']['KMAX'] = ku * gcd_to_use
1059+
verified.append(corrected)
1060+
found = True
1061+
break
1062+
if found:
1063+
break
1064+
1065+
if not found:
1066+
orig = face_matches[idx]
1067+
b1_orig = orig['block1']
1068+
b2_orig = orig['block2']
1069+
print(f"verify_periodicity: MISMATCH at face_match index {idx}")
1070+
print(f" block1 (block_index={b1_orig['block_index']}): "
1071+
f"lower=({b1_orig['IMIN']},{b1_orig['JMIN']},{b1_orig['KMIN']}) "
1072+
f"upper=({b1_orig['IMAX']},{b1_orig['JMAX']},{b1_orig['KMAX']})")
1073+
print(f" block2 (block_index={b2_orig['block_index']}): "
1074+
f"lower=({b2_orig['IMIN']},{b2_orig['JMIN']},{b2_orig['KMIN']}) "
1075+
f"upper=({b2_orig['IMAX']},{b2_orig['JMAX']},{b2_orig['KMAX']})")
1076+
print(f" Closest rotated block1 corner dist to block2 lower: {best_d_lower:.6e}")
1077+
print(f" Closest rotated block1 corner dist to block2 upper: {best_d_upper:.6e}")
1078+
mismatched.append(face_matches[idx])
1079+
1080+
return verified, mismatched

0 commit comments

Comments
 (0)