Thank you, you were absolutely right, and it was user (my) error in creating my mesh. I used the following code in readers/base.py to pinpoint the problematic faces and fixed my bug.
def _dump_unpaired_faces(self, resid, prefix="debug_unpaired_faces"):
"""
Write two artifacts in the current working directory when unpaired faces
are detected:
- {prefix}.csv : one row per unpaired face with element/face ids, nodes,
centroid and bbox
- {prefix}.vtk : legacy VTK POLYDATA with the faces as lines/tri/quads
"""
import os, csv
from collections import OrderedDict
csv_path = f"{prefix}.csv"
vtk_path = f"{prefix}.vtk"
# Build rows for CSV and collect unique point ids for VTK
rows = []
unique_node_ids = OrderedDict() # preserves insertion order
idx_of = {} # global node id -> local sequential id for VTK
def _add_point(pid):
if pid not in unique_node_ids:
unique_node_ids[pid] = len(unique_node_ids)
return unique_node_ids[pid]
# Compose rows
for nodes_tuple, fdesc in resid.items():
petype, eidx, fidx, flags = fdesc # ('tri'|'quad'|...), element index, face index, flags
node_ids = list(nodes_tuple)
pts = self._nodepts[node_ids, :]
centroid = pts.mean(axis=0)
xmin, ymin, zmin = pts.min(axis=0)
xmax, ymax, zmax = pts.max(axis=0)
rows.append({
"petype": petype,
"elem_index": int(eidx),
"face_index": int(fidx),
"flags": int(flags),
"nnodes": len(node_ids),
"node_ids": " ".join(str(int(n)) for n in node_ids),
"centroid_x": float(centroid[0]),
"centroid_y": float(centroid[1]),
"centroid_z": float(centroid[2]) if pts.shape[1] == 3 else 0.0,
"bbox_xmin": float(xmin),
"bbox_ymin": float(ymin),
"bbox_zmin": float(zmin) if pts.shape[1] == 3 else 0.0,
"bbox_xmax": float(xmax),
"bbox_ymax": float(ymax),
"bbox_zmax": float(zmax) if pts.shape[1] == 3 else 0.0,
})
for gid in node_ids:
_add_point(int(gid))
# Write CSV
fieldnames = [
"petype","elem_index","face_index","flags","nnodes","node_ids",
"centroid_x","centroid_y","centroid_z",
"bbox_xmin","bbox_ymin","bbox_zmin","bbox_xmax","bbox_ymax","bbox_zmax",
]
with open(csv_path, "w", newline="") as cf:
wr = csv.DictWriter(cf, fieldnames=fieldnames)
wr.writeheader()
wr.writerows(rows)
# Prepare VTK data
# Map global node ids -> local ids
for gid, lid in unique_node_ids.items():
idx_of[int(gid)] = int(lid)
# Points array
pts_arr = self._nodepts[list(unique_node_ids.keys()), :]
# Ensure 3 components for legacy VTK
if pts_arr.shape[1] == 2:
import numpy as np
pts_arr = np.c_[pts_arr, np.zeros((pts_arr.shape[0], 1))]
# Split faces by arity
line_faces = []
tri_faces = []
quad_faces = []
other_faces = []
for nodes_tuple, _ in resid.items():
lids = [idx_of[int(n)] for n in nodes_tuple]
if len(lids) == 2:
line_faces.append(lids)
elif len(lids) == 3:
tri_faces.append(lids)
elif len(lids) == 4:
quad_faces.append(lids)
else:
other_faces.append(lids)
# Write legacy VTK POLYDATA
with open(vtk_path, "w") as vf:
vf.write("# vtk DataFile Version 3.0\n")
vf.write("Unpaired faces\n")
vf.write("ASCII\n")
vf.write("DATASET POLYDATA\n")
vf.write(f"POINTS {pts_arr.shape[0]} float\n")
for p in pts_arr:
vf.write(f"{p[0]} {p[1]} {p[2]}\n")
# LINES block
if line_faces:
total_ints = sum(1 + len(f) for f in line_faces)
vf.write(f"LINES {len(line_faces)} {total_ints}\n")
for f in line_faces:
vf.write(f"{len(f)} " + " ".join(str(i) for i in f) + "\n")
# POLYGONS block for tri and quad
n_polys = len(tri_faces) + len(quad_faces)
if n_polys:
total_ints = sum(1 + len(f) for f in tri_faces) + sum(1 + len(f) for f in quad_faces)
vf.write(f"POLYGONS {n_polys} {total_ints}\n")
for f in tri_faces + quad_faces:
vf.write(f"{len(f)} " + " ".join(str(i) for i in f) + "\n")
# If any other arity slipped in, write them as vertices so you still see them
if other_faces:
vf.write(f"VERTICES {len(other_faces)} {len(other_faces)*2}\n")
for f in other_faces:
# take first node as a marker
vf.write(f"1 {f[0]}\n")
# Helpful console summary
print(f"[PyFR] Unpaired faces: {len(rows)}")
print(f"[PyFR] Wrote centroid list to {os.path.abspath(csv_path)}")
print(f"[PyFR] Wrote face geometry to {os.path.abspath(vtk_path)}")
if any(resid.values()):
try:
self._dump_unpaired_faces(resid, prefix="debug_unpaired_faces")
except Exception as e:
print(f"[PyFR] Failed to dump unpaired faces: {e}")
raise ValueError('Unpaired faces in mesh. See debug_unpaired_faces.csv and debug_unpaired_faces.vtk')
