pyNastran
pyNastran copied to clipboard
SteveDoyle2
Coordinate transformation of stresses/strains using CSTM, TRMBU, TRMBD result tables
Hello Steve,
I would like to thank you for your great work on pyNastran. It's a wonderful package.
I'm writing a python script that can convert OP2 results (stresses/strains) from the output coordinate system to other coordinate systems (basic or material). The software I use for pre-processing is Simcenter 3D.
versions in use:
- pyNastran: latest commit to main (I downloaded the zip file).
- python: 3.8
My plan was to use the CSTM, TRMBU and TRMBD result tables from the OP2 file to create the necessary transformation matrices T and get the transformed stress results through matrix multiplication: S' = T @ S @ T^T. Simcenter's documentation mentions this about the tables:
- CSTM: Coordinate system transformation matrices to transform from the global to the basic coordinate system.
- TRMBD: Contains euler angles to transform from the material to the basic coordinate system, in the 'deformed' configuration.
- TRMBU: Contains euler angles to transform from the material to the basic coordinate system, in the 'undeformed' configuration.
There's some code added to pyNastran in order to read the TRMBU/D tables. I've attached the edited files to this issue. For now, I've been able to get the stress results in the absolute coordinate system for a linear solution by using the transformation matrices of CSTM, but not for a non-linear one (the results start out close, but diverge with increasing time step). In the attached zip file, I added a small example code which performs this transformation on a given op2 file. I also attached the op2 files that I used while testing. The code is copied below:
import numpy as np
from numpy import pi, cos, sin
from pyNastran.op2.op2 import OP2
from pyNastran.op2.tables.oes_stressStrain.oes import RealSolidStressArrayNx
from pyNastran.op2.result_objects.op2_results import CSTM, TRMBU, TRMBD
def get_stresses(op2: OP2, nonlinear: bool):
# Get stress results
subcase_id = 1
chexa_stress = op2.chexa_stress[subcase_id] # type: RealSolidStressArrayNx
ntimes = chexa_stress.ntimes
nelements = chexa_stress.nelements
# Reshape results (Time, Element, Node, Stress component)
if nonlinear is True:
reshaped_results = chexa_stress.data.reshape(ntimes, nelements, 8, 7)
else:
# Remove elemental stresses
reshaped_results = np.delete(chexa_stress.data, np.arange(0, chexa_stress.data.shape[1], 9), axis=1)
reshaped_results = reshaped_results.reshape(1, nelements, 8, 10)
return reshaped_results
def test_cstm(op2: OP2, stresses):
ntimes = stresses.shape[0]
nelements = stresses.shape[1]
# Create 3x3 stress matrix
stress_mat = np.empty([ntimes, nelements, 8, 3, 3], dtype="float64")
stress_mat[:, :, :, 0, 0] = stresses[:, :, :, 0] # XX
stress_mat[:, :, :, 0, 1] = stresses[:, :, :, 3] # XY
stress_mat[:, :, :, 0, 2] = stresses[:, :, :, 5] # XZ
stress_mat[:, :, :, 1, 0] = stresses[:, :, :, 3] # XY
stress_mat[:, :, :, 1, 1] = stresses[:, :, :, 1] # YY
stress_mat[:, :, :, 1, 2] = stresses[:, :, :, 4] # YZ
stress_mat[:, :, :, 2, 0] = stresses[:, :, :, 5] # XZ
stress_mat[:, :, :, 2, 1] = stresses[:, :, :, 4] # YZ
stress_mat[:, :, :, 2, 2] = stresses[:, :, :, 2] # ZZ
# Create Global -> Basic T matrices
cstm = op2.cstm # type: CSTM
gbt = np.empty([nelements, 3, 3], dtype="float64")
gbt[:, 0, 0] = cstm.data[:, cstm.headers["T11"]]
gbt[:, 0, 1] = cstm.data[:, cstm.headers["T12"]]
gbt[:, 0, 2] = cstm.data[:, cstm.headers["T13"]]
gbt[:, 1, 0] = cstm.data[:, cstm.headers["T21"]]
gbt[:, 1, 1] = cstm.data[:, cstm.headers["T22"]]
gbt[:, 1, 2] = cstm.data[:, cstm.headers["T23"]]
gbt[:, 2, 0] = cstm.data[:, cstm.headers["T31"]]
gbt[:, 2, 1] = cstm.data[:, cstm.headers["T32"]]
gbt[:, 2, 2] = cstm.data[:, cstm.headers["T33"]]
# Create transpose of T matrices
gbt_T = np.transpose(gbt, (0, 2, 1))
# Transform stresses to basic
stress_mat_b = np.empty([ntimes, nelements, 8, 3, 3])
for time_i in range(0, ntimes):
for node_i in range(0, 8):
stress_mat_b[time_i, :, node_i, :, :] = gbt @ stress_mat[time_i, :, node_i, :, :] @ gbt_T
# Print transformed stresses for eID=1, first node (gID=1094) at time 0
print("eID = 1, nodeID=1094, time=0")
print(stress_mat_b[0, 0, 0, :, :])
# Print transformed stresses for eID=1, first node (gID=1094) at time -1
print("eID = 1, nodeID=1094, time=-1")
print(stress_mat_b[-1, 0, 0, :, :])
if __name__ == "__main__":
# Read OP2 results
op2_results = OP2(mode="nx", debug=True)
op2_results.read_op2(op2_full_path)
# Get stresses from OP2
nonlinear = True
stress = get_stresses(op2=op2_results, nonlinear=nonlinear)
# Test CSTM
test_cstm(op2_results, stress)
Results
# non-linear results Simcenter [MPa]
# XX YY ZZ XY YZ ZX
# 1.37270241e+01 5.33785641e-01 4.34055746e-01 -1.85074583e-01 -1.03073418e-01 3.19195569e-01
# non-linear results pyNastran [kPa]
# [[13727.03243347 -185.39733719 318.81154467]
# [ -185.39733719 533.79855125 -103.07372103]
# [ 318.81154467 -103.07372103 434.0334562 ]]
# non-linear results Simcenter [MPa]
# 1.36058472e+03 5.19877281e+01 4.26953926e+01 -1.84598961e+01 -9.89693642e+00 3.23777771e+01
# non-linear results pyNastran [kPa]
# [[1360679.96533602 -20678.74034289 28902.34315818]
# [ -20678.74034289 52088.9989944 -9878.25658968]
# [ 28902.34315818 -9878.25658968 42498.86770083]]
Presumably the stresses' coordinate systems change in orientation through time, but I'm not sure what extra information we need in that case to perform the correct transformations. I wrote a similar post on Reddit, and your responses prompted me to post this issue here.
I would like to thank you for taking a look at this. If you'd need any extra information or files, I'd be happy to provide them!
Best regards,
Ben Van Bavel coordinate_transformations.zip
I should mention that I also edited the cstm_reader in pyNastran. It reads the CSTM table and directly stores the extracted values in an array instead of creating CORD2X cards. That allows me to directly retrieve the transformation matrices in the code above.
The edited pyNastran file (op2_reader) is included in the zip-file of the previous comment.
The CSTM is for the undeformed state. I think you need to use the TRMBD data. As far as I can tell, you're just referencing CSTM in the deformed calculation.
I see you have a trmbu_example that calculates a transform based on the TRMBU Euler angles. Are you sure the Euler angle transform is ZYX? I would try and get the TRMBU data block (element based) to calculate the CSTM matrices (coord based). Once that's right swap out the TRMBU datablock with the TRMBD datablock.
It's probably good to start with a model with simpler loading as well.
Hi Steve,
Thanks for explaining! I'm following your advice and have been experimenting on a single CHEXA-element that is rotated through time with a simple load case. In the pre-processor I can specify a rotated material coordinate system through three subsequent euler angles.
I'm able to check that the CSTM table contains the rotation matrix produced by these three euler angles (R1 @ R2 @ R3).
import numpy as np
from numpy import pi, cos, sin
# Pre-processor euler angles: 45 deg about X, 30 deg about Y, 12 deg about Z
eul1 = 45 * pi / 180
eul2 = 30 * pi / 180
eul3 = 12 * pi / 180
c1 = cos(eul1)
s1 = sin(eul1)
c2 = cos(eul2)
s2 = sin(eul2)
c3 = cos(eul3)
s3 = sin(eul3)
unit_t_mats = {
"Z3": np.array([[c3, -s3, 0],
[s3, c3, 0],
[0, 0, 1]]),
"Y2": np.array([[c2, 0, s2],
[0, 1, 0],
[-s2, 0, c2]]),
"X1": np.array([[1, 0, 0],
[0, c1, -s1],
[0, s1, c1]])}
# Create transformation matrix global to basic
rotation_matrix = unit_t_mats["X1"] @ unit_t_mats["Y2"] @ unit_t_mats["Z3"]
print(f"CSTM Transformation matrix:\n\n{gbt}\n")
print(f"Euler angle transformation matrix:\n\n{rotation_matrix}\n")
Results
CSTM Transformation matrix:
[[ 0.84710067 -0.18005681 0.5 ]
[ 0.49284317 0.61814692 -0.61237244]
[-0.19881163 0.76516268 0.61237244]]
Euler angle transformation matrix:
[[ 0.84710067 -0.18005681 0.5 ]
[ 0.49284317 0.61814692 -0.61237244]
[-0.19881163 0.76516268 0.61237244]]
I'm also able to create the CSTM matrix from the euler angles provided in TRMBU if there's only one non-zero euler angle (example: eul1=45 deg, eul2=eul3=0). In case there are multiple non-zero euler angles, the angles that I retrieve from TRMBU are not equal to those provided in the pre-processor.
I tried a brute-force method to find the correct euler combination, to see if it was possible to create the CSTM table through some other combination:
import numpy as np
from numpy import pi, cos, sin
subcase_id = 0
eltype = "CHEXA8"
eindex = 0
eul1 = op2.trmbu.eulers[subcase_id][eltype][0, eindex, 1] # * pi / 180.0
eul2 = op2.trmbu.eulers[subcase_id][eltype][0, eindex, 2] # * pi / 180.0
eul3 = op2.trmbu.eulers[subcase_id][eltype][0, eindex, 3] # * pi / 180.0
eul1d = eul1 * 180 / pi
eul2d = eul2 * 180 / pi
eul3d = eul3 * 180 / pi
c1 = cos(eul1)
s1 = sin(eul1)
c2 = cos(eul2)
s2 = sin(eul2)
c3 = cos(eul3)
s3 = sin(eul3)
unit_t_mats = {"Z1": np.array([[c1, -s1, 0],
[s1, c1, 0],
[0, 0, 1]]),
"Z2": np.array([[c2, -s2, 0],
[s2, c2, 0],
[0, 0, 1]]),
"Z3": np.array([[c3, -s3, 0],
[s3, c3, 0],
[0, 0, 1]]),
"Y1": np.array([[c1, 0, s1],
[0, 1, 0],
[-s1, 0, c1]]),
"Y2": np.array([[c2, 0, s2],
[0, 1, 0],
[-s2, 0, c2]]),
"Y3": np.array([[c3, 0, s3],
[0, 1, 0],
[-s3, 0, c3]]),
"X1": np.array([[1, 0, 0],
[0, c1, -s1],
[0, s1, c1]]),
"X2": np.array([[1, 0, 0],
[0, c2, -s2],
[0, s2, c2]]),
"X3": np.array([[1, 0, 0],
[0, c3, -s3],
[0, s3, c3]])}
unit_t_mats["X1T"] = unit_t_mats["X1"].T
unit_t_mats["X2T"] = unit_t_mats["X2"].T
unit_t_mats["X3T"] = unit_t_mats["X3"].T
unit_t_mats["Y1T"] = unit_t_mats["Y1"].T
unit_t_mats["Y2T"] = unit_t_mats["Y2"].T
unit_t_mats["Y3T"] = unit_t_mats["Y3"].T
unit_t_mats["Z1T"] = unit_t_mats["Z1"].T
unit_t_mats["Z2T"] = unit_t_mats["Z2"].T
unit_t_mats["Z3T"] = unit_t_mats["Z3"].T
# Create all possible transformation matrix combinations
t_mats = {}
for key1, mat1 in unit_t_mats.items():
for key2, mat2 in unit_t_mats.items():
for key3, mat3 in unit_t_mats.items():
t_mats[" @ ".join((key1, key2, key3))] = mat1 @ mat2 @ mat3
# Find closest transformation matrix to CSTM (gbt)
best_mat = 0
best_key = 0
best_value = 1e10
for key, t_mat in t_mats.items():
diff = np.sum(np.abs((gbt - t_mat)))
if diff < best_value:
best_value = diff
best_key = key
best_mat = t_mat
# Print best result
print(f"CSTM Transformation matrix:\n\n{gbt}\n")
print(f"Euler angle transformation matrix:\n\n{best_mat}\n")
print(f"Best combination: {best_key}\n")
print(f"Euler angles: {eul1d:.3f} deg, {eul2d:.3f} deg, {eul3d:.3f} deg")
Results
CSTM Transformation matrix:
[[ 0.84710067 -0.18005681 0.5 ]
[ 0.49284317 0.61814692 -0.61237244]
[-0.19881163 0.76516268 0.61237244]]
Euler angle transformation matrix:
[[ 0.82329392 -0.19347388 0.53362438]
[ 0.49115223 0.71407289 -0.49886812]
[-0.28452875 0.67280589 0.68291699]]
Best combination: Z1T @ X1T @ Z3
Euler angles: -46.928 deg, -23.806 deg, -22.923 deg
That doesn't seem to work for the moment. I attached the test model to this comment in case you're interested.
I'll keep you updated should some 'breakthrough' happen. :)
Best regards,
Ben test.zip
Small note: I'm quite sure that the provided euler angles in TRMBU are for rotations about the X, Y and Z axis, respectively (and not ZXZ or anything else). If I use a single non-zero euler angle in the pre-processor, it always turns up in the TRMBU data as a single non-zero angle along the respective axis.
ZXZ is definitely a possibility. I’d just try the n! case and I guess 2n! if you include the transpose..
Can you attach a sample model. I want to go poke at the DMAP and see if there’s anything obvious. You can add DIAG,14 after the SOL x to dump the core code. It’s cryptic, but you can search for TRMBU and look at the code around it and print certain matrices.
Also dumb question? Is there a specific theory guide for the solution? Rotations seem like a really basic question.
On Wed, Feb 16, 2022 at 3:31 AM ben.vanbavel @.***> wrote:
Small note: I'm quite sure that the provided euler angles in TRMBU are for rotations about the X, Y and Z axis, respectively (and not ZXZ or anything else). If I use a single non-zero euler angle in the pre-processor, it always turns up in the TRMBU data as a single non-zero angle along the respective axis.
— Reply to this email directly, view it on GitHub https://github.com/SteveDoyle2/pyNastran/issues/684#issuecomment-1041393168, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAICUWNPAXOJ3JNERT33DNTU3ODH3ANCNFSM5OHLMEKQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.
You are receiving this because you commented.Message ID: @.***>
Hi Steve,
I think you might've missed my previous comment (unless I'm confused), there I elaborate a bit about brute-forcing all possible euler combinations. I also attached a small model to that comment (test.zip).
I didn't try the trick with DIAG,14 yet. I'll definitely try that as well, thanks!
From the documentation it seems the TRMBU/TRMBD tables were added as of NX Nastran 12 (in the NX Nastran Release Guide they're listed under the section 'new data blocks') but not much else that elaborates about the theory behind it. There's a user's guide for SOL 401/402, but it doesn't explain the TRMBU/TRMBD tables. I'll try searching more in depth this week and I'll let you know if I find anything!
Best regards,
Ben
In addition to my last comment: Here's a copy of the previous test model (a single CHEXA8 element that rotates through time with a simple load), but with the DIAG,14 line added to the .dat file. I also included the output f06 and op2 files.
It's indeed pretty cryptic, but I'll try looking for clues.
Hi Steve,
A small update: I believe TRMBU/TRMBD contains the data in the axis-angle representation instead of the classical XYZ euler angles. That would explain the difference in angles between the pre-processor and the TRMBU/D tables. I've been able to confirm that I can retrieve the same values as those in TRMBU/D that way.
I still have to write a small example, but I'll share it once I have done so.
Best regards,
Ben
That’s good news. The DMAP was cryptic even for DMAP. I was only ever referenced in the restart, so I think it was called something else in the main code.
I ran some tests and got a more general reader. To add a new element, you just need the number of nodes. For a CQUAD8, obviously you use 4 nodes.
Just for future reference, it’s probably just a linear interpolation for the mudslide nodes either before or after the rotations are applied (and that question is probably not solvable without a one element problem cause the values are going to be close unless your strains are very large).
On Thu, Feb 17, 2022 at 5:29 AM ben.vanbavel @.***> wrote:
Hi Steve,
A small update: I believe TRMBU/TRMBD contains the data in the axis-angle representation instead of the classical XYZ euler angles. That would explain the difference in angles between the pre-processor and the TRMBU/D tables. I've been able to confirm that I can retrieve the same values as those in TRMBU/D that way.
I still have to write a small example, but I'll share it once I have done so.
Best regards,
Ben
— Reply to this email directly, view it on GitHub https://github.com/SteveDoyle2/pyNastran/issues/684#issuecomment-1042950066, or unsubscribe https://github.com/notifications/unsubscribe-auth/AAICUWIXZANKHI3NU7MLUVDU3TZ2FANCNFSM5OHLMEKQ . Triage notifications on the go with GitHub Mobile for iOS https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675 or Android https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub.
You are receiving this because you commented.Message ID: @.***>
Hi Steve,
Thanks! Here's a small example code that I hope might be of interest to you, and others that might stumble upon this issue in the future while working with the TRMBU/D tables. It shows how to transform the stresses of a single node of a single element from the output coordinate system to the basic/absolute coordinate system. Results were checked with the test model. The TRMBU/D data is indeed using the axis-angle representation, instead of three XYZ euler angles (as it's written somewhat confusingly in the documentation).
In order to be complete, these are the edited pyNastran files (I'll use your more general reader in the future): editedPynastran.zip
Once again, a big thank you for all the help!
Best regards,
Ben
import numpy as np
from numpy import pi, cos, sin
from pyNastran.op2.op2 import OP2
from pyNastran.op2.tables.oes_stressStrain.oes import RealSolidStressArrayNx
from pyNastran.op2.result_objects.op2_results import CSTM, TRMBU, TRMBD
def get_stresses(op2: OP2, nonlinear: bool, time_id, e_id, grid_id):
"""
Get stress results of a single node of a single element at some point in time
"""
# Get stress results
subcase_id = 1
chexa_stress = op2.chexa_stress[subcase_id] # type: RealSolidStressArrayNx
ntimes = chexa_stress.ntimes
nelements = chexa_stress.nelements
# Reshape results (Time, Element, Node, Stress component)
if nonlinear is True:
reshaped_results = chexa_stress.data.reshape(ntimes, nelements, 8, 7)
else:
# Remove elemental stresses
reshaped_results = np.delete(chexa_stress.data, np.arange(0, chexa_stress.data.shape[1], 9), axis=1)
reshaped_results = reshaped_results.reshape(1, nelements, 8, 10)
stresses = reshaped_results[time_id, e_id, grid_id]
# Create 3x3 stress matrix
stress_mat = np.empty([3, 3], dtype="float64")
stress_mat[0, 0] = stresses[0] # XX
stress_mat[0, 1] = stresses[3] # XY
stress_mat[0, 2] = stresses[5] # XZ
stress_mat[1, 0] = stresses[3] # XY
stress_mat[1, 1] = stresses[1] # YY
stress_mat[1, 2] = stresses[4] # YZ
stress_mat[2, 0] = stresses[5] # XZ
stress_mat[2, 1] = stresses[4] # YZ
stress_mat[2, 2] = stresses[2] # ZZ
return stress_mat
def test_cstm(op2: OP2, csys_id):
"""Create global to basic csys transformation matrix using CSTM table"""
cstm = op2.cstm # type: CSTM
gbt = np.empty([3, 3], dtype="float64")
gbt[0, 0] = cstm.data[csys_id, cstm.headers["T11"]]
gbt[0, 1] = cstm.data[csys_id, cstm.headers["T12"]]
gbt[0, 2] = cstm.data[csys_id, cstm.headers["T13"]]
gbt[1, 0] = cstm.data[csys_id, cstm.headers["T21"]]
gbt[1, 1] = cstm.data[csys_id, cstm.headers["T22"]]
gbt[1, 2] = cstm.data[csys_id, cstm.headers["T23"]]
gbt[2, 0] = cstm.data[csys_id, cstm.headers["T31"]]
gbt[2, 1] = cstm.data[csys_id, cstm.headers["T32"]]
gbt[2, 2] = cstm.data[csys_id, cstm.headers["T33"]]
print("CSTM")
print(gbt)
def test_trmbu(op2: OP2, e_id):
"""
Shows how to create transformation matrix from information in TRMBU
"""
#
subcase_id = 0
eltype = "CHEXA8"
# Retrieve axis-angle from TRMBU table
aa = op2.trmbu.eulers[subcase_id][eltype][0, e_id, 1:4]
# Normalize rotation vector with rotation angle (magnitude == rotation angle)
angle = np.linalg.norm(aa)
aa = aa / angle
# Create rotation matrix
x = aa[0]
y = aa[1]
z = aa[2]
ca = cos(angle)
sa = sin(angle)
t = 1 - ca
T = np.empty([3, 3])
T[0, 0] = t*x*x + ca
T[0, 1] = t*x*y - z*sa
T[0, 2] = t*x*z + y*sa
T[1, 0] = t*x*y + z*sa
T[1, 1] = t*y*y + ca
T[1, 2] = t*y*z - x*sa
T[2, 0] = t*x*z - y*sa
T[2, 1] = t*y*z + x*sa
T[2, 2] = t*z*z + ca
# Compare result with CSTM table
print("TRMBU Transformation matrix")
print(T.T)
def test_trmbd(op2: OP2, stresses, time_id, e_id, grid_id):
"""
Shows how to create transformation matrix from information in TRMBD
"""
#
subcase_id = 1
eltype = "CHEXA8"
# Retrieve axis-angle from TRMBU table
eulx = op2.trmbd.eulersx[subcase_id][eltype][time_id, e_id, grid_id]
euly = op2.trmbd.eulersy[subcase_id][eltype][time_id, e_id, grid_id]
eulz = op2.trmbd.eulersz[subcase_id][eltype][time_id, e_id, grid_id]
# axis-angle in rads
aa = np.array([eulx, euly, eulz])
# Normalize rotation vector with rotation angle (magnitude == rotation angle)
angle = np.linalg.norm(aa)
aa = aa / angle
# Create rotation matrix
x = aa[0]
y = aa[1]
z = aa[2]
ca = cos(angle)
sa = sin(angle)
t = 1 - ca
T = np.empty([3, 3])
T[0, 0] = t * x * x + ca
T[0, 1] = t * x * y - z * sa
T[0, 2] = t * x * z + y * sa
T[1, 0] = t * x * y + z * sa
T[1, 1] = t * y * y + ca
T[1, 2] = t * y * z - x * sa
T[2, 0] = t * x * z - y * sa
T[2, 1] = t * y * z + x * sa
T[2, 2] = t * z * z + ca
# Print results
print("TRMBD Transformation matrix")
print(T.T)
# Print transformed stresses
stresses_absolute = T.T @ stresses @ T
print("Transformed stresses")
print(stresses_absolute)
if __name__ == "__main__":
# Create OP2 path
op2_filepath = ""
op2_filename = ""
op2_full_path = "/".join((op2_filepath, op2_filename))
# Read OP2 results
op2_results = OP2(mode="nx", debug=True)
op2_results.read_op2(op2_full_path)
# Get stresses from OP2
nonlinear = True
timeindex = 1
eindex = 0
nodeindex = 0
stress = get_stresses(op2=op2_results, nonlinear=nonlinear, time_id=timeindex, e_id=eindex, grid_id=nodeindex)
# Test CSTM
test_cstm(op2_results, csys_id=1)
# Test TRMBU/D
if nonlinear:
test_trmbu(op2_results, e_id=eindex)
test_trmbd(op2_results, stress, time_id=timeindex, e_id=eindex, grid_id=nodeindex)
Hi Steve, I noticed that the code for TRMBD/TRMBU tables has been implemented in the op2 reader, but that it was still commented out (and hence these tables are not read by the op2 reader). I'll send a pull request soon that uncomments the code alongside some small changes. Hope it helps!