Gold Standard Test: Shaw et al. (2009)¶
This Jupyter notebook is part of the HBM4VT Validation Catalogue.
© 2024, OpenVT Organization OVTO
Available openly under Creative Commons Attribution 4.0 International License 
Experiments by Shaw et al. (2009)¶
Summary¶
The "Gold Standard" sled test was first published by Shaw et al. (2009).
Shaw, G., Parent, D., Purtsezov, S., Lessley, D., Crandall, J., Kent, R., Guillemot, H., Ridella, S. A., Takhounts, E., & Martin, P. (2009). Impact response of restrained PMHS in frontal sled tests: Skeletal deformation patterns under seat belt loading. Stapp Car Crash Journal, 53, 1–48. https://doi.org/10.4271/2009-22-0001
However, variants using the same sled setup using different speeds and belt load limiters have been performed. For details, please refer to the README.md and the documentation.
Example animations illustrating the sled test in a simulation with VIVA+ 50F are provided below. Note that while the sled moves in the simulations, a viewpoint moving with the sled was chosen for the animations.
When sharing your simulation results (TU Graz cloud, uploaded data is only accessible for the TU Graz developer group), please add animations (e.g.in *.avi format) of the views shown above to the zip-folder. These can be very helpful to analyze potential differences between different HBMs.
Experiments¶
Information on the subjects/specimens¶
Please refer to the "PMHS reference data overview" section in the README.md for detailed info on the tested PMHS for each of the three Gold Standard test setups.
Responses recorded¶
The reference values from the tests were obtained from the NHTSA biomechanics database of the tests listed in the README.md. There you can also find the links to the original data.
How to use this Notebook¶
- Create a csv file including the raw data and a settings file according to the format of the files provided on https://openvt.eu/EuroNCAP/hbm4vt-validation-catalogue/ls-dyna/shaw_2009/-/tree/main/data/processed/sim_results?ref_type=heads
- Make sure the file arc_length_ellipse_score.py is provided in the folder the notebook is executed from (download it from https://openvt.eu/EuroNCAP/hbm4vt-validation-catalogue/ls-dyna/shaw_2009/-/tree/main/data/metadata?ref_type=heads)
- Define the paths and file names in the following block
- Run the notebook and check your results for plausibility.
In [ ]:
Copied!
# Define paths and file names
result_directiory = "Insert the path to your csv file here"
csv_file_name = "Insert the name of your csv file here" # e.g.: 'Gold_Standard_2_Dynasaur_diagram_output_THUMS_v4.1_50M.csv'
setting_file_name = "Insert the name of your settings file here" # e.g.: 'Gold_Standard_2_settings_THUMS_v4.1_50M.txt'
experimental_data_dir = "Insert the path to the experimental data here" # Download the experimental data from https://openvt.eu/EuroNCAP/hbm4vt-validation-catalogue/ls-dyna/shaw_2009/-/tree/main/data/experiment?ref_type=heads
reference_data_dir = "Insert the path to the experimental data here" # Download the experimental data from https://openvt.eu/EuroNCAP/hbm4vt-validation-catalogue/ls-dyna/shaw_2009/-/tree/main/data/reference?ref_type=heads
# Define rib criterion
rib_criterion = 'Larsson2021'
# Define paths and file names
result_directiory = "Insert the path to your csv file here"
csv_file_name = "Insert the name of your csv file here" # e.g.: 'Gold_Standard_2_Dynasaur_diagram_output_THUMS_v4.1_50M.csv'
setting_file_name = "Insert the name of your settings file here" # e.g.: 'Gold_Standard_2_settings_THUMS_v4.1_50M.txt'
experimental_data_dir = "Insert the path to the experimental data here" # Download the experimental data from https://openvt.eu/EuroNCAP/hbm4vt-validation-catalogue/ls-dyna/shaw_2009/-/tree/main/data/experiment?ref_type=heads
reference_data_dir = "Insert the path to the experimental data here" # Download the experimental data from https://openvt.eu/EuroNCAP/hbm4vt-validation-catalogue/ls-dyna/shaw_2009/-/tree/main/data/reference?ref_type=heads
# Define rib criterion
rib_criterion = 'Larsson2021'
In [ ]:
Copied!
setup_dict={'40kph_notForceLimited': 'Gold Standard 1',
'30kph_3kN': 'Gold Standard 2',
'30kph_2kN': 'Gold Standard 2.1'}
belt_angle_range_dict={'40kph_notForceLimited': {'low': 24, 'high': 29},
'30kph_3kN': {'low': 25, 'high': 26},
'30kph_2kN': {'low': 23, 'high': 30}}
file_settings = [os.path.join(result_directiory, setting_file_name)]
# Read the content of the first settings file
with open(file_settings[0], 'r') as fp:
lines = fp.readlines()
unit_line = lines[2]
if not unit_line.startswith(' "units"'):
raise ValueError(f'Line does not start with "units:" but {unit_line}.')
unit = unit_line.split('":')[1].strip().replace('"','').replace(',','')
HBM_line = lines[1]
if not HBM_line.startswith(' "HBM"'):
raise ValueError(f'Line does not start with "units:" but {HBM_line}.')
HBM = HBM_line.split('":')[1].strip().replace('"','').replace(',','')
t_settle_line = lines[6]
if not t_settle_line.startswith(' "t_settle"'):
raise ValueError(f'Line does not start with "t_settle:" but {t_settle_line}.')
t_settle1 = float(t_settle_line.split('":')[1].strip().replace('"','').replace(',',''))
setup_type_line = lines[10]
if not setup_type_line.startswith(' "setup_type"'):
raise ValueError(f'Line does not start with "setup_type:" but {setup_type_line}.')
setup_type = setup_type_line.split('":')[1].strip().replace('"','').replace(',','')
unit_system_to_time = {"s_mm_ton": "second", "ms_mm_kg": "millisecond"}
t_settle = (float(t_settle1) * ureg("seconds")).to(unit_system_to_time[unit]).magnitude
setup_dict={'40kph_notForceLimited': 'Gold Standard 1',
'30kph_3kN': 'Gold Standard 2',
'30kph_2kN': 'Gold Standard 2.1'}
belt_angle_range_dict={'40kph_notForceLimited': {'low': 24, 'high': 29},
'30kph_3kN': {'low': 25, 'high': 26},
'30kph_2kN': {'low': 23, 'high': 30}}
file_settings = [os.path.join(result_directiory, setting_file_name)]
# Read the content of the first settings file
with open(file_settings[0], 'r') as fp:
lines = fp.readlines()
unit_line = lines[2]
if not unit_line.startswith(' "units"'):
raise ValueError(f'Line does not start with "units:" but {unit_line}.')
unit = unit_line.split('":')[1].strip().replace('"','').replace(',','')
HBM_line = lines[1]
if not HBM_line.startswith(' "HBM"'):
raise ValueError(f'Line does not start with "units:" but {HBM_line}.')
HBM = HBM_line.split('":')[1].strip().replace('"','').replace(',','')
t_settle_line = lines[6]
if not t_settle_line.startswith(' "t_settle"'):
raise ValueError(f'Line does not start with "t_settle:" but {t_settle_line}.')
t_settle1 = float(t_settle_line.split('":')[1].strip().replace('"','').replace(',',''))
setup_type_line = lines[10]
if not setup_type_line.startswith(' "setup_type"'):
raise ValueError(f'Line does not start with "setup_type:" but {setup_type_line}.')
setup_type = setup_type_line.split('":')[1].strip().replace('"','').replace(',','')
unit_system_to_time = {"s_mm_ton": "second", "ms_mm_kg": "millisecond"}
t_settle = (float(t_settle1) * ureg("seconds")).to(unit_system_to_time[unit]).magnitude
In [ ]:
Copied!
def get_unit_from_column_name(pd_df):
unit = list(pd_df.columns)
assert len(unit) == 1
assert type(unit[0]) == str
return unit[0]
def convert_to_unit(pd_df, unit, offset=None):
data = pd_df.dropna().to_numpy()
# Offset is done before conversion.
if offset is not None:
data = data[offset:] - data[offset]
data_unit = get_unit_from_column_name(pd_df)
data = data * ureg(data_unit)
return data.to(ureg(unit)).magnitude
if setup_type == '30kph_2kN':
sim_data_no_occupant = pd.read_csv(os.path.join(reference_data_dir, 'load_cell_forces_without_occupant_G2.1.csv'), delimiter=';', na_values='-', header = [0,1,2,3])
sim_data_no_occupant.sort_index(axis=1, inplace=True)
def subtract_no_occupant_data(idx_1, idx_2, idx_3, offset):
time = convert_to_unit(sim_data_results[idx_1][idx_2]['time'], unit_time_plot, offset=offset)
# No offset for simulation data with no occupant.
time_no_occupant = convert_to_unit(sim_data_no_occupant[idx_1][idx_2]['time'], unit_time_plot)
force_no_occupant = convert_to_unit(sim_data_no_occupant[idx_1][idx_2][idx_3], unit_force_plot)
# Extrapolate with first and last value respectively (if necessary).
f = interpolate.interp1d(time_no_occupant[:, 0], force_no_occupant[:, 0], bounds_error=False, fill_value=(force_no_occupant[0, 0], force_no_occupant[-1, 0]))
force_subtracted = convert_to_unit(sim_data_results[idx_1][idx_2][idx_3], unit_force_plot, offset=offset) - f(time)
return time, force_subtracted
def get_unit_from_column_name(pd_df):
unit = list(pd_df.columns)
assert len(unit) == 1
assert type(unit[0]) == str
return unit[0]
def convert_to_unit(pd_df, unit, offset=None):
data = pd_df.dropna().to_numpy()
# Offset is done before conversion.
if offset is not None:
data = data[offset:] - data[offset]
data_unit = get_unit_from_column_name(pd_df)
data = data * ureg(data_unit)
return data.to(ureg(unit)).magnitude
if setup_type == '30kph_2kN':
sim_data_no_occupant = pd.read_csv(os.path.join(reference_data_dir, 'load_cell_forces_without_occupant_G2.1.csv'), delimiter=';', na_values='-', header = [0,1,2,3])
sim_data_no_occupant.sort_index(axis=1, inplace=True)
def subtract_no_occupant_data(idx_1, idx_2, idx_3, offset):
time = convert_to_unit(sim_data_results[idx_1][idx_2]['time'], unit_time_plot, offset=offset)
# No offset for simulation data with no occupant.
time_no_occupant = convert_to_unit(sim_data_no_occupant[idx_1][idx_2]['time'], unit_time_plot)
force_no_occupant = convert_to_unit(sim_data_no_occupant[idx_1][idx_2][idx_3], unit_force_plot)
# Extrapolate with first and last value respectively (if necessary).
f = interpolate.interp1d(time_no_occupant[:, 0], force_no_occupant[:, 0], bounds_error=False, fill_value=(force_no_occupant[0, 0], force_no_occupant[-1, 0]))
force_subtracted = convert_to_unit(sim_data_results[idx_1][idx_2][idx_3], unit_force_plot, offset=offset) - f(time)
return time, force_subtracted
Results¶
Shoulder belt angle¶
In [ ]:
Copied!
offset = np.argmax(sim_data_results['TESTBED']['End_of_2D_shoulder_belt_X_coordinate_over_time']['time'] >= t_settle)
assert offset > 0
D_ring_point = [sim_data_results['TESTBED']['Slipring_shoulder_belt_X_coordinate_over_time']['x-coordinate'].iloc[offset][0], sim_data_results['TESTBED']['Slipring_shoulder_belt_Z_coordinate_over_time']['z-coordinate'].iloc[offset][0]]
shoulder_belt_point = [sim_data_results['TESTBED']['End_of_2D_shoulder_belt_X_coordinate_over_time']['x-coordinate'].iloc[offset][0], sim_data_results['TESTBED']['End_of_2D_shoulder_belt_Z_coordinate_over_time']['z-coordinate'].iloc[offset][0]]
horizontal_distance = math.dist([D_ring_point[0]], [shoulder_belt_point[0]])
vertical_distance = math.dist([D_ring_point[1]], [shoulder_belt_point[1]])
xz_distance = math.dist(D_ring_point, shoulder_belt_point)
shoulder_belt_angle = math.degrees(math.acos(horizontal_distance/xz_distance))
print("The shoulder belt angle after settling is "+"{:.1f}".format((shoulder_belt_angle))+"°.")
print("It should be between "+str(belt_angle_range_dict[setup_type]['low'])+"° and "+str(belt_angle_range_dict[setup_type]['high'])+"° for the "+str(setup_dict[setup_type])+" load case.")
assert shoulder_belt_angle >= 24, "Shoulder belt angle is too small. Adjust your setup to increase it (cf. Step 8 in the documentation)."
assert shoulder_belt_angle <= 29, "Shoulder belt angle is too large. Adjust your setup to reduce it (cf. Step 8 in the documentation)."
offset = np.argmax(sim_data_results['TESTBED']['End_of_2D_shoulder_belt_X_coordinate_over_time']['time'] >= t_settle)
assert offset > 0
D_ring_point = [sim_data_results['TESTBED']['Slipring_shoulder_belt_X_coordinate_over_time']['x-coordinate'].iloc[offset][0], sim_data_results['TESTBED']['Slipring_shoulder_belt_Z_coordinate_over_time']['z-coordinate'].iloc[offset][0]]
shoulder_belt_point = [sim_data_results['TESTBED']['End_of_2D_shoulder_belt_X_coordinate_over_time']['x-coordinate'].iloc[offset][0], sim_data_results['TESTBED']['End_of_2D_shoulder_belt_Z_coordinate_over_time']['z-coordinate'].iloc[offset][0]]
horizontal_distance = math.dist([D_ring_point[0]], [shoulder_belt_point[0]])
vertical_distance = math.dist([D_ring_point[1]], [shoulder_belt_point[1]])
xz_distance = math.dist(D_ring_point, shoulder_belt_point)
shoulder_belt_angle = math.degrees(math.acos(horizontal_distance/xz_distance))
print("The shoulder belt angle after settling is "+"{:.1f}".format((shoulder_belt_angle))+"°.")
print("It should be between "+str(belt_angle_range_dict[setup_type]['low'])+"° and "+str(belt_angle_range_dict[setup_type]['high'])+"° for the "+str(setup_dict[setup_type])+" load case.")
assert shoulder_belt_angle >= 24, "Shoulder belt angle is too small. Adjust your setup to increase it (cf. Step 8 in the documentation)."
assert shoulder_belt_angle <= 29, "Shoulder belt angle is too large. Adjust your setup to reduce it (cf. Step 8 in the documentation)."
Head displacement plots¶
T1 displacement plots¶
In [ ]:
Copied!
fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(18, 9))
ax[1, 2].remove()
fig.suptitle(setup_dict[setup_type])
offset = np.argmax(sim_data_results['HBM']['T1_X_coordinate_over_time']['time'] >= t_settle)
assert offset > 0
T1_x_displacement = convert_to_unit(sim_data_results['HBM']['T1_X_coordinate_over_time']['x-coordinate'], unit_deflection_plot, offset=offset) - convert_to_unit(sim_data_results['TESTBED']['Testbed_X_coordinate_over_time']['x-coordinate'], unit_deflection_plot, offset=offset)
T1_y_displacement = convert_to_unit(sim_data_results['HBM']['T1_Y_coordinate_over_time']['y-coordinate'], unit_deflection_plot, offset=offset) - convert_to_unit(sim_data_results['TESTBED']['Testbed_Y_coordinate_over_time']['y-coordinate'], unit_deflection_plot, offset=offset)
T1_z_displacement = convert_to_unit(sim_data_results['HBM']['T1_Z_coordinate_over_time']['z-coordinate'], unit_deflection_plot, offset=offset) - convert_to_unit(sim_data_results['TESTBED']['Testbed_Z_coordinate_over_time']['z-coordinate'], unit_deflection_plot, offset=offset)
time = convert_to_unit(sim_data_results['HBM']['T1_X_coordinate_over_time']['time'], unit_time_plot, offset=offset)
sim_results_for_score['t1'] = [time, T1_x_displacement, T1_y_displacement, T1_z_displacement]
ax[0, 0].plot(convert_to_unit(sim_data_results['HBM']['T1_X_coordinate_over_time']['time'], unit_time_plot, offset=offset), T1_x_displacement, label=HBM)
for data_exp in data_exp_T1_dx:
ax[0, 0].plot(data_exp[:, 0], -(data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[0, 1].plot(convert_to_unit(sim_data_results['HBM']['T1_Y_coordinate_over_time']['time'], unit_time_plot, offset=offset), T1_y_displacement)
for data_exp in data_exp_T1_dy:
ax[0, 1].plot(data_exp[:, 0], (data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[0, 2].plot(convert_to_unit(sim_data_results['HBM']['T1_Z_coordinate_over_time']['time'], unit_time_plot, offset=offset), T1_z_displacement)
for data_exp in data_exp_T1_dz:
ax[0, 2].plot(data_exp[:, 0], -(data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[1, 0].plot(T1_x_displacement, T1_z_displacement)
for data_exp in data_exp_T1_dx_dz:
ax[1, 0].plot(-(data_exp[:, 0] - data_exp[0, 0]), -(data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[1, 1].plot(T1_y_displacement, T1_z_displacement)
for data_exp in data_exp_T1_dy_dz:
ax[1, 1].plot((data_exp[:, 0] - data_exp[0, 0]), -(data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[0, 0].set(title="T1 X-Displacement", xlabel="Time / ms", ylabel="Displacement / mm")
ax[0, 1].set(title="T1 Y-Displacement", xlabel="Time / ms", ylabel="Displacement / mm")
ax[0, 2].set(title="T1 Z-Displacement", xlabel="Time / ms", ylabel="Displacement / mm")
ax[1, 0].set(title="T1 Trajectory", xlabel="X Displacement / mm", ylabel="Z Displacement / mm")
ax[1, 1].set(title="T1 Trajectory", xlabel="Y Displacement / mm", ylabel="Z Displacement / mm")
fig.tight_layout()
fig.legend(loc='lower right', bbox_to_anchor=(0.9, 0.25))
plt.show()
fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(18, 9))
ax[1, 2].remove()
fig.suptitle(setup_dict[setup_type])
offset = np.argmax(sim_data_results['HBM']['T1_X_coordinate_over_time']['time'] >= t_settle)
assert offset > 0
T1_x_displacement = convert_to_unit(sim_data_results['HBM']['T1_X_coordinate_over_time']['x-coordinate'], unit_deflection_plot, offset=offset) - convert_to_unit(sim_data_results['TESTBED']['Testbed_X_coordinate_over_time']['x-coordinate'], unit_deflection_plot, offset=offset)
T1_y_displacement = convert_to_unit(sim_data_results['HBM']['T1_Y_coordinate_over_time']['y-coordinate'], unit_deflection_plot, offset=offset) - convert_to_unit(sim_data_results['TESTBED']['Testbed_Y_coordinate_over_time']['y-coordinate'], unit_deflection_plot, offset=offset)
T1_z_displacement = convert_to_unit(sim_data_results['HBM']['T1_Z_coordinate_over_time']['z-coordinate'], unit_deflection_plot, offset=offset) - convert_to_unit(sim_data_results['TESTBED']['Testbed_Z_coordinate_over_time']['z-coordinate'], unit_deflection_plot, offset=offset)
time = convert_to_unit(sim_data_results['HBM']['T1_X_coordinate_over_time']['time'], unit_time_plot, offset=offset)
sim_results_for_score['t1'] = [time, T1_x_displacement, T1_y_displacement, T1_z_displacement]
ax[0, 0].plot(convert_to_unit(sim_data_results['HBM']['T1_X_coordinate_over_time']['time'], unit_time_plot, offset=offset), T1_x_displacement, label=HBM)
for data_exp in data_exp_T1_dx:
ax[0, 0].plot(data_exp[:, 0], -(data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[0, 1].plot(convert_to_unit(sim_data_results['HBM']['T1_Y_coordinate_over_time']['time'], unit_time_plot, offset=offset), T1_y_displacement)
for data_exp in data_exp_T1_dy:
ax[0, 1].plot(data_exp[:, 0], (data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[0, 2].plot(convert_to_unit(sim_data_results['HBM']['T1_Z_coordinate_over_time']['time'], unit_time_plot, offset=offset), T1_z_displacement)
for data_exp in data_exp_T1_dz:
ax[0, 2].plot(data_exp[:, 0], -(data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[1, 0].plot(T1_x_displacement, T1_z_displacement)
for data_exp in data_exp_T1_dx_dz:
ax[1, 0].plot(-(data_exp[:, 0] - data_exp[0, 0]), -(data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[1, 1].plot(T1_y_displacement, T1_z_displacement)
for data_exp in data_exp_T1_dy_dz:
ax[1, 1].plot((data_exp[:, 0] - data_exp[0, 0]), -(data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[0, 0].set(title="T1 X-Displacement", xlabel="Time / ms", ylabel="Displacement / mm")
ax[0, 1].set(title="T1 Y-Displacement", xlabel="Time / ms", ylabel="Displacement / mm")
ax[0, 2].set(title="T1 Z-Displacement", xlabel="Time / ms", ylabel="Displacement / mm")
ax[1, 0].set(title="T1 Trajectory", xlabel="X Displacement / mm", ylabel="Z Displacement / mm")
ax[1, 1].set(title="T1 Trajectory", xlabel="Y Displacement / mm", ylabel="Z Displacement / mm")
fig.tight_layout()
fig.legend(loc='lower right', bbox_to_anchor=(0.9, 0.25))
plt.show()
T8 displacement plots¶
In [ ]:
Copied!
fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(18, 9))
ax[1, 2].remove()
fig.suptitle(setup_dict[setup_type])
offset = np.argmax(sim_data_results['HBM']['T8_X_coordinate_over_time']['time'] >= t_settle)
assert offset > 0
T8_x_displacement = convert_to_unit(sim_data_results['HBM']['T8_X_coordinate_over_time']['x-coordinate'], unit_deflection_plot, offset=offset) - convert_to_unit(sim_data_results['TESTBED']['Testbed_X_coordinate_over_time']['x-coordinate'], unit_deflection_plot, offset=offset)
T8_y_displacement = convert_to_unit(sim_data_results['HBM']['T8_Y_coordinate_over_time']['y-coordinate'], unit_deflection_plot, offset=offset) - convert_to_unit(sim_data_results['TESTBED']['Testbed_Y_coordinate_over_time']['y-coordinate'], unit_deflection_plot, offset=offset)
T8_z_displacement = convert_to_unit(sim_data_results['HBM']['T8_Z_coordinate_over_time']['z-coordinate'], unit_deflection_plot, offset=offset) - convert_to_unit(sim_data_results['TESTBED']['Testbed_Z_coordinate_over_time']['z-coordinate'], unit_deflection_plot, offset=offset)
time = convert_to_unit(sim_data_results['HBM']['T8_X_coordinate_over_time']['time'], unit_time_plot, offset=offset)
sim_results_for_score['t8'] = [time, T8_x_displacement, T8_y_displacement, T8_z_displacement]
ax[0, 0].plot(convert_to_unit(sim_data_results['HBM']['T8_X_coordinate_over_time']['time'], unit_time_plot, offset=offset), T8_x_displacement, label=HBM)
for data_exp in data_exp_T8_dx:
ax[0, 0].plot(data_exp[:, 0], -(data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[0, 1].plot(convert_to_unit(sim_data_results['HBM']['T8_Y_coordinate_over_time']['time'], unit_time_plot, offset=offset), T8_y_displacement)
for data_exp in data_exp_T8_dy:
ax[0, 1].plot(data_exp[:, 0], (data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[0, 2].plot(convert_to_unit(sim_data_results['HBM']['T8_Z_coordinate_over_time']['time'], unit_time_plot, offset=offset), T8_z_displacement)
for data_exp in data_exp_T8_dz:
ax[0, 2].plot(data_exp[:, 0], -(data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[1, 0].plot(T8_x_displacement, T8_z_displacement)
for data_exp in data_exp_T8_dx_dz:
ax[1, 0].plot(-(data_exp[:, 0] - data_exp[0, 0]), -(data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[1, 1].plot(T8_y_displacement, T8_z_displacement)
for data_exp in data_exp_T8_dy_dz:
ax[1, 1].plot((data_exp[:, 0] - data_exp[0, 0]), -(data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[0, 0].set(title="T8 X-Displacement", xlabel="Time / ms", ylabel="Displacement / mm")
ax[0, 1].set(title="T8 Y-Displacement", xlabel="Time / ms", ylabel="Displacement / mm")
ax[0, 2].set(title="T8 Z-Displacement", xlabel="Time / ms", ylabel="Displacement / mm")
ax[1, 0].set(title="T8 Trajectory", xlabel="X Displacement / mm", ylabel="Z Displacement / mm")
ax[1, 1].set(title="T8 Trajectory", xlabel="Y Displacement / mm", ylabel="Z Displacement / mm")
fig.tight_layout()
fig.legend(loc='lower right', bbox_to_anchor=(0.9, 0.25))
plt.show()
fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(18, 9))
ax[1, 2].remove()
fig.suptitle(setup_dict[setup_type])
offset = np.argmax(sim_data_results['HBM']['T8_X_coordinate_over_time']['time'] >= t_settle)
assert offset > 0
T8_x_displacement = convert_to_unit(sim_data_results['HBM']['T8_X_coordinate_over_time']['x-coordinate'], unit_deflection_plot, offset=offset) - convert_to_unit(sim_data_results['TESTBED']['Testbed_X_coordinate_over_time']['x-coordinate'], unit_deflection_plot, offset=offset)
T8_y_displacement = convert_to_unit(sim_data_results['HBM']['T8_Y_coordinate_over_time']['y-coordinate'], unit_deflection_plot, offset=offset) - convert_to_unit(sim_data_results['TESTBED']['Testbed_Y_coordinate_over_time']['y-coordinate'], unit_deflection_plot, offset=offset)
T8_z_displacement = convert_to_unit(sim_data_results['HBM']['T8_Z_coordinate_over_time']['z-coordinate'], unit_deflection_plot, offset=offset) - convert_to_unit(sim_data_results['TESTBED']['Testbed_Z_coordinate_over_time']['z-coordinate'], unit_deflection_plot, offset=offset)
time = convert_to_unit(sim_data_results['HBM']['T8_X_coordinate_over_time']['time'], unit_time_plot, offset=offset)
sim_results_for_score['t8'] = [time, T8_x_displacement, T8_y_displacement, T8_z_displacement]
ax[0, 0].plot(convert_to_unit(sim_data_results['HBM']['T8_X_coordinate_over_time']['time'], unit_time_plot, offset=offset), T8_x_displacement, label=HBM)
for data_exp in data_exp_T8_dx:
ax[0, 0].plot(data_exp[:, 0], -(data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[0, 1].plot(convert_to_unit(sim_data_results['HBM']['T8_Y_coordinate_over_time']['time'], unit_time_plot, offset=offset), T8_y_displacement)
for data_exp in data_exp_T8_dy:
ax[0, 1].plot(data_exp[:, 0], (data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[0, 2].plot(convert_to_unit(sim_data_results['HBM']['T8_Z_coordinate_over_time']['time'], unit_time_plot, offset=offset), T8_z_displacement)
for data_exp in data_exp_T8_dz:
ax[0, 2].plot(data_exp[:, 0], -(data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[1, 0].plot(T8_x_displacement, T8_z_displacement)
for data_exp in data_exp_T8_dx_dz:
ax[1, 0].plot(-(data_exp[:, 0] - data_exp[0, 0]), -(data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[1, 1].plot(T8_y_displacement, T8_z_displacement)
for data_exp in data_exp_T8_dy_dz:
ax[1, 1].plot((data_exp[:, 0] - data_exp[0, 0]), -(data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[0, 0].set(title="T8 X-Displacement", xlabel="Time / ms", ylabel="Displacement / mm")
ax[0, 1].set(title="T8 Y-Displacement", xlabel="Time / ms", ylabel="Displacement / mm")
ax[0, 2].set(title="T8 Z-Displacement", xlabel="Time / ms", ylabel="Displacement / mm")
ax[1, 0].set(title="T8 Trajectory", xlabel="X Displacement / mm", ylabel="Z Displacement / mm")
ax[1, 1].set(title="T8 Trajectory", xlabel="Y Displacement / mm", ylabel="Z Displacement / mm")
fig.tight_layout()
fig.legend(loc='lower right', bbox_to_anchor=(0.9, 0.25))
plt.show()
L2 displacement plots¶
In [ ]:
Copied!
fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(18, 9))
ax[1, 2].remove()
fig.suptitle(setup_dict[setup_type])
offset = np.argmax(sim_data_results['HBM']['L2_X_coordinate_over_time']['time'] >= t_settle)
assert offset > 0
L2_x_displacement = convert_to_unit(sim_data_results['HBM']['L2_X_coordinate_over_time']['x-coordinate'], unit_deflection_plot, offset=offset) - convert_to_unit(sim_data_results['TESTBED']['Testbed_X_coordinate_over_time']['x-coordinate'], unit_deflection_plot, offset=offset)
L2_y_displacement = convert_to_unit(sim_data_results['HBM']['L2_Y_coordinate_over_time']['y-coordinate'], unit_deflection_plot, offset=offset) - convert_to_unit(sim_data_results['TESTBED']['Testbed_Y_coordinate_over_time']['y-coordinate'], unit_deflection_plot, offset=offset)
L2_z_displacement = convert_to_unit(sim_data_results['HBM']['L2_Z_coordinate_over_time']['z-coordinate'], unit_deflection_plot, offset=offset) - convert_to_unit(sim_data_results['TESTBED']['Testbed_Z_coordinate_over_time']['z-coordinate'], unit_deflection_plot, offset=offset)
ax[0, 0].plot(convert_to_unit(sim_data_results['HBM']['L2_X_coordinate_over_time']['time'], unit_time_plot, offset=offset), L2_x_displacement, label=HBM)
for data_exp in data_exp_L2_dx:
ax[0, 0].plot(data_exp[:, 0], -(data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[0, 1].plot(convert_to_unit(sim_data_results['HBM']['L2_Y_coordinate_over_time']['time'], unit_time_plot, offset=offset), L2_y_displacement)
for data_exp in data_exp_L2_dy:
ax[0, 1].plot(data_exp[:, 0], (data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[0, 2].plot(convert_to_unit(sim_data_results['HBM']['L2_Z_coordinate_over_time']['time'], unit_time_plot, offset=offset), L2_z_displacement)
for data_exp in data_exp_L2_dz:
ax[0, 2].plot(data_exp[:, 0], -(data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[1, 0].plot(L2_x_displacement, L2_z_displacement)
for data_exp in data_exp_L2_dx_dz:
ax[1, 0].plot(-(data_exp[:, 0] - data_exp[0, 0]), -(data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[1, 1].plot(L2_y_displacement, L2_z_displacement)
for data_exp in data_exp_L2_dy_dz:
ax[1, 1].plot((data_exp[:, 0] - data_exp[0, 0]), -(data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[0, 0].set(title="L2 X-Displacement", xlabel="Time [ms]", ylabel="Displacement [mm]")
ax[0, 1].set(title="L2 Y-Displacement", xlabel="Time [ms]", ylabel="Displacement [mm]")
ax[0, 2].set(title="L2 Z-Displacement", xlabel="Time [ms]", ylabel="Displacement [mm]")
ax[1, 0].set(title="L2 Trajectory", xlabel="X Displacement [mm]", ylabel="Z Displacement [mm]")
ax[1, 1].set(title="L2 Trajectory", xlabel="Y Displacement [mm]", ylabel="Z Displacement [mm]")
fig.tight_layout()
fig.legend(loc='lower right', bbox_to_anchor=(0.9, 0.25))
plt.show()
fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(18, 9))
ax[1, 2].remove()
fig.suptitle(setup_dict[setup_type])
offset = np.argmax(sim_data_results['HBM']['L2_X_coordinate_over_time']['time'] >= t_settle)
assert offset > 0
L2_x_displacement = convert_to_unit(sim_data_results['HBM']['L2_X_coordinate_over_time']['x-coordinate'], unit_deflection_plot, offset=offset) - convert_to_unit(sim_data_results['TESTBED']['Testbed_X_coordinate_over_time']['x-coordinate'], unit_deflection_plot, offset=offset)
L2_y_displacement = convert_to_unit(sim_data_results['HBM']['L2_Y_coordinate_over_time']['y-coordinate'], unit_deflection_plot, offset=offset) - convert_to_unit(sim_data_results['TESTBED']['Testbed_Y_coordinate_over_time']['y-coordinate'], unit_deflection_plot, offset=offset)
L2_z_displacement = convert_to_unit(sim_data_results['HBM']['L2_Z_coordinate_over_time']['z-coordinate'], unit_deflection_plot, offset=offset) - convert_to_unit(sim_data_results['TESTBED']['Testbed_Z_coordinate_over_time']['z-coordinate'], unit_deflection_plot, offset=offset)
ax[0, 0].plot(convert_to_unit(sim_data_results['HBM']['L2_X_coordinate_over_time']['time'], unit_time_plot, offset=offset), L2_x_displacement, label=HBM)
for data_exp in data_exp_L2_dx:
ax[0, 0].plot(data_exp[:, 0], -(data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[0, 1].plot(convert_to_unit(sim_data_results['HBM']['L2_Y_coordinate_over_time']['time'], unit_time_plot, offset=offset), L2_y_displacement)
for data_exp in data_exp_L2_dy:
ax[0, 1].plot(data_exp[:, 0], (data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[0, 2].plot(convert_to_unit(sim_data_results['HBM']['L2_Z_coordinate_over_time']['time'], unit_time_plot, offset=offset), L2_z_displacement)
for data_exp in data_exp_L2_dz:
ax[0, 2].plot(data_exp[:, 0], -(data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[1, 0].plot(L2_x_displacement, L2_z_displacement)
for data_exp in data_exp_L2_dx_dz:
ax[1, 0].plot(-(data_exp[:, 0] - data_exp[0, 0]), -(data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[1, 1].plot(L2_y_displacement, L2_z_displacement)
for data_exp in data_exp_L2_dy_dz:
ax[1, 1].plot((data_exp[:, 0] - data_exp[0, 0]), -(data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[0, 0].set(title="L2 X-Displacement", xlabel="Time [ms]", ylabel="Displacement [mm]")
ax[0, 1].set(title="L2 Y-Displacement", xlabel="Time [ms]", ylabel="Displacement [mm]")
ax[0, 2].set(title="L2 Z-Displacement", xlabel="Time [ms]", ylabel="Displacement [mm]")
ax[1, 0].set(title="L2 Trajectory", xlabel="X Displacement [mm]", ylabel="Z Displacement [mm]")
ax[1, 1].set(title="L2 Trajectory", xlabel="Y Displacement [mm]", ylabel="Z Displacement [mm]")
fig.tight_layout()
fig.legend(loc='lower right', bbox_to_anchor=(0.9, 0.25))
plt.show()
Pelvis displacement plots¶
In [ ]:
Copied!
fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(18, 9))
ax[1, 2].remove()
fig.suptitle(setup_dict[setup_type])
offset = np.argmax(sim_data_results['HBM']['Pelvis_X_coordinate_over_time']['time'] >= t_settle)
assert offset > 0
Pelvis_x_displacement = convert_to_unit(sim_data_results['HBM']['Pelvis_X_coordinate_over_time']['x-coordinate'], unit_deflection_plot, offset=offset) - convert_to_unit(sim_data_results['TESTBED']['Testbed_X_coordinate_over_time']['x-coordinate'], unit_deflection_plot, offset=offset)
Pelvis_y_displacement = convert_to_unit(sim_data_results['HBM']['Pelvis_Y_coordinate_over_time']['y-coordinate'], unit_deflection_plot, offset=offset) - convert_to_unit(sim_data_results['TESTBED']['Testbed_Y_coordinate_over_time']['y-coordinate'], unit_deflection_plot, offset=offset)
Pelvis_z_displacement = convert_to_unit(sim_data_results['HBM']['Pelvis_Z_coordinate_over_time']['z-coordinate'], unit_deflection_plot, offset=offset) - convert_to_unit(sim_data_results['TESTBED']['Testbed_Z_coordinate_over_time']['z-coordinate'], unit_deflection_plot, offset=offset)
time = convert_to_unit(sim_data_results['HBM']['Pelvis_X_coordinate_over_time']['time'], unit_time_plot, offset=offset)
sim_results_for_score['pelvis'] = [time, Pelvis_x_displacement, Pelvis_y_displacement, Pelvis_z_displacement]
ax[0, 0].plot(convert_to_unit(sim_data_results['HBM']['Pelvis_X_coordinate_over_time']['time'], unit_time_plot, offset=offset), Pelvis_x_displacement, label=HBM)
for data_exp in data_exp_Pelvis_dx:
ax[0, 0].plot(data_exp[:, 0], -(data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[0, 1].plot(convert_to_unit(sim_data_results['HBM']['Pelvis_Y_coordinate_over_time']['time'], unit_time_plot, offset=offset), Pelvis_y_displacement)
for data_exp in data_exp_Pelvis_dy:
ax[0, 1].plot(data_exp[:, 0], (data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[0, 2].plot(convert_to_unit(sim_data_results['HBM']['Pelvis_Z_coordinate_over_time']['time'], unit_time_plot, offset=offset), Pelvis_z_displacement)
for data_exp in data_exp_Pelvis_dz:
ax[0, 2].plot(data_exp[:, 0], -(data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[1, 0].plot(Pelvis_x_displacement, Pelvis_z_displacement)
for data_exp in data_exp_Pelvis_dx_dz:
ax[1, 0].plot(-(data_exp[:, 0] - data_exp[0, 0]), -(data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[1, 1].plot(Pelvis_y_displacement, Pelvis_z_displacement)
for data_exp in data_exp_Pelvis_dy_dz:
ax[1, 1].plot((data_exp[:, 0] - data_exp[0, 0]), -(data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[0, 0].set(title="Pelvis X-Displacement", xlabel="Time / ms", ylabel="Displacement / mm")
ax[0, 1].set(title="Pelvis Y-Displacement", xlabel="Time / ms", ylabel="Displacement / mm")
ax[0, 2].set(title="Pelvis Z-Displacement", xlabel="Time / ms", ylabel="Displacement / mm")
ax[1, 0].set(title="Pelvis Trajectory", xlabel="X Displacement / mm", ylabel="Z Displacement / mm")
ax[1, 1].set(title="Pelvis Trajectory", xlabel="Y Displacement / mm", ylabel="Z Displacement / mm")
fig.tight_layout()
fig.legend(loc='lower right', bbox_to_anchor=(0.9, 0.25))
plt.show()
fig, ax = plt.subplots(nrows=2, ncols=3, figsize=(18, 9))
ax[1, 2].remove()
fig.suptitle(setup_dict[setup_type])
offset = np.argmax(sim_data_results['HBM']['Pelvis_X_coordinate_over_time']['time'] >= t_settle)
assert offset > 0
Pelvis_x_displacement = convert_to_unit(sim_data_results['HBM']['Pelvis_X_coordinate_over_time']['x-coordinate'], unit_deflection_plot, offset=offset) - convert_to_unit(sim_data_results['TESTBED']['Testbed_X_coordinate_over_time']['x-coordinate'], unit_deflection_plot, offset=offset)
Pelvis_y_displacement = convert_to_unit(sim_data_results['HBM']['Pelvis_Y_coordinate_over_time']['y-coordinate'], unit_deflection_plot, offset=offset) - convert_to_unit(sim_data_results['TESTBED']['Testbed_Y_coordinate_over_time']['y-coordinate'], unit_deflection_plot, offset=offset)
Pelvis_z_displacement = convert_to_unit(sim_data_results['HBM']['Pelvis_Z_coordinate_over_time']['z-coordinate'], unit_deflection_plot, offset=offset) - convert_to_unit(sim_data_results['TESTBED']['Testbed_Z_coordinate_over_time']['z-coordinate'], unit_deflection_plot, offset=offset)
time = convert_to_unit(sim_data_results['HBM']['Pelvis_X_coordinate_over_time']['time'], unit_time_plot, offset=offset)
sim_results_for_score['pelvis'] = [time, Pelvis_x_displacement, Pelvis_y_displacement, Pelvis_z_displacement]
ax[0, 0].plot(convert_to_unit(sim_data_results['HBM']['Pelvis_X_coordinate_over_time']['time'], unit_time_plot, offset=offset), Pelvis_x_displacement, label=HBM)
for data_exp in data_exp_Pelvis_dx:
ax[0, 0].plot(data_exp[:, 0], -(data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[0, 1].plot(convert_to_unit(sim_data_results['HBM']['Pelvis_Y_coordinate_over_time']['time'], unit_time_plot, offset=offset), Pelvis_y_displacement)
for data_exp in data_exp_Pelvis_dy:
ax[0, 1].plot(data_exp[:, 0], (data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[0, 2].plot(convert_to_unit(sim_data_results['HBM']['Pelvis_Z_coordinate_over_time']['time'], unit_time_plot, offset=offset), Pelvis_z_displacement)
for data_exp in data_exp_Pelvis_dz:
ax[0, 2].plot(data_exp[:, 0], -(data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[1, 0].plot(Pelvis_x_displacement, Pelvis_z_displacement)
for data_exp in data_exp_Pelvis_dx_dz:
ax[1, 0].plot(-(data_exp[:, 0] - data_exp[0, 0]), -(data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[1, 1].plot(Pelvis_y_displacement, Pelvis_z_displacement)
for data_exp in data_exp_Pelvis_dy_dz:
ax[1, 1].plot((data_exp[:, 0] - data_exp[0, 0]), -(data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax[0, 0].set(title="Pelvis X-Displacement", xlabel="Time / ms", ylabel="Displacement / mm")
ax[0, 1].set(title="Pelvis Y-Displacement", xlabel="Time / ms", ylabel="Displacement / mm")
ax[0, 2].set(title="Pelvis Z-Displacement", xlabel="Time / ms", ylabel="Displacement / mm")
ax[1, 0].set(title="Pelvis Trajectory", xlabel="X Displacement / mm", ylabel="Z Displacement / mm")
ax[1, 1].set(title="Pelvis Trajectory", xlabel="Y Displacement / mm", ylabel="Z Displacement / mm")
fig.tight_layout()
fig.legend(loc='lower right', bbox_to_anchor=(0.9, 0.25))
plt.show()
Chest deflection plot¶
In [ ]:
Copied!
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(5, 6))
fig.suptitle(setup_dict[setup_type])
offset = np.argmax(sim_data_results['HBM']['Chest_Sternum_mid_X_coordinate_over_time']['time'] >= t_settle)
assert offset > 0
Time = convert_to_unit(sim_data_results['HBM']['Chest_defl_Sternum_mid_over_time']['time'], unit_time_plot, offset=offset)
Chest_defl = convert_to_unit(sim_data_results['HBM']['Chest_defl_Sternum_mid_over_time']['x-displacement'], unit_deflection_plot, offset=offset)
sim_results_for_score['chest deflection'] = [Time, -Chest_defl]
ax.plot(Time, -Chest_defl, label=HBM)
for data_exp in data_exp_chest_defl:
ax.plot(data_exp[:, 0], -(data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax.set(title="Chest deflection", xlabel="Time / ms", ylabel="Displacement / mm")
fig.tight_layout()
fig.legend(loc='center left', bbox_to_anchor=(1.0, 0.5))
plt.show()
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(5, 6))
fig.suptitle(setup_dict[setup_type])
offset = np.argmax(sim_data_results['HBM']['Chest_Sternum_mid_X_coordinate_over_time']['time'] >= t_settle)
assert offset > 0
Time = convert_to_unit(sim_data_results['HBM']['Chest_defl_Sternum_mid_over_time']['time'], unit_time_plot, offset=offset)
Chest_defl = convert_to_unit(sim_data_results['HBM']['Chest_defl_Sternum_mid_over_time']['x-displacement'], unit_deflection_plot, offset=offset)
sim_results_for_score['chest deflection'] = [Time, -Chest_defl]
ax.plot(Time, -Chest_defl, label=HBM)
for data_exp in data_exp_chest_defl:
ax.plot(data_exp[:, 0], -(data_exp[:, 1] - data_exp[0, 1]), **plotExperiment)
ax.set(title="Chest deflection", xlabel="Time / ms", ylabel="Displacement / mm")
fig.tight_layout()
fig.legend(loc='center left', bbox_to_anchor=(1.0, 0.5))
plt.show()
Seat force plots¶
In [ ]:
Copied!
fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(12, 5))
fig.suptitle(setup_dict[setup_type])
offset = np.argmax(sim_data_results['LOADCELL']['Loadcell_X_force_Seat_CFC60']['time'] >= t_settle)
assert offset > 0
for data_exp in data_exp_seat_fx:
ax1.plot(data_exp[:, 0], data_exp[:, 1], **plotExperiment)
if setup_type == '30kph_2kN':
time, force = subtract_no_occupant_data('LOADCELL', 'Loadcell_X_force_Seat_CFC60', 'x-force', offset)
ax1.plot(time, -force, label=settings["HBM"])
else:
ax1.plot(convert_to_unit(sim_data_results['LOADCELL']['Loadcell_X_force_Seat_CFC60']['time'], unit_time_plot, offset=offset),
-convert_to_unit(sim_data_results['LOADCELL']['Loadcell_X_force_Seat_CFC60']['x-force'], unit_force_plot, offset=offset),
label=HBM)
for data_exp in data_exp_seat_fz:
ax2.plot(data_exp[:, 0], data_exp[:, 1], **plotExperiment)
ax2.plot(convert_to_unit(sim_data_results['LOADCELL']['Loadcell_Z_force_Seat_CFC60']['time'], unit_time_plot, offset=offset),
-convert_to_unit(sim_data_results['LOADCELL']['Loadcell_Z_force_Seat_CFC60']['z-force'], unit_force_plot, offset=offset))
ax1.set(title="Seat X-Force", xlabel ="Time [ms]", ylabel="Force [kN]")
ax2.set(title="Seat Z-Force", xlabel ="Time [ms]", ylabel="Force [kN]")
ax1.set_xlim(right=250)
ax2.set_xlim(right=250)
fig.legend(loc='center left', bbox_to_anchor=(0.9, 0.5))
plt.show()
fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(12, 5))
fig.suptitle(setup_dict[setup_type])
offset = np.argmax(sim_data_results['LOADCELL']['Loadcell_X_force_Seat_CFC60']['time'] >= t_settle)
assert offset > 0
for data_exp in data_exp_seat_fx:
ax1.plot(data_exp[:, 0], data_exp[:, 1], **plotExperiment)
if setup_type == '30kph_2kN':
time, force = subtract_no_occupant_data('LOADCELL', 'Loadcell_X_force_Seat_CFC60', 'x-force', offset)
ax1.plot(time, -force, label=settings["HBM"])
else:
ax1.plot(convert_to_unit(sim_data_results['LOADCELL']['Loadcell_X_force_Seat_CFC60']['time'], unit_time_plot, offset=offset),
-convert_to_unit(sim_data_results['LOADCELL']['Loadcell_X_force_Seat_CFC60']['x-force'], unit_force_plot, offset=offset),
label=HBM)
for data_exp in data_exp_seat_fz:
ax2.plot(data_exp[:, 0], data_exp[:, 1], **plotExperiment)
ax2.plot(convert_to_unit(sim_data_results['LOADCELL']['Loadcell_Z_force_Seat_CFC60']['time'], unit_time_plot, offset=offset),
-convert_to_unit(sim_data_results['LOADCELL']['Loadcell_Z_force_Seat_CFC60']['z-force'], unit_force_plot, offset=offset))
ax1.set(title="Seat X-Force", xlabel ="Time [ms]", ylabel="Force [kN]")
ax2.set(title="Seat Z-Force", xlabel ="Time [ms]", ylabel="Force [kN]")
ax1.set_xlim(right=250)
ax2.set_xlim(right=250)
fig.legend(loc='center left', bbox_to_anchor=(0.9, 0.5))
plt.show()
Knee support force plots¶
In [ ]:
Copied!
fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(12, 5))
fig.suptitle(setup_dict[setup_type])
offset = np.argmax(sim_data_results['LOADCELL']['Loadcell_X_force_Kneesupport_left_CFC60']['time'] >= t_settle)
assert offset > 0
for data_exp in data_exp_knee_left_fx:
ax1.plot(data_exp[:, 0], data_exp[:, 1], **plotExperiment)
if setup_type == '30kph_2kN':
time, force = subtract_no_occupant_data('LOADCELL', 'Loadcell_X_force_Kneesupport_left_CFC60', 'x-force', offset)
ax1.plot(time, -force, label=HBM)
else:
ax1.plot(convert_to_unit(sim_data_results['LOADCELL']['Loadcell_X_force_Kneesupport_left_CFC60']['time'], unit_time_plot, offset=offset),
-convert_to_unit(sim_data_results['LOADCELL']['Loadcell_X_force_Kneesupport_left_CFC60']['x-force'], unit_force_plot, offset=offset),
label=HBM)
for data_exp in data_exp_knee_right_fx:
ax2.plot(data_exp[:, 0], data_exp[:, 1], **plotExperiment)
if setup_type == '30kph_2kN':
time, force = subtract_no_occupant_data('LOADCELL', 'Loadcell_X_force_Kneesupport_right_CFC60', 'x-force', offset)
ax2.plot(time, -force)
else:
ax2.plot(convert_to_unit(sim_data_results['LOADCELL']['Loadcell_X_force_Kneesupport_right_CFC60']['time'], unit_time_plot, offset=offset),
-convert_to_unit(sim_data_results['LOADCELL']['Loadcell_X_force_Kneesupport_right_CFC60']['x-force'], unit_force_plot, offset=offset))
ax1.set(title="Knee Support Left X-Force", xlabel ="Time [ms]", ylabel="Force [kN]")
ax2.set(title="Knee Support Right X-Force", xlabel ="Time [ms]", ylabel="Force [kN]")
ax1.set_xlim(right=250)
ax2.set_xlim(right=250)
fig.legend(loc='center left', bbox_to_anchor=(0.9, 0.5))
plt.show()
fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(12, 5))
fig.suptitle(setup_dict[setup_type])
offset = np.argmax(sim_data_results['LOADCELL']['Loadcell_X_force_Kneesupport_left_CFC60']['time'] >= t_settle)
assert offset > 0
for data_exp in data_exp_knee_left_fx:
ax1.plot(data_exp[:, 0], data_exp[:, 1], **plotExperiment)
if setup_type == '30kph_2kN':
time, force = subtract_no_occupant_data('LOADCELL', 'Loadcell_X_force_Kneesupport_left_CFC60', 'x-force', offset)
ax1.plot(time, -force, label=HBM)
else:
ax1.plot(convert_to_unit(sim_data_results['LOADCELL']['Loadcell_X_force_Kneesupport_left_CFC60']['time'], unit_time_plot, offset=offset),
-convert_to_unit(sim_data_results['LOADCELL']['Loadcell_X_force_Kneesupport_left_CFC60']['x-force'], unit_force_plot, offset=offset),
label=HBM)
for data_exp in data_exp_knee_right_fx:
ax2.plot(data_exp[:, 0], data_exp[:, 1], **plotExperiment)
if setup_type == '30kph_2kN':
time, force = subtract_no_occupant_data('LOADCELL', 'Loadcell_X_force_Kneesupport_right_CFC60', 'x-force', offset)
ax2.plot(time, -force)
else:
ax2.plot(convert_to_unit(sim_data_results['LOADCELL']['Loadcell_X_force_Kneesupport_right_CFC60']['time'], unit_time_plot, offset=offset),
-convert_to_unit(sim_data_results['LOADCELL']['Loadcell_X_force_Kneesupport_right_CFC60']['x-force'], unit_force_plot, offset=offset))
ax1.set(title="Knee Support Left X-Force", xlabel ="Time [ms]", ylabel="Force [kN]")
ax2.set(title="Knee Support Right X-Force", xlabel ="Time [ms]", ylabel="Force [kN]")
ax1.set_xlim(right=250)
ax2.set_xlim(right=250)
fig.legend(loc='center left', bbox_to_anchor=(0.9, 0.5))
plt.show()
Foot support force plots¶
In [ ]:
Copied!
fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(12, 5))
fig.suptitle(setup_dict[setup_type])
offset = np.argmax(sim_data_results['LOADCELL']['Loadcell_X_force_Footrest_left_CFC60']['time'] >= t_settle)
assert offset > 0
for data_exp in data_exp_feet_fx:
ax1.plot(data_exp[:, 0], data_exp[:, 1], **plotExperiment)
if setup_type == '30kph_2kN':
time_left, force_left = subtract_no_occupant_data('LOADCELL', 'Loadcell_X_force_Footrest_left_CFC60', 'x-force', offset)
time_right, force_right = subtract_no_occupant_data('LOADCELL', 'Loadcell_X_force_Footrest_right_CFC60', 'x-force', offset)
ax1.plot(time_left, -(force_left + force_right), label=organisation + " - " + HBM)
else:
ax1.plot(convert_to_unit(sim_data_results['LOADCELL']['Loadcell_X_force_Footrest_left_CFC60']['time'], unit_time_plot, offset=offset),
-(convert_to_unit(sim_data_results['LOADCELL']['Loadcell_X_force_Footrest_left_CFC60']['x-force'], unit_force_plot, offset=offset) +
convert_to_unit(sim_data_results['LOADCELL']['Loadcell_X_force_Footrest_right_CFC60']['x-force'], unit_force_plot, offset=offset)),
label=HBM)
for data_exp in data_exp_feet_fz:
ax2.plot(data_exp[:, 0], data_exp[:, 1], **plotExperiment)
ax1.set(title="Foot Support X-Force", xlabel ="Time [ms]", ylabel="Force [kN]")
ax2.set(title="Foot Support Z-Force", xlabel ="Time [ms]", ylabel="Force [kN]")
if setup_type == '30kph_2kN':
time_left, force_left = subtract_no_occupant_data('LOADCELL', 'Loadcell_Z_force_Footrest_left_CFC60', 'z-force', offset)
time_right, force_right = subtract_no_occupant_data('LOADCELL', 'Loadcell_Z_force_Footrest_right_CFC60', 'z-force', offset)
ax2.plot(time_left, -(force_left + force_right))
else:
ax2.plot(convert_to_unit(sim_data_results['LOADCELL']['Loadcell_Z_force_Footrest_left_CFC60']['time'], unit_time_plot, offset=offset),
-(convert_to_unit(sim_data_results['LOADCELL']['Loadcell_Z_force_Footrest_left_CFC60']['z-force'], unit_force_plot, offset=offset) +
convert_to_unit(sim_data_results['LOADCELL']['Loadcell_Z_force_Footrest_right_CFC60']['z-force'], unit_force_plot, offset=offset)))
ax1.set_xlim(right=250)
ax2.set_xlim(right=250)
fig.legend(loc='center left', bbox_to_anchor=(0.9, 0.5))
plt.show()
fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(12, 5))
fig.suptitle(setup_dict[setup_type])
offset = np.argmax(sim_data_results['LOADCELL']['Loadcell_X_force_Footrest_left_CFC60']['time'] >= t_settle)
assert offset > 0
for data_exp in data_exp_feet_fx:
ax1.plot(data_exp[:, 0], data_exp[:, 1], **plotExperiment)
if setup_type == '30kph_2kN':
time_left, force_left = subtract_no_occupant_data('LOADCELL', 'Loadcell_X_force_Footrest_left_CFC60', 'x-force', offset)
time_right, force_right = subtract_no_occupant_data('LOADCELL', 'Loadcell_X_force_Footrest_right_CFC60', 'x-force', offset)
ax1.plot(time_left, -(force_left + force_right), label=organisation + " - " + HBM)
else:
ax1.plot(convert_to_unit(sim_data_results['LOADCELL']['Loadcell_X_force_Footrest_left_CFC60']['time'], unit_time_plot, offset=offset),
-(convert_to_unit(sim_data_results['LOADCELL']['Loadcell_X_force_Footrest_left_CFC60']['x-force'], unit_force_plot, offset=offset) +
convert_to_unit(sim_data_results['LOADCELL']['Loadcell_X_force_Footrest_right_CFC60']['x-force'], unit_force_plot, offset=offset)),
label=HBM)
for data_exp in data_exp_feet_fz:
ax2.plot(data_exp[:, 0], data_exp[:, 1], **plotExperiment)
ax1.set(title="Foot Support X-Force", xlabel ="Time [ms]", ylabel="Force [kN]")
ax2.set(title="Foot Support Z-Force", xlabel ="Time [ms]", ylabel="Force [kN]")
if setup_type == '30kph_2kN':
time_left, force_left = subtract_no_occupant_data('LOADCELL', 'Loadcell_Z_force_Footrest_left_CFC60', 'z-force', offset)
time_right, force_right = subtract_no_occupant_data('LOADCELL', 'Loadcell_Z_force_Footrest_right_CFC60', 'z-force', offset)
ax2.plot(time_left, -(force_left + force_right))
else:
ax2.plot(convert_to_unit(sim_data_results['LOADCELL']['Loadcell_Z_force_Footrest_left_CFC60']['time'], unit_time_plot, offset=offset),
-(convert_to_unit(sim_data_results['LOADCELL']['Loadcell_Z_force_Footrest_left_CFC60']['z-force'], unit_force_plot, offset=offset) +
convert_to_unit(sim_data_results['LOADCELL']['Loadcell_Z_force_Footrest_right_CFC60']['z-force'], unit_force_plot, offset=offset)))
ax1.set_xlim(right=250)
ax2.set_xlim(right=250)
fig.legend(loc='center left', bbox_to_anchor=(0.9, 0.5))
plt.show()
Belt force plots¶
In [ ]:
Copied!
fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(12, 5))
fig.suptitle(setup_dict[setup_type])
offset = np.argmax(sim_data_results['TESTBED']['Shoulder_belt_res_force_over_time_CFC600']['time'] >= t_settle)
assert offset > 0
ax1.plot(convert_to_unit(sim_data_results['TESTBED']['Shoulder_belt_res_force_over_time_CFC600']['time'], unit_time_plot, offset=offset),
convert_to_unit(sim_data_results['TESTBED']['Shoulder_belt_res_force_over_time_CFC600']['res-force'], unit_force_plot, offset=offset),
label=HBM)
for data_exp in data_exp_shoulder_belt:
ax1.plot(data_exp[:, 0], data_exp[:, 1], **plotExperiment)
ax2.plot(convert_to_unit(sim_data_results['TESTBED']['Lap_belt_res_force_over_time_CFC600']['time'], unit_time_plot, offset=offset),
convert_to_unit(sim_data_results['TESTBED']['Lap_belt_res_force_over_time_CFC600']['res-force'], unit_force_plot, offset=offset))
for data_exp in data_exp_lap_belt:
ax2.plot(data_exp[:, 0], data_exp[:, 1], **plotExperiment)
ax1.set(title="Shoulder Belt Force", xlabel ="Time [ms]", ylabel="Force [kN]")
ax2.set(title="Lap Belt Force", xlabel ="Time [ms]", ylabel="Force [kN]")
ax1.set_xlim(right=250)
ax2.set_xlim(right=250)
fig.legend(loc='center left', bbox_to_anchor=(0.9, 0.5))
plt.show()
fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(12, 5))
fig.suptitle(setup_dict[setup_type])
offset = np.argmax(sim_data_results['TESTBED']['Shoulder_belt_res_force_over_time_CFC600']['time'] >= t_settle)
assert offset > 0
ax1.plot(convert_to_unit(sim_data_results['TESTBED']['Shoulder_belt_res_force_over_time_CFC600']['time'], unit_time_plot, offset=offset),
convert_to_unit(sim_data_results['TESTBED']['Shoulder_belt_res_force_over_time_CFC600']['res-force'], unit_force_plot, offset=offset),
label=HBM)
for data_exp in data_exp_shoulder_belt:
ax1.plot(data_exp[:, 0], data_exp[:, 1], **plotExperiment)
ax2.plot(convert_to_unit(sim_data_results['TESTBED']['Lap_belt_res_force_over_time_CFC600']['time'], unit_time_plot, offset=offset),
convert_to_unit(sim_data_results['TESTBED']['Lap_belt_res_force_over_time_CFC600']['res-force'], unit_force_plot, offset=offset))
for data_exp in data_exp_lap_belt:
ax2.plot(data_exp[:, 0], data_exp[:, 1], **plotExperiment)
ax1.set(title="Shoulder Belt Force", xlabel ="Time [ms]", ylabel="Force [kN]")
ax2.set(title="Lap Belt Force", xlabel ="Time [ms]", ylabel="Force [kN]")
ax1.set_xlim(right=250)
ax2.set_xlim(right=250)
fig.legend(loc='center left', bbox_to_anchor=(0.9, 0.5))
plt.show()
Trajectory plots using head, spine and pelvis markers¶
In [ ]:
Copied!
fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(11, 5))
fig.suptitle(setup_dict[setup_type])
offset = np.argmax(sim_data_results['HBM']['Head_CoG_Y_displacement_over_time']['time'] >= t_settle)
assert offset > 0
head_x_displacement = convert_to_unit(sim_data_results['HBM']['Head_CoG_X_coordinate_over_time']['x-coordinate'], unit_deflection_plot, offset=offset) + sim_data_results['HBM']['Head_CoG_X_coordinate_over_time']['x-coordinate'].iloc[0,0]
head_y_displacement = convert_to_unit(sim_data_results['HBM']['Head_CoG_Y_coordinate_over_time']['y-coordinate'], unit_deflection_plot, offset=offset) + sim_data_results['HBM']['Head_CoG_Y_coordinate_over_time']['y-coordinate'].iloc[0,0]
head_z_displacement = convert_to_unit(sim_data_results['HBM']['Head_CoG_Z_coordinate_over_time']['z-coordinate'], unit_deflection_plot, offset=offset) + sim_data_results['HBM']['Head_CoG_Z_coordinate_over_time']['z-coordinate'].iloc[0,0]
T1_x_displacement = convert_to_unit(sim_data_results['HBM']['T1_X_coordinate_over_time']['x-coordinate'], unit_deflection_plot, offset=offset) + sim_data_results['HBM']['T1_X_coordinate_over_time']['x-coordinate'].iloc[0,0]
T1_y_displacement = convert_to_unit(sim_data_results['HBM']['T1_Y_coordinate_over_time']['y-coordinate'], unit_deflection_plot, offset=offset) + sim_data_results['HBM']['T1_Y_coordinate_over_time']['y-coordinate'].iloc[0,0]
T1_z_displacement = convert_to_unit(sim_data_results['HBM']['T1_Z_coordinate_over_time']['z-coordinate'], unit_deflection_plot, offset=offset) + sim_data_results['HBM']['T1_Z_coordinate_over_time']['z-coordinate'].iloc[0,0]
T8_x_displacement = convert_to_unit(sim_data_results['HBM']['T8_X_coordinate_over_time']['x-coordinate'], unit_deflection_plot, offset=offset) + sim_data_results['HBM']['T8_X_coordinate_over_time']['x-coordinate'].iloc[0,0]
T8_y_displacement = convert_to_unit(sim_data_results['HBM']['T8_Y_coordinate_over_time']['y-coordinate'], unit_deflection_plot, offset=offset) + sim_data_results['HBM']['T8_Y_coordinate_over_time']['y-coordinate'].iloc[0,0]
T8_z_displacement = convert_to_unit(sim_data_results['HBM']['T8_Z_coordinate_over_time']['z-coordinate'], unit_deflection_plot, offset=offset) + sim_data_results['HBM']['T8_Z_coordinate_over_time']['z-coordinate'].iloc[0,0]
L2_x_displacement = convert_to_unit(sim_data_results['HBM']['L2_X_coordinate_over_time']['x-coordinate'], unit_deflection_plot, offset=offset) + sim_data_results['HBM']['L2_X_coordinate_over_time']['x-coordinate'].iloc[0,0]
L2_y_displacement = convert_to_unit(sim_data_results['HBM']['L2_Y_coordinate_over_time']['y-coordinate'], unit_deflection_plot, offset=offset) + sim_data_results['HBM']['L2_Y_coordinate_over_time']['y-coordinate'].iloc[0,0]
L2_z_displacement = convert_to_unit(sim_data_results['HBM']['L2_Z_coordinate_over_time']['z-coordinate'], unit_deflection_plot, offset=offset) + sim_data_results['HBM']['L2_Z_coordinate_over_time']['z-coordinate'].iloc[0,0]
Pelvis_x_displacement = convert_to_unit(sim_data_results['HBM']['Pelvis_X_coordinate_over_time']['x-coordinate'], unit_deflection_plot, offset=offset) + sim_data_results['HBM']['Pelvis_X_coordinate_over_time']['x-coordinate'].iloc[0,0]
Pelvis_y_displacement = convert_to_unit(sim_data_results['HBM']['Pelvis_Y_coordinate_over_time']['y-coordinate'], unit_deflection_plot, offset=offset) + sim_data_results['HBM']['Pelvis_Y_coordinate_over_time']['y-coordinate'].iloc[0,0]
Pelvis_z_displacement = convert_to_unit(sim_data_results['HBM']['Pelvis_Z_coordinate_over_time']['z-coordinate'], unit_deflection_plot, offset=offset) + sim_data_results['HBM']['Pelvis_Z_coordinate_over_time']['z-coordinate'].iloc[0,0]
head_x_displacement_rel = head_x_displacement - convert_to_unit(sim_data_results['TESTBED']['Testbed_X_coordinate_over_time']['x-coordinate'], unit_deflection_plot, offset=offset)
head_y_displacement_rel = head_y_displacement - convert_to_unit(sim_data_results['TESTBED']['Testbed_Y_coordinate_over_time']['y-coordinate'], unit_deflection_plot, offset=offset)
head_z_displacement_rel = head_z_displacement - convert_to_unit(sim_data_results['TESTBED']['Testbed_Z_coordinate_over_time']['z-coordinate'], unit_deflection_plot, offset=offset)
T1_x_displacement_rel = T1_x_displacement - convert_to_unit(sim_data_results['TESTBED']['Testbed_X_coordinate_over_time']['x-coordinate'], unit_deflection_plot, offset=offset)
T1_y_displacement_rel = T1_y_displacement - convert_to_unit(sim_data_results['TESTBED']['Testbed_Y_coordinate_over_time']['y-coordinate'], unit_deflection_plot, offset=offset)
T1_z_displacement_rel = T1_z_displacement - convert_to_unit(sim_data_results['TESTBED']['Testbed_Z_coordinate_over_time']['z-coordinate'], unit_deflection_plot, offset=offset)
T8_x_displacement_rel = T8_x_displacement - convert_to_unit(sim_data_results['TESTBED']['Testbed_X_coordinate_over_time']['x-coordinate'], unit_deflection_plot, offset=offset)
T8_y_displacement_rel = T8_y_displacement - convert_to_unit(sim_data_results['TESTBED']['Testbed_Y_coordinate_over_time']['y-coordinate'], unit_deflection_plot, offset=offset)
T8_z_displacement_rel = T8_z_displacement - convert_to_unit(sim_data_results['TESTBED']['Testbed_Z_coordinate_over_time']['z-coordinate'], unit_deflection_plot, offset=offset)
L2_x_displacement_rel = L2_x_displacement - convert_to_unit(sim_data_results['TESTBED']['Testbed_X_coordinate_over_time']['x-coordinate'], unit_deflection_plot, offset=offset)
L2_y_displacement_rel = L2_y_displacement - convert_to_unit(sim_data_results['TESTBED']['Testbed_Y_coordinate_over_time']['y-coordinate'], unit_deflection_plot, offset=offset)
L2_z_displacement_rel = L2_z_displacement - convert_to_unit(sim_data_results['TESTBED']['Testbed_Z_coordinate_over_time']['z-coordinate'], unit_deflection_plot, offset=offset)
Pelvis_x_displacement_rel = Pelvis_x_displacement - convert_to_unit(sim_data_results['TESTBED']['Testbed_X_coordinate_over_time']['x-coordinate'], unit_deflection_plot, offset=offset)
Pelvis_y_displacement_rel = Pelvis_y_displacement - convert_to_unit(sim_data_results['TESTBED']['Testbed_Y_coordinate_over_time']['y-coordinate'], unit_deflection_plot, offset=offset)
Pelvis_z_displacement_rel = Pelvis_z_displacement - convert_to_unit(sim_data_results['TESTBED']['Testbed_Z_coordinate_over_time']['z-coordinate'], unit_deflection_plot, offset=offset)
def get_indices(total_size, number_indices):
assert number_indices > 1
assert total_size >= number_indices
last_idx = total_size - 1
q, r = divmod(last_idx, number_indices - 1)
result = [0]
idx = 0
for i in range(1, number_indices):
idx += (q + 1) if (i <= r) else q
result.append(idx)
return result
for idx in get_indices(len(head_y_displacement_rel), 3):
ax2.plot([-head_x_displacement_rel[idx], -T1_x_displacement_rel[idx], -T8_x_displacement_rel[idx], -L2_x_displacement_rel[idx], -Pelvis_x_displacement_rel[idx]],
[head_z_displacement_rel[idx], T1_z_displacement_rel[idx], T8_z_displacement_rel[idx], L2_z_displacement_rel[idx], Pelvis_z_displacement_rel[idx]], marker='o', linestyle='--', color='grey', alpha=0.7)
ax1.plot(head_x_displacement, head_z_displacement, label="Head")
ax1.plot(T1_x_displacement, T1_z_displacement, label="T1")
ax1.plot(T8_x_displacement, T8_z_displacement, label="T8")
ax1.plot(L2_x_displacement, L2_z_displacement, label="L2")
ax1.plot(Pelvis_x_displacement, Pelvis_z_displacement, label="Pelvis")
ax2.plot(-head_x_displacement_rel, head_z_displacement_rel)
ax2.plot(-T1_x_displacement_rel, T1_z_displacement_rel)
ax2.plot(-T8_x_displacement_rel, T8_z_displacement_rel)
ax2.plot(-L2_x_displacement_rel, L2_z_displacement_rel)
ax2.plot(-Pelvis_x_displacement_rel, Pelvis_z_displacement_rel)
ax1.set(title="Absolute trajectories", xlabel="X-Displacement [mm]", ylabel="Z-Displacement [mm]")
ax2.set(title="Trajectories relative to sled", xlabel="X-Displacement [mm]", ylabel="Z-Displacement [mm]")
ax2.set_xlim(left=-250, right=500)
ax2.set_ylim(top=900, bottom=-100)
fig.legend(loc='center left', bbox_to_anchor=(0.9, 0.5))
plt.show()
fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(11, 5))
fig.suptitle(setup_dict[setup_type])
offset = np.argmax(sim_data_results['HBM']['Head_CoG_Y_displacement_over_time']['time'] >= t_settle)
assert offset > 0
head_x_displacement = convert_to_unit(sim_data_results['HBM']['Head_CoG_X_coordinate_over_time']['x-coordinate'], unit_deflection_plot, offset=offset) + sim_data_results['HBM']['Head_CoG_X_coordinate_over_time']['x-coordinate'].iloc[0,0]
head_y_displacement = convert_to_unit(sim_data_results['HBM']['Head_CoG_Y_coordinate_over_time']['y-coordinate'], unit_deflection_plot, offset=offset) + sim_data_results['HBM']['Head_CoG_Y_coordinate_over_time']['y-coordinate'].iloc[0,0]
head_z_displacement = convert_to_unit(sim_data_results['HBM']['Head_CoG_Z_coordinate_over_time']['z-coordinate'], unit_deflection_plot, offset=offset) + sim_data_results['HBM']['Head_CoG_Z_coordinate_over_time']['z-coordinate'].iloc[0,0]
T1_x_displacement = convert_to_unit(sim_data_results['HBM']['T1_X_coordinate_over_time']['x-coordinate'], unit_deflection_plot, offset=offset) + sim_data_results['HBM']['T1_X_coordinate_over_time']['x-coordinate'].iloc[0,0]
T1_y_displacement = convert_to_unit(sim_data_results['HBM']['T1_Y_coordinate_over_time']['y-coordinate'], unit_deflection_plot, offset=offset) + sim_data_results['HBM']['T1_Y_coordinate_over_time']['y-coordinate'].iloc[0,0]
T1_z_displacement = convert_to_unit(sim_data_results['HBM']['T1_Z_coordinate_over_time']['z-coordinate'], unit_deflection_plot, offset=offset) + sim_data_results['HBM']['T1_Z_coordinate_over_time']['z-coordinate'].iloc[0,0]
T8_x_displacement = convert_to_unit(sim_data_results['HBM']['T8_X_coordinate_over_time']['x-coordinate'], unit_deflection_plot, offset=offset) + sim_data_results['HBM']['T8_X_coordinate_over_time']['x-coordinate'].iloc[0,0]
T8_y_displacement = convert_to_unit(sim_data_results['HBM']['T8_Y_coordinate_over_time']['y-coordinate'], unit_deflection_plot, offset=offset) + sim_data_results['HBM']['T8_Y_coordinate_over_time']['y-coordinate'].iloc[0,0]
T8_z_displacement = convert_to_unit(sim_data_results['HBM']['T8_Z_coordinate_over_time']['z-coordinate'], unit_deflection_plot, offset=offset) + sim_data_results['HBM']['T8_Z_coordinate_over_time']['z-coordinate'].iloc[0,0]
L2_x_displacement = convert_to_unit(sim_data_results['HBM']['L2_X_coordinate_over_time']['x-coordinate'], unit_deflection_plot, offset=offset) + sim_data_results['HBM']['L2_X_coordinate_over_time']['x-coordinate'].iloc[0,0]
L2_y_displacement = convert_to_unit(sim_data_results['HBM']['L2_Y_coordinate_over_time']['y-coordinate'], unit_deflection_plot, offset=offset) + sim_data_results['HBM']['L2_Y_coordinate_over_time']['y-coordinate'].iloc[0,0]
L2_z_displacement = convert_to_unit(sim_data_results['HBM']['L2_Z_coordinate_over_time']['z-coordinate'], unit_deflection_plot, offset=offset) + sim_data_results['HBM']['L2_Z_coordinate_over_time']['z-coordinate'].iloc[0,0]
Pelvis_x_displacement = convert_to_unit(sim_data_results['HBM']['Pelvis_X_coordinate_over_time']['x-coordinate'], unit_deflection_plot, offset=offset) + sim_data_results['HBM']['Pelvis_X_coordinate_over_time']['x-coordinate'].iloc[0,0]
Pelvis_y_displacement = convert_to_unit(sim_data_results['HBM']['Pelvis_Y_coordinate_over_time']['y-coordinate'], unit_deflection_plot, offset=offset) + sim_data_results['HBM']['Pelvis_Y_coordinate_over_time']['y-coordinate'].iloc[0,0]
Pelvis_z_displacement = convert_to_unit(sim_data_results['HBM']['Pelvis_Z_coordinate_over_time']['z-coordinate'], unit_deflection_plot, offset=offset) + sim_data_results['HBM']['Pelvis_Z_coordinate_over_time']['z-coordinate'].iloc[0,0]
head_x_displacement_rel = head_x_displacement - convert_to_unit(sim_data_results['TESTBED']['Testbed_X_coordinate_over_time']['x-coordinate'], unit_deflection_plot, offset=offset)
head_y_displacement_rel = head_y_displacement - convert_to_unit(sim_data_results['TESTBED']['Testbed_Y_coordinate_over_time']['y-coordinate'], unit_deflection_plot, offset=offset)
head_z_displacement_rel = head_z_displacement - convert_to_unit(sim_data_results['TESTBED']['Testbed_Z_coordinate_over_time']['z-coordinate'], unit_deflection_plot, offset=offset)
T1_x_displacement_rel = T1_x_displacement - convert_to_unit(sim_data_results['TESTBED']['Testbed_X_coordinate_over_time']['x-coordinate'], unit_deflection_plot, offset=offset)
T1_y_displacement_rel = T1_y_displacement - convert_to_unit(sim_data_results['TESTBED']['Testbed_Y_coordinate_over_time']['y-coordinate'], unit_deflection_plot, offset=offset)
T1_z_displacement_rel = T1_z_displacement - convert_to_unit(sim_data_results['TESTBED']['Testbed_Z_coordinate_over_time']['z-coordinate'], unit_deflection_plot, offset=offset)
T8_x_displacement_rel = T8_x_displacement - convert_to_unit(sim_data_results['TESTBED']['Testbed_X_coordinate_over_time']['x-coordinate'], unit_deflection_plot, offset=offset)
T8_y_displacement_rel = T8_y_displacement - convert_to_unit(sim_data_results['TESTBED']['Testbed_Y_coordinate_over_time']['y-coordinate'], unit_deflection_plot, offset=offset)
T8_z_displacement_rel = T8_z_displacement - convert_to_unit(sim_data_results['TESTBED']['Testbed_Z_coordinate_over_time']['z-coordinate'], unit_deflection_plot, offset=offset)
L2_x_displacement_rel = L2_x_displacement - convert_to_unit(sim_data_results['TESTBED']['Testbed_X_coordinate_over_time']['x-coordinate'], unit_deflection_plot, offset=offset)
L2_y_displacement_rel = L2_y_displacement - convert_to_unit(sim_data_results['TESTBED']['Testbed_Y_coordinate_over_time']['y-coordinate'], unit_deflection_plot, offset=offset)
L2_z_displacement_rel = L2_z_displacement - convert_to_unit(sim_data_results['TESTBED']['Testbed_Z_coordinate_over_time']['z-coordinate'], unit_deflection_plot, offset=offset)
Pelvis_x_displacement_rel = Pelvis_x_displacement - convert_to_unit(sim_data_results['TESTBED']['Testbed_X_coordinate_over_time']['x-coordinate'], unit_deflection_plot, offset=offset)
Pelvis_y_displacement_rel = Pelvis_y_displacement - convert_to_unit(sim_data_results['TESTBED']['Testbed_Y_coordinate_over_time']['y-coordinate'], unit_deflection_plot, offset=offset)
Pelvis_z_displacement_rel = Pelvis_z_displacement - convert_to_unit(sim_data_results['TESTBED']['Testbed_Z_coordinate_over_time']['z-coordinate'], unit_deflection_plot, offset=offset)
def get_indices(total_size, number_indices):
assert number_indices > 1
assert total_size >= number_indices
last_idx = total_size - 1
q, r = divmod(last_idx, number_indices - 1)
result = [0]
idx = 0
for i in range(1, number_indices):
idx += (q + 1) if (i <= r) else q
result.append(idx)
return result
for idx in get_indices(len(head_y_displacement_rel), 3):
ax2.plot([-head_x_displacement_rel[idx], -T1_x_displacement_rel[idx], -T8_x_displacement_rel[idx], -L2_x_displacement_rel[idx], -Pelvis_x_displacement_rel[idx]],
[head_z_displacement_rel[idx], T1_z_displacement_rel[idx], T8_z_displacement_rel[idx], L2_z_displacement_rel[idx], Pelvis_z_displacement_rel[idx]], marker='o', linestyle='--', color='grey', alpha=0.7)
ax1.plot(head_x_displacement, head_z_displacement, label="Head")
ax1.plot(T1_x_displacement, T1_z_displacement, label="T1")
ax1.plot(T8_x_displacement, T8_z_displacement, label="T8")
ax1.plot(L2_x_displacement, L2_z_displacement, label="L2")
ax1.plot(Pelvis_x_displacement, Pelvis_z_displacement, label="Pelvis")
ax2.plot(-head_x_displacement_rel, head_z_displacement_rel)
ax2.plot(-T1_x_displacement_rel, T1_z_displacement_rel)
ax2.plot(-T8_x_displacement_rel, T8_z_displacement_rel)
ax2.plot(-L2_x_displacement_rel, L2_z_displacement_rel)
ax2.plot(-Pelvis_x_displacement_rel, Pelvis_z_displacement_rel)
ax1.set(title="Absolute trajectories", xlabel="X-Displacement [mm]", ylabel="Z-Displacement [mm]")
ax2.set(title="Trajectories relative to sled", xlabel="X-Displacement [mm]", ylabel="Z-Displacement [mm]")
ax2.set_xlim(left=-250, right=500)
ax2.set_ylim(top=900, bottom=-100)
fig.legend(loc='center left', bbox_to_anchor=(0.9, 0.5))
plt.show()
Rib criterion plots¶
In [ ]:
Copied!
# This cell prints the maximum rib fracture risk.
NFR1 = ("NFR=0", "NFR=1", "NFR=2", "NFR=3+")
NFR2 = ("NFR=0", "NFR=1+", "NFR=2+", "NFR=3+")
rib_risks = {
'25 yo': (np.round(np.squeeze(100*sim_data_results['RIBS']['Overall_Rib_Risk_25YO_' + rib_criterion]['risk'][0:4].values),1)),
'50 yo': (np.round(np.squeeze(100*sim_data_results['RIBS']['Overall_Rib_Risk_50YO_' + rib_criterion]['risk'][0:4].values),1)),
'75 yo': (np.round(np.squeeze(100*sim_data_results['RIBS']['Overall_Rib_Risk_75YO_' + rib_criterion]['risk'][0:4].values),1))
}
x = np.arange(len(NFR1)) # the label locations
width = 0.25 # the width of the bars
multiplier1 = 0
multiplier2 = 0
fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, constrained_layout=True, figsize=(10,5))
for attribute, measurement in rib_risks.items():
offset = width * multiplier1
rects = ax1.bar(x + offset, measurement, width, label=attribute)
ax1.bar_label(rects, padding=3)
multiplier1 += 1
ax1.set_ylabel('Risk / %')
ax1.set_title('Rib fracture risk max')
ax1.set_xticks(x + width)
ax1.set_xticklabels(NFR1)
ax1.legend(loc='best')
ax1.set_ylim(0, 110)
for attribute, measurement in rib_risks.items():
measurement[2] += measurement[3]
measurement[1] += measurement[2]
offset = width * multiplier2
rects = ax2.bar(x + offset, measurement, width, label=attribute)
ax2.bar_label(rects, padding=3)
multiplier2 += 1
ax2.set_ylabel('Risk / %')
ax2.set_title('Rib fracture risk max')
ax2.set_xticks(x + width)
ax2.set_xticklabels(NFR2)
ax2.legend(loc='best')
ax2.set_ylim(0, 110)
# This cell prints the maximum rib fracture risk.
NFR1 = ("NFR=0", "NFR=1", "NFR=2", "NFR=3+")
NFR2 = ("NFR=0", "NFR=1+", "NFR=2+", "NFR=3+")
rib_risks = {
'25 yo': (np.round(np.squeeze(100*sim_data_results['RIBS']['Overall_Rib_Risk_25YO_' + rib_criterion]['risk'][0:4].values),1)),
'50 yo': (np.round(np.squeeze(100*sim_data_results['RIBS']['Overall_Rib_Risk_50YO_' + rib_criterion]['risk'][0:4].values),1)),
'75 yo': (np.round(np.squeeze(100*sim_data_results['RIBS']['Overall_Rib_Risk_75YO_' + rib_criterion]['risk'][0:4].values),1))
}
x = np.arange(len(NFR1)) # the label locations
width = 0.25 # the width of the bars
multiplier1 = 0
multiplier2 = 0
fig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, constrained_layout=True, figsize=(10,5))
for attribute, measurement in rib_risks.items():
offset = width * multiplier1
rects = ax1.bar(x + offset, measurement, width, label=attribute)
ax1.bar_label(rects, padding=3)
multiplier1 += 1
ax1.set_ylabel('Risk / %')
ax1.set_title('Rib fracture risk max')
ax1.set_xticks(x + width)
ax1.set_xticklabels(NFR1)
ax1.legend(loc='best')
ax1.set_ylim(0, 110)
for attribute, measurement in rib_risks.items():
measurement[2] += measurement[3]
measurement[1] += measurement[2]
offset = width * multiplier2
rects = ax2.bar(x + offset, measurement, width, label=attribute)
ax2.bar_label(rects, padding=3)
multiplier2 += 1
ax2.set_ylabel('Risk / %')
ax2.set_title('Rib fracture risk max')
ax2.set_xticks(x + width)
ax2.set_xticklabels(NFR2)
ax2.legend(loc='best')
ax2.set_ylim(0, 110)
Energy plots¶
In [ ]:
Copied!
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(5,6))
fig.suptitle(setup_dict[setup_type] + ' - Model Energies')
ax.plot(convert_to_unit(sim_data_results.MODEL.Hourglass_Energy_time.time, unit_time_plot),
convert_to_unit(sim_data_results.MODEL.Hourglass_Energy_time.hourglass_energy, unit_energy_plot),
label = "Hourglass Energy")
ax.plot(convert_to_unit(sim_data_results.MODEL.Internal_Energy_time.time, unit_time_plot),
convert_to_unit(sim_data_results.MODEL.Internal_Energy_time.internal_energy, unit_energy_plot),
label = "Internal Energy")
ax.plot(convert_to_unit(sim_data_results.MODEL.Kinetic_Energy_time.time, unit_time_plot),
convert_to_unit(sim_data_results.MODEL.Kinetic_Energy_time.kinetic_energy, unit_energy_plot),
label = "Kinetic Energy")
ax.plot(convert_to_unit(sim_data_results.MODEL.Total_Energy_time.time, unit_time_plot),
convert_to_unit(sim_data_results.MODEL.Total_Energy_time.total_energy, unit_energy_plot),
label = "Total Energy")
ax.set(xlabel ="Time / ms", ylabel="Energy / J")
fig.legend(loc='center left', bbox_to_anchor=(1, 0.5))
fig.tight_layout(pad=0.5)
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(5,6))
fig.suptitle(setup_dict[setup_type] + ' - Model Energies')
ax.plot(convert_to_unit(sim_data_results.MODEL.Hourglass_Energy_time.time, unit_time_plot),
convert_to_unit(sim_data_results.MODEL.Hourglass_Energy_time.hourglass_energy, unit_energy_plot),
label = "Hourglass Energy")
ax.plot(convert_to_unit(sim_data_results.MODEL.Internal_Energy_time.time, unit_time_plot),
convert_to_unit(sim_data_results.MODEL.Internal_Energy_time.internal_energy, unit_energy_plot),
label = "Internal Energy")
ax.plot(convert_to_unit(sim_data_results.MODEL.Kinetic_Energy_time.time, unit_time_plot),
convert_to_unit(sim_data_results.MODEL.Kinetic_Energy_time.kinetic_energy, unit_energy_plot),
label = "Kinetic Energy")
ax.plot(convert_to_unit(sim_data_results.MODEL.Total_Energy_time.time, unit_time_plot),
convert_to_unit(sim_data_results.MODEL.Total_Energy_time.total_energy, unit_energy_plot),
label = "Total Energy")
ax.set(xlabel ="Time / ms", ylabel="Energy / J")
fig.legend(loc='center left', bbox_to_anchor=(1, 0.5))
fig.tight_layout(pad=0.5)
In [ ]:
Copied!
offset = np.argmax(sim_data_results.MODEL.Hourglass_Energy_time.time > t_settle) + 1
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(5,6))
fig.suptitle(setup_dict[setup_type] + ' - Model Energies')
ax.plot(convert_to_unit(sim_data_results.MODEL.Hourglass_Energy_time.time, unit_time_plot, offset=offset),
convert_to_unit(sim_data_results.MODEL.Hourglass_Energy_time.hourglass_energy, unit_energy_plot, offset=offset),
label = "Hourglass Energy")
ax.plot(convert_to_unit(sim_data_results.MODEL.Internal_Energy_time.time, unit_time_plot, offset=offset),
convert_to_unit(sim_data_results.MODEL.Internal_Energy_time.internal_energy, unit_energy_plot, offset=offset),
label = "Internal Energy")
ax.plot(convert_to_unit(sim_data_results.MODEL.Kinetic_Energy_time.time, unit_time_plot, offset=offset),
convert_to_unit(sim_data_results.MODEL.Kinetic_Energy_time.kinetic_energy, unit_energy_plot, offset=offset),
label = "Kinetic Energy")
ax.plot(convert_to_unit(sim_data_results.MODEL.Total_Energy_time.time, unit_time_plot, offset=offset),
convert_to_unit(sim_data_results.MODEL.Total_Energy_time.total_energy, unit_energy_plot, offset=offset),
label = "Total Energy")
ax.set(xlabel ="Time / ms", ylabel="Energy / J")
fig.legend(loc='center left', bbox_to_anchor=(1, 0.5))
fig.tight_layout(pad=0.5)
offset = np.argmax(sim_data_results.MODEL.Hourglass_Energy_time.time > t_settle) + 1
fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(5,6))
fig.suptitle(setup_dict[setup_type] + ' - Model Energies')
ax.plot(convert_to_unit(sim_data_results.MODEL.Hourglass_Energy_time.time, unit_time_plot, offset=offset),
convert_to_unit(sim_data_results.MODEL.Hourglass_Energy_time.hourglass_energy, unit_energy_plot, offset=offset),
label = "Hourglass Energy")
ax.plot(convert_to_unit(sim_data_results.MODEL.Internal_Energy_time.time, unit_time_plot, offset=offset),
convert_to_unit(sim_data_results.MODEL.Internal_Energy_time.internal_energy, unit_energy_plot, offset=offset),
label = "Internal Energy")
ax.plot(convert_to_unit(sim_data_results.MODEL.Kinetic_Energy_time.time, unit_time_plot, offset=offset),
convert_to_unit(sim_data_results.MODEL.Kinetic_Energy_time.kinetic_energy, unit_energy_plot, offset=offset),
label = "Kinetic Energy")
ax.plot(convert_to_unit(sim_data_results.MODEL.Total_Energy_time.time, unit_time_plot, offset=offset),
convert_to_unit(sim_data_results.MODEL.Total_Energy_time.total_energy, unit_energy_plot, offset=offset),
label = "Total Energy")
ax.set(xlabel ="Time / ms", ylabel="Energy / J")
fig.legend(loc='center left', bbox_to_anchor=(1, 0.5))
fig.tight_layout(pad=0.5)
Added mass of HBM¶
The added mass calculation will be available for version v1.0
In [ ]:
Copied!
from scipy.interpolate import interp1d
from arc_length_ellipse_score import ArcLengthEllipseScore
from scipy.interpolate import interp1d
from arc_length_ellipse_score import ArcLengthEllipseScore
In [ ]:
Copied!
def cut_at_time(data, time_max, interpolate=True, cut_negative_time=True):
# If interpolate is True, linear interpolate time step where cut is performed.
if cut_negative_time:
idx_time_geq_0 = np.argmax(data[:, 0] >= 0.0)
data = data[idx_time_geq_0:, :]
idx_time_larger = np.argmax(data[:, 0] > time_max)
if idx_time_larger == 0:
return data
else:
if interpolate:
assert np.shape(data)[1] == 2
assert data[idx_time_larger - 1, 0] <= time_max
assert data[idx_time_larger, 0] > time_max
if np.isclose(data[idx_time_larger - 1, 0], time_max):
return data[:idx_time_larger, :]
f_time_max = interp1d(data[idx_time_larger - 1: idx_time_larger + 1, 0], data[idx_time_larger - 1: idx_time_larger + 1, 1])(time_max)
return np.append(data[:idx_time_larger, :], [[time_max, f_time_max]], axis=0)
else:
return data[:idx_time_larger, :]
def remove_nan_rows(data):
return data[~np.isnan(data).any(axis=1), :]
def cut_at_time(data, time_max, interpolate=True, cut_negative_time=True):
# If interpolate is True, linear interpolate time step where cut is performed.
if cut_negative_time:
idx_time_geq_0 = np.argmax(data[:, 0] >= 0.0)
data = data[idx_time_geq_0:, :]
idx_time_larger = np.argmax(data[:, 0] > time_max)
if idx_time_larger == 0:
return data
else:
if interpolate:
assert np.shape(data)[1] == 2
assert data[idx_time_larger - 1, 0] <= time_max
assert data[idx_time_larger, 0] > time_max
if np.isclose(data[idx_time_larger - 1, 0], time_max):
return data[:idx_time_larger, :]
f_time_max = interp1d(data[idx_time_larger - 1: idx_time_larger + 1, 0], data[idx_time_larger - 1: idx_time_larger + 1, 1])(time_max)
return np.append(data[:idx_time_larger, :], [[time_max, f_time_max]], axis=0)
else:
return data[:idx_time_larger, :]
def remove_nan_rows(data):
return data[~np.isnan(data).any(axis=1), :]
In [ ]:
Copied!
if setup_type == '30kph_3kN':
cut_at_arc_length = False
num_warp_ctrl_points_arcgen = 1
data_exp_dict = {}
data_exp_dict['head'] = [data_exp_head_dx, data_exp_head_dy, data_exp_head_dz]
data_exp_dict['t1'] = [data_exp_T1_dx, data_exp_T1_dy, data_exp_T1_dz]
data_exp_dict['t8'] = [data_exp_T8_dx, data_exp_T8_dy, data_exp_T8_dz]
data_exp_dict['pelvis'] = [data_exp_Pelvis_dx, data_exp_Pelvis_dy, data_exp_Pelvis_dz]
data_exp_dict['chest deflection'] = [data_exp_chest_defl]
time_cut_ms_dict = {'head': 240, 't1': 240, 't8': 240, 'pelvis': 240, 'chest deflection': 150}
direction_dict = {'head': ['x', 'y', 'z'], 't1': ['x', 'y', 'z'], 't8': ['x', 'y', 'z'], 'pelvis': ['x', 'y', 'z'], 'chest deflection': [None]}
signs_dict = {'head': [-1, 1, -1], 't1': [-1, 1, -1], 't8': [-1, 1, -1], 'pelvis': [-1, 1, -1], 'chest deflection': [-1]}
for case in ['head', 't1', 't8', 'pelvis', 'chest deflection']:
for data_exp_, displacement, direction, sign in zip(data_exp_dict[case], sim_results_for_score[case][1:], direction_dict[case], signs_dict[case]):
if direction is None:
title = f'{case.capitalize()}'
else:
title = f'{case.capitalize()} {direction.upper()}-Displacement'
print(title)
time_cut_ms = time_cut_ms_dict[case]
# Prepare test data.
test_data_list = []
for test_data in data_exp_:
test_data = remove_nan_rows(test_data)
assert not np.any(np.isnan(test_data))
assert test_data[-1, 0] >= time_cut_ms
test_data = cut_at_time(test_data, time_cut_ms)
test_data = np.column_stack((test_data[:, 0], sign * (test_data[:, 1] - test_data[0, 1])))
assert test_data[0, 0] >= 0.0
test_data_list.append(test_data)
# Prepare simulation data.
sim_data_list = [cut_at_time(np.column_stack((sim_results_for_score[case][0], displacement)), time_cut_ms)]
# Calculate objective metric scores.
ales = ArcLengthEllipseScore(test_data_list, sim_data_list, cut_at_arc_length, num_warp_ctrl_points_arcgen, use_dtw=True)
ales.arc_length_score_for_simulation_data()
score = np.round(ales.scores_simulation_list[0], decimals=2)
score_dtw = np.round(ales.scores_dtw_simulation_list[0], decimals=2)
print(f'Objective metric score: {score}')
print(f'Objective metric score dtw: {score_dtw}')
# Plot arc length score over normalized arc length.
# ales.plot_score_over_arc_length()
label_simulation = f'{HBM}\nScore: {score}\nScore (dtw): {score_dtw}'
fig_list = ales.plot_data_and_results(title_fig=setup_dict[setup_type], indices_plot_ellipses=[150, 300, 499], title_ax=f'Objective Metric Score - {title}', xlabel='Time / ms', ylabel='Displacement / mm',
label_sim=label_simulation, figsize=(7,5), bbox_to_anchor=(1.1, 0.5), loc='center left', tight_layout=False, min_corridor_width=0.02)
else:
print('No Objective Metric Score was calculated, since data for AM50 anthropometry is available only for GS2 load case.')
if setup_type == '30kph_3kN':
cut_at_arc_length = False
num_warp_ctrl_points_arcgen = 1
data_exp_dict = {}
data_exp_dict['head'] = [data_exp_head_dx, data_exp_head_dy, data_exp_head_dz]
data_exp_dict['t1'] = [data_exp_T1_dx, data_exp_T1_dy, data_exp_T1_dz]
data_exp_dict['t8'] = [data_exp_T8_dx, data_exp_T8_dy, data_exp_T8_dz]
data_exp_dict['pelvis'] = [data_exp_Pelvis_dx, data_exp_Pelvis_dy, data_exp_Pelvis_dz]
data_exp_dict['chest deflection'] = [data_exp_chest_defl]
time_cut_ms_dict = {'head': 240, 't1': 240, 't8': 240, 'pelvis': 240, 'chest deflection': 150}
direction_dict = {'head': ['x', 'y', 'z'], 't1': ['x', 'y', 'z'], 't8': ['x', 'y', 'z'], 'pelvis': ['x', 'y', 'z'], 'chest deflection': [None]}
signs_dict = {'head': [-1, 1, -1], 't1': [-1, 1, -1], 't8': [-1, 1, -1], 'pelvis': [-1, 1, -1], 'chest deflection': [-1]}
for case in ['head', 't1', 't8', 'pelvis', 'chest deflection']:
for data_exp_, displacement, direction, sign in zip(data_exp_dict[case], sim_results_for_score[case][1:], direction_dict[case], signs_dict[case]):
if direction is None:
title = f'{case.capitalize()}'
else:
title = f'{case.capitalize()} {direction.upper()}-Displacement'
print(title)
time_cut_ms = time_cut_ms_dict[case]
# Prepare test data.
test_data_list = []
for test_data in data_exp_:
test_data = remove_nan_rows(test_data)
assert not np.any(np.isnan(test_data))
assert test_data[-1, 0] >= time_cut_ms
test_data = cut_at_time(test_data, time_cut_ms)
test_data = np.column_stack((test_data[:, 0], sign * (test_data[:, 1] - test_data[0, 1])))
assert test_data[0, 0] >= 0.0
test_data_list.append(test_data)
# Prepare simulation data.
sim_data_list = [cut_at_time(np.column_stack((sim_results_for_score[case][0], displacement)), time_cut_ms)]
# Calculate objective metric scores.
ales = ArcLengthEllipseScore(test_data_list, sim_data_list, cut_at_arc_length, num_warp_ctrl_points_arcgen, use_dtw=True)
ales.arc_length_score_for_simulation_data()
score = np.round(ales.scores_simulation_list[0], decimals=2)
score_dtw = np.round(ales.scores_dtw_simulation_list[0], decimals=2)
print(f'Objective metric score: {score}')
print(f'Objective metric score dtw: {score_dtw}')
# Plot arc length score over normalized arc length.
# ales.plot_score_over_arc_length()
label_simulation = f'{HBM}\nScore: {score}\nScore (dtw): {score_dtw}'
fig_list = ales.plot_data_and_results(title_fig=setup_dict[setup_type], indices_plot_ellipses=[150, 300, 499], title_ax=f'Objective Metric Score - {title}', xlabel='Time / ms', ylabel='Displacement / mm',
label_sim=label_simulation, figsize=(7,5), bbox_to_anchor=(1.1, 0.5), loc='center left', tight_layout=False, min_corridor_width=0.02)
else:
print('No Objective Metric Score was calculated, since data for AM50 anthropometry is available only for GS2 load case.')