Atomizer/atomizer-field/check_actual_values.py

"""
Check actual values from OP2 to understand what's correct
"""
from pyNastran.op2.op2 import OP2
from pyNastran.bdf.bdf import BDF
import numpy as np

print("="*60)
print("CHECKING ACTUAL FEA VALUES")
print("="*60)

# Load OP2
op2_file = 'test_case_beam/output/model.op2'
print(f"\nLoading OP2: {op2_file}")
op2 = OP2()
op2.read_op2(op2_file)

# Load BDF
bdf_file = 'test_case_beam/input/model.bdf'
print(f"Loading BDF: {bdf_file}")
bdf = BDF()
bdf.read_bdf(bdf_file)

print("\n1. UNIT SYSTEM:")
print("-"*60)
if hasattr(bdf, 'params') and 'UNITSYS' in bdf.params:
    unitsys = str(bdf.params['UNITSYS'].values[0])
    print(f"  PARAM UNITSYS: {unitsys}")
    if 'MN' in unitsys and 'MM' in unitsys:
        print("  This is MegaNewton-Millimeter system:")
        print("    - Length: mm")
        print("    - Force: MN (MegaNewton)")
        print("    - Stress: Pa (base SI)")
        print("    - Mass: tonne")
else:
    print("  No UNITSYS parameter found")

print("\n2. MATERIAL PROPERTIES:")
print("-"*60)
if hasattr(bdf, 'materials') and bdf.materials:
    mat = list(bdf.materials.values())[0]
    print(f"  Material ID: {mat.mid}")
    print(f"  Type: {mat.type}")
    if hasattr(mat, 'e') and mat.e:
        print(f"  Young's modulus E: {mat.e:.2e}")
        print(f"    -> E = {mat.e/1e9:.1f} GPa (if units are Pa)")
        print(f"    -> E = {mat.e/1e6:.1f} GPa (if units are kPa)")
        print(f"    -> E = {mat.e/1e3:.1f} GPa (if units are MPa)")
        print(f"    Steel E ~= 200 GPa, so units must be Pa")

print("\n3. APPLIED FORCES FROM BDF:")
print("-"*60)
total_force = 0
n_forces = 0
if hasattr(bdf, 'loads') and bdf.loads:
    for load_id, load_list in bdf.loads.items():
        for load in load_list:
            if hasattr(load, 'mag'):
                print(f"  Load ID {load_id}, Node {load.node}: {load.mag:.2e} (unit depends on UNITSYS)")
                total_force += abs(load.mag)
                n_forces += 1
                if n_forces >= 3:
                    break
        if n_forces >= 3:
            break
    print(f"  Total applied force (first 3): {total_force:.2e}")
    print(f"  In MN: {total_force:.2e} MN")
    print(f"  In N: {total_force*1e6:.2e} N")

print("\n4. DISPLACEMENT FROM OP2:")
print("-"*60)
if hasattr(op2, 'displacements') and op2.displacements:
    for subcase_id, disp in op2.displacements.items():
        # Get translation only (DOF 1-3)
        translations = disp.data[0, :, :3]
        max_trans = np.max(np.abs(translations))
        max_idx = np.unravel_index(np.argmax(np.abs(translations)), translations.shape)

        print(f"  Subcase {subcase_id}:")
        print(f"    Max translation: {max_trans:.6f} mm")
        print(f"    Location: node index {max_idx[0]}, DOF {max_idx[1]}")

print("\n5. STRESS FROM OP2 (RAW VALUES):")
print("-"*60)
# Try new API
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:
        pass

if stress_dict:
    for subcase_id, stress in stress_dict.items():
        # Stress data columns depend on element type
        # For CQUAD4: typically [fiber_distance, oxx, oyy, txy, angle, major, minor, von_mises]
        stress_data = stress.data[0, :, :]

        print(f"  Subcase {subcase_id}:")
        print(f"    Data shape: {stress_data.shape} (elements × stress_components)")
        print(f"    Stress components: {stress_data.shape[1]}")

        # Von Mises is usually last column
        von_mises = stress_data[:, -1]
        print(f"\n  Von Mises stress (column {stress_data.shape[1]-1}):")
        print(f"    Min: {np.min(von_mises):.2e}")
        print(f"    Max: {np.max(von_mises):.2e}")
        print(f"    Mean: {np.mean(von_mises):.2e}")
        print(f"    Median: {np.median(von_mises):.2e}")

        print(f"\n  Principal stresses (columns 5-6 typically):")
        if stress_data.shape[1] >= 7:
            major = stress_data[:, -3]
            minor = stress_data[:, -2]
            print(f"    Major max: {np.max(major):.2e}")
            print(f"    Minor min: {np.min(minor):.2e}")

        print(f"\n  Direct stresses (columns 1-2 typically):")
        if stress_data.shape[1] >= 3:
            sxx = stress_data[:, 1]
            syy = stress_data[:, 2]
            print(f"    sigmaxx range: {np.min(sxx):.2e} to {np.max(sxx):.2e}")
            print(f"    sigmayy range: {np.min(syy):.2e} to {np.max(syy):.2e}")

print("\n6. REACTIONS FROM OP2:")
print("-"*60)
if hasattr(op2, 'grid_point_forces') and op2.grid_point_forces:
    for subcase_id, gpf in op2.grid_point_forces.items():
        forces = gpf.data[0]
        print(f"  Subcase {subcase_id}:")
        print(f"    Data shape: {forces.shape}")
        print(f"    Max reaction (all DOF): {np.max(np.abs(forces)):.2e}")

        # Get force components (first 3 columns usually Fx, Fy, Fz)
        if forces.shape[1] >= 3:
            fx = forces[:, 0]
            fy = forces[:, 1]
            fz = forces[:, 2]
            print(f"    Fx range: {np.min(fx):.2e} to {np.max(fx):.2e}")
            print(f"    Fy range: {np.min(fy):.2e} to {np.max(fy):.2e}")
            print(f"    Fz range: {np.min(fz):.2e} to {np.max(fz):.2e}")

print("\n7. YOUR STATED VALUES:")
print("-"*60)
print("  You said:")
print("    - Stress around 117 MPa")
print("    - Force around 152,200 N")
print("\n  From OP2 raw data above:")
print("    - If Von Mises max = 1.17e+05, this is:")
print("      -> 117,000 Pa = 117 kPa = 0.117 MPa (if UNITSYS=MN-MM)")
print("      -> OR 117,000 MPa (if somehow in MPa already)")
print("\n  For force 152,200 N:")
print("    - If reactions max = 1.52e+08 from OP2:")
print("      -> 152,000,000 Pa or N·mm⁻² (if MN-MM system)")
print("      -> 152 MN = 152,000,000 N (conversion)")
print("      -> OR your expected value is 0.1522 MN = 152,200 N")

print("\n8. DIRECTIONS AND TENSORS:")
print("-"*60)
print("  Stress tensor (symmetric 3×3):")
print("    sigma = [sigmaxx  tauxy  tauxz]")
print("        [tauxy  sigmayy  tauyz]")
print("        [tauxz  tauyz  sigmazz]")
print("\n  Stored in OP2 for shells (CQUAD4) as:")
print("    [fiber_dist, sigmaxx, sigmayy, tauxy, angle, sigma_major, sigma_minor, von_mises]")
print("\n  Displacement vector (6 DOF per node):")
print("    [Tx, Ty, Tz, Rx, Ry, Rz]")
print("\n  Force/Reaction vector (6 DOF per node):")
print("    [Fx, Fy, Fz, Mx, My, Mz]")

print("\n" + "="*60)