OP2 Overview
Introduction
This is meant as a tutorial on how to use the pyNastran pyNastran.op2.op2.OP2
class
This page runs through examples relating to the OP2. The OP2 is preferred as it is much faster and easier to parse. How fast? You can read a 2 GB OP2 in 4 seconds, an 8 GB file in 15 seconds, and a 60 GB file in 1-2 minutes.
Note that a static model is a SOL 101 or SOL 144. A dynamic/”transient” solution is any transient/modal/load step/frequency based solution (e.g. 103, 109, 145).
The head/tail/file_slice methods can be found at:
These examples can be found at:
Example 1: Read Write
This example will demonstate:
reading the OP2
getting some basic information
writing the F06
our model
>>> import pyNastran
>>> pkg_path = pyNastran.__path__[0]
>>> test_path = os.path.join(pkg_path, '..', 'models', 'solid_bending')
>>> op2_filename = os.path.join(test_path, 'solid_bending.op2')
>>> f06_filename = os.path.join(test_path, 'solid_bending_out.f06')
instantiate the model
>>> from pyNastran.op2.op2 import OP2
>>> model = OP2()
>>> model.read_op2(op2_filename)
>>> print(model.get_op2_stats())
op2.displacements[1]
type=RealDisplacementArray nnodes=72
data: [t1, t2, t3, r1, r2, r3] shape=[1, 72, 6] dtype=float32
gridTypes
lsdvmns = [1]
op2.spc_forces[1]
type=RealSPCForcesArray nnodes=72
data: [t1, t2, t3, r1, r2, r3] shape=[1, 72, 6] dtype=float32
gridTypes
lsdvmns = [1]
op2.ctetra_stress[1]
type=RealSolidStressArray nelements=186 nnodes=930
nodes_per_element=5 (including centroid)
eType, cid
data: [1, nnodes, 10] where 10=[oxx, oyy, ozz, txy, tyz, txz, o1, o2, o3, von_mises]
data.shape = (1, 930, 10)
element types: CTETRA
lsdvmns = [1]
>>> model.write_f06(f06_filename)
F06:
RealDisplacementArray SUBCASE=1
RealSPCForcesArray SUBCASE=1
RealSolidStressArray SUBCASE=1 - CTETRA
>>> print(tail(f06_filename, 21))
0 186 0GRID CS 4 GP
0 CENTER X 9.658666E+02 XY -2.978357E+01 A 2.559537E+04 LX-0.02 0.20 0.98 -1.094517E+04 2.288671E+04
Y 7.329372E+03 YZ 5.895411E+02 B -7.168877E+01 LY-1.00-0.03-0.01
Z 2.454026E+04 ZX -5.050599E+03 C 7.311813E+03 LZ 0.03-0.98 0.20
0 8 X 9.658666E+02 XY -2.978357E+01 A 2.559537E+04 LX-0.02 0.20 0.98 -1.094517E+04 2.288671E+04
Y 7.329372E+03 YZ 5.895411E+02 B -7.168877E+01 LY-1.00-0.03-0.01
Z 2.454026E+04 ZX -5.050599E+03 C 7.311813E+03 LZ 0.03-0.98 0.20
0 62 X 9.658666E+02 XY -2.978357E+01 A 2.559537E+04 LX-0.02 0.20 0.98 -1.094517E+04 2.288671E+04
Y 7.329372E+03 YZ 5.895411E+02 B -7.168877E+01 LY-1.00-0.03-0.01
Z 2.454026E+04 ZX -5.050599E+03 C 7.311813E+03 LZ 0.03-0.98 0.20
0 4 X 9.658666E+02 XY -2.978357E+01 A 2.559537E+04 LX-0.02 0.20 0.98 -1.094517E+04 2.288671E+04
Y 7.329372E+03 YZ 5.895411E+02 B -7.168877E+01 LY-1.00-0.03-0.01
Z 2.454026E+04 ZX -5.050599E+03 C 7.311813E+03 LZ 0.03-0.98 0.20
0 58 X 9.658666E+02 XY -2.978357E+01 A 2.559537E+04 LX-0.02 0.20 0.98 -1.094517E+04 2.288671E+04
Y 7.329372E+03 YZ 5.895411E+02 B -7.168877E+01 LY-1.00-0.03-0.01
Z 2.454026E+04 ZX -5.050599E+03 C 7.311813E+03 LZ 0.03-0.98 0.20
1 MSC.NASTRAN JOB CREATED ON 28-JAN-12 AT 12:52:32 JANUARY 28, 2012 pyNastran v0.7.1 PAGE 3
1 * * * END OF JOB * * *
Example 2: Displacement (static)
This example will demonstate:
calculating total deflection of the nodes for a static case for an OP2
calculate von mises stress and max shear
our model
>>> import pyNastran
>>> pkg_path = pyNastran.__path__[0]
>>> test_path = os.path.join(pkg_path, '..', 'models', 'solid_bending')
>>> op2_filename = os.path.join(test_path, 'solid_bending.op2')
>>> out_filename = os.path.join(test_path, 'solid_bending.out')
instantiate the model
>>> from pyNastran.op2.op2 import OP2
>>> model = OP2()
>>> model.read_op2(op2_filename)
>>> print(model.get_op2_stats())
we’re analyzing a static problem, so itime=0
we’re also assuming subcase 1
>>> itime = 0
>>> isubcase = 1
get the displacement object
>>> disp = model.displacements[isubcase]
displacement is an array
# data = [tx, ty, tz, rx, ry, rz]
# for some itime
# all the nodes -> :
# get [tx, ty, tz] -> :3
>>> txyz = disp.data[itime, :, :3]
calculate the total deflection of the vector
>>> from numpy.linalg import norm
>>> total_xyz = norm(txyz, axis=1)
since norm’s axis parameter can be tricky, we’ll double check the length
>>> nnodes = disp.data.shape[1]
>>> assert len(total_xyz) == nnodes
we could also have found nnodes by using the attribute.
It has an underscore because the object is also used for elements. The underscore name is the SORT1 name, but your data may be in SORT2 format.
>>> nnodes2 = disp._nnodes
>>> assert nnodes == nnodes2
>>> assert nnodes == 72
Additionally, we know we have 72 nodes from the shape:
op2.displacements[1]
type=RealDisplacementArray nnodes=72
data: [t1, t2, t3, r1, r2, r3] shape=[1, 72, 6] dtype=float32
gridTypes
lsdvmns = [1]
now we’ll loop over the nodes and print the total deflection
>>> msg = 'nid, gridtype, tx, ty, tz, txyz'
>>> print(msg)
>>> for (nid, grid_type), txyz, total_xyzi in zip(disp.node_gridtype, txyz, total_xyz):
>>> msg = '%s, %s, %s, %s, %s, %s' % (nid, grid_type, txyz[0], txyz[1], txyz[2], total_xyzi)
>>> print(msg)
nid, gridtype, tx, ty, tz, txyz
1, 1, 0.00764469, 4.01389e-05, 0.000111137, 0.00764561
2, 1, 0.00762899, 5.29171e-05, 0.000142154, 0.0076305
3, 1, 0.00944763, 6.38675e-05, 7.66179e-05, 0.00944816
4, 1, 0.00427092, 2.62277e-05, 7.27848e-05, 0.00427162
5, 1, 0.00152884, 1.71054e-05, -3.47525e-06, 0.00152894
...
Example 3: Eigenvector (transient)
NA
Example 4: Solid Stress (static)
This example will demonstate:
calculate von mises stress and max shear for solid elements for a static case for an OP2
our model
>>> import pyNastran
>>> pkg_path = pyNastran.__path__[0]
>>> test_path = os.path.join(pkg_path, '..', 'models', 'solid_bending')
>>> op2_filename = os.path.join(test_path, 'solid_bending.op2')
>>> out_filename = os.path.join(test_path, 'solid_bending.out')
instantiate the model
>>> from pyNastran.op2.op2 import OP2
>>> model = OP2()
>>> model.read_op2(op2_filename)
>>> print(model.get_op2_stats())
op2.ctetra_stress[1]
type=RealSolidStressArray nelements=186 nnodes=930
nodes_per_element=5 (including centroid)
eType, cid
data: [1, nnodes, 10] where 10=[oxx, oyy, ozz, txy, tyz, txz, o1, o2, o3, von_mises]
data.shape = (1, 930, 10)
element types: CTETRA
lsdvmns = [1]
we’re analyzing a static problem, so itime=0
we’re also assuming subcase 1
>>> itime = 0
>>> isubcase = 1
get the stress object (there is also cpenta_stress and chexa_stress as well as ctetra_strain/cpenta_strain/chexa_strain)
>>> stress = model.ctetra_stress[isubcase]
The stress/strain data can often be von_mises/max_shear (same for fiber_distance/curvature), so check!
#data = [oxx, oyy, ozz, txy, tyz, txz, o1, o2, o3, von_mises]
>>> o1 = stress.data[itime, :, 6]
>>> o3 = stress.data[itime, :, 8]
>>> if stress.is_von_mises():
>>> max_shear = (o1 - o3) / 2.
>>> von_mises = stress.data[itime, :, 9]
>>> else:
>>> from numpy import sqrt
>>> o2 = data[itime, :, 8]
>>> von_mises = sqrt(0.5*((o1-o2)**2 + (o2-o3)**2, (o3-o1)**2))
>>> max_shear = stress.data[itime, :, 9]
>>> for (eid, node), vm, ms in zip(stress.element_node, von_mises, max_shear):
>>> print(eid, 'CEN/4' if node == 0 else node, vm, ms)
1 CEN/4 15900.2 2957.35
1 8 15900.2 2957.35
1 13 15900.2 2957.35
1 67 15900.2 2957.35
1 33 15900.2 2957.35
2 CEN/4 16272.3 6326.18
2 8 16272.3 6326.18
2 7 16272.3 6326.18
2 62 16272.3 6326.18
2 59 16272.3 6326.18
Note that because element_node is an integer array, the centroid is 0. We renamed it to CEN/4 when we wrote it
Example 4: Solid Stress (static)
This example will demonstate:
calculating total deflection of the nodes for a dynamic case for an OP2
our model
>>> import pyNastran
>>> pkg_path = pyNastran.__path__[0]
>>> test_path = os.path.join(pkg_path, '..', 'models', 'plate_py')
>>> op2_filename = os.path.join(test_path, 'plate_py.op2')
ut_filename = os.path.join(test_path, ‘solid_bending.out’)
instantiate the model
>>> from pyNastran.op2.op2 import OP2
>>> model = OP2()
>>> model.read_op2(op2_filename)
>>> print(model.get_op2_stats())
op2.eigenvectors[1]
type=RealEigenvectorArray ntimes=10 nnodes=231
data: [t1, t2, t3, r1, r2, r3] shape=[10, 231, 6] dtype=float32
gridTypes
modes = [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]
eigrs = [-0.00037413835525512695, -0.00022113323211669922, -0.0001882314682006836, -0.00010025501251220703, 0.0001621246337890625, 0.00
07478296756744385, 1583362560.0, 2217974016.0, 10409966592.0, 11627085824.0]
mode_cycles = [0, 0, 0, 0, 0, 0, 0, 0, 0, 0]
>>> isubcase = 1
>>> eigenvector = model.eigenvectors[isubcase]
“time/mode/frequency are stored by id, so to get mode 5:
>>> modes = eigenvector._times # it may not be "time" so we don't use the name "time"
>>> from numpy import where
>>> imode5 = where(modes == 5)[0]
>>> txyz = eigenvector.data[imode5, :, :3]
calculate the total deflection of the vector
>>> from numpy.linalg import norm
>>> total_xyz = norm(txyz, axis=1)
get the eigenvalue
>>> print('eigr5 = %s' % eigenvector.eigrs[imode5])
eigr5 = 0.000162124633789
Example 5: Isotropic Plate Stress (static)
This example will demonstate:
print the fiber distance and the max principal stress for a static case for an OP2
our model
>>> import pyNastran
>>> pkg_path = pyNastran.__path__[0]
>>> test_path = os.path.join(pkg_path, '..', 'models', 'sol_101_elements')
>>> op2_filename = os.path.join(test_path, 'static_solid_shell_bar.op2')
instantiate the model
>>> from pyNastran.op2.op2 import OP2
>>> model = OP2()
>>> model.read_op2(op2_filename)
>>> print(model.get_op2_stats())
op2.cquad4_stress[1]
type=RealPlateStressArray nelements=2 nnodes_per_element=5 nlayers=2 ntotal=20
data: [1, ntotal, 8] where 8=[fiber_distance, oxx, oyy, txy, angle, omax, omin, von_mises]
data.shape=(1L, 20L, 8L)
element types: CQUAD4
lsdvmns = [1]
>>> isubcase = 1
>>> itime = 0 # this is a static case
>>> stress = model.cquad4_stress[isubcase]
>>> assert stress.nnodes == 5, 'this is a bilinear quad'
write the data
#[fiber_dist, oxx, oyy, txy, angle, majorP, minorP, ovm]
>>> eids = stress.element_node[:, 0]
>>> nids = stress.element_node[:, 1]
>>> if stress.is_fiber_distance():
>>> fiber_dist = stress.data[itime, :, 0]
>>> else:
>>> raise RuntimeError('found fiber curvature; expected fiber distance')
>>> maxp = stress.data[itime, :, 5]
>>> for (eid, nid, fdi, maxpi) in zip(eids, nids, fiber_dist, maxp):
>>> print(eid, 'CEN/4' if nid == 0 else nid, fdi, maxpi)
6 CEN/4 -0.125 8022.26
6 CEN/4 0.125 12015.9
6 4 -0.125 7580.84
6 4 0.125 11872.9
6 1 -0.125 8463.42
6 1 0.125 12158.9
6 14 -0.125 8463.69
6 14 0.125 12158.9
6 15 -0.125 7581.17
6 15 0.125 11872.9
7 CEN/4 -0.125 10016.3
7 CEN/4 0.125 10019.5
7 3 -0.125 10307.1
7 3 0.125 10311.0
7 2 -0.125 9725.54
7 2 0.125 9727.9
7 17 -0.125 9725.54
7 17 0.125 9728.06
7 16 -0.125 10307.1
7 16 0.125 10311.1
note we have 2 layers (upper and lower surface) for any PSHELL-based elements
Example 6: Composite Plate Stress (static)
This example will demonstate:
print the fiber distance and the max principal stress for a static case for an OP2
our model
>>> import pyNastran
>>> pkg_path = pyNastran.__path__[0]
>>> test_path = os.path.join(pkg_path, '..', 'models', 'sol_101_elements')
>>> op2_filename = os.path.join(test_path, 'static_solid_shell_bar.op2')
instantiate the model
>>> from pyNastran.op2.op2 import OP2
>>> model = OP2()
>>> model.read_op2(op2_filename)
>>> print(model.get_op2_stats())
op2.ctria3_composite_stress[1]
type=RealCompositePlateStressArray nelements=4 ntotal=18
data: [1, ntotal, 9] where 9=[o11, o22, t12, t1z, t2z, angle, major, minor, max_shear]
data.shape = (1, 18, 9)
element types: CTRIA3
lsdvmns = [1]
>>> isubcase = 1
>>> itime = 0 # this is a static case
>>> stress = model.ctria3_composite_stress[isubcase]
In the previous example, we had an option for a variable number of nodes for the CQUAD4s (1/5), but only nnodes=1 for the CTRIA3s.
In this example, we have 4 layers on one element and 5 on another, but they’re all at the centroid.
#[o11, o22, t12, t1z, t2z, angle, major, minor, ovm]
>>> eids = stress.element_layer[:, 0]
>>> layers = stress.element_layer[:, 1]
>>> maxp = stress.data[itime, :, 6]
>>> if stress.is_fiber_distance():
>>> fiber_dist = stress.data[itime, :, 0]
>>> else:
>>> raise RuntimeError('found fiber curvature; expected fiber distance')
>>> maxp = stress.data[itime, :, 5]
>>> for (eid, layer, maxpi) in zip(eids, layers, maxp):
>>> print(eid, 'CEN/4', layer, maxpi)
7 CEN/4 1 89.3406
7 CEN/4 2 89.3745
7 CEN/4 3 89.4313
7 CEN/4 4 89.5115
8 CEN/4 1 -85.6691
8 CEN/4 2 -85.6121
8 CEN/4 3 -85.5193
8 CEN/4 4 -85.3937
8 CEN/4 5 -85.2394
9 CEN/4 1 86.3663
9 CEN/4 2 86.6389
9 CEN/4 3 87.0977
9 CEN/4 4 87.7489
10 CEN/4 1 -87.6962
10 CEN/4 2 -87.4949
10 CEN/4 3 -87.1543
10 CEN/4 4 -86.6662
10 CEN/4 5 -86.0192