Files
Atomizer/atomizer-field/check_op2_units.py

123 lines
4.0 KiB
Python
Raw Normal View History

"""
Check actual units from pyNastran OP2 reader
"""
from pyNastran.op2.op2 import OP2
import numpy as np
print("="*60)
print("CHECKING OP2 FILE UNITS")
print("="*60)
# Load OP2 file
op2_file = 'test_case_beam/output/model.op2'
print(f"\nLoading: {op2_file}")
op2 = OP2()
op2.read_op2(op2_file)
print("\n1. DISPLACEMENT DATA:")
print("-"*60)
if hasattr(op2, 'displacements') and op2.displacements:
for subcase_id, disp in op2.displacements.items():
print(f" Subcase {subcase_id}:")
print(f" Data shape: {disp.data.shape}")
print(f" Max displacement (all DOF): {np.max(np.abs(disp.data)):.6f}")
print(f" Translation max (DOF 1-3): {np.max(np.abs(disp.data[0, :, :3])):.6f}")
# Check if pyNastran has unit info
if hasattr(disp, 'units'):
print(f" Units: {disp.units}")
else:
print(f" Units: Not specified by pyNastran")
print("\n2. STRESS DATA:")
print("-"*60)
# Try new API first
stress_dict = None
if hasattr(op2, 'op2_results') and hasattr(op2.op2_results, 'stress'):
if hasattr(op2.op2_results.stress, 'cquad4_stress'):
stress_dict = op2.op2_results.stress.cquad4_stress
elif hasattr(op2, 'cquad4_stress'):
try:
stress_dict = op2.cquad4_stress
except:
stress_dict = None
if stress_dict:
for subcase_id, stress in stress_dict.items():
print(f" Subcase {subcase_id}:")
print(f" Data shape: {stress.data.shape}")
# Get von Mises stress (last column)
von_mises = stress.data[0, :, -1]
print(f" Von Mises min: {np.min(von_mises):.2e}")
print(f" Von Mises max: {np.max(von_mises):.2e}")
print(f" Von Mises mean: {np.mean(von_mises):.2e}")
# Check if pyNastran has unit info
if hasattr(stress, 'units'):
print(f" Units: {stress.units}")
else:
print(f" Units: Not specified by pyNastran")
# Check data type names
if hasattr(stress, 'data_names'):
print(f" Data names: {stress.data_names}")
else:
print(" No CQUAD4 stress data found")
print("\n3. REACTION FORCES:")
print("-"*60)
if hasattr(op2, 'grid_point_forces') and op2.grid_point_forces:
for subcase_id, forces in op2.grid_point_forces.items():
print(f" Subcase {subcase_id}:")
print(f" Data shape: {forces.data.shape}")
print(f" Max force: {np.max(np.abs(forces.data)):.2e}")
if hasattr(forces, 'units'):
print(f" Units: {forces.units}")
else:
print(f" Units: Not specified by pyNastran")
print("\n4. BDF FILE UNITS:")
print("-"*60)
# Try to read BDF to check units
from pyNastran.bdf.bdf import BDF
bdf_file = 'test_case_beam/input/model.bdf'
print(f"Loading: {bdf_file}")
bdf = BDF()
bdf.read_bdf(bdf_file)
# Check for PARAM cards that define units
if hasattr(bdf, 'params') and bdf.params:
print(" PARAM cards found:")
for param_name, param in bdf.params.items():
if 'UNIT' in param_name.upper() or 'WTMASS' in param_name.upper():
print(f" {param_name}: {param}")
# Check material properties (can infer units from magnitude)
if hasattr(bdf, 'materials') and bdf.materials:
print("\n Material properties (first material):")
mat = list(bdf.materials.values())[0]
print(f" Type: {mat.type}")
if hasattr(mat, 'e') and mat.e:
print(f" Young's modulus E: {mat.e:.2e}")
print(f" If E~200,000 then units are MPa (steel)")
print(f" If E~200,000,000 then units are Pa (steel)")
print(f" If E~29,000,000 then units are psi (steel)")
print("\n5. NASTRAN UNIT SYSTEM ANALYSIS:")
print("-"*60)
print(" Common Nastran unit systems:")
print(" - SI: m, kg, N, Pa, s")
print(" - mm-tonne-N: mm, tonne, N, MPa, s")
print(" - mm-kg-N: mm, kg, N, MPa, s")
print(" - in-lb-s: in, lb, lbf, psi, s")
print("")
print(" Our metadata claims: mm, kg, N, MPa")
print(" But pyNastran might return stress in Pa (base SI)")
print(" This would explain the 1000× error!")
print("\n" + "="*60)