Errors: field-based position/angle locking

This example demonstrates how to create field errors on 3D geometries where we want to perturb sensor locations and orientation but we want our sensors to stay attached to the surfaces they are on. To demonstrate this we will use a thermo-mechanical simulation of a 10mm cube where we will place sensors on each of the 6 faces of the cube.

from pathlib import Path
import numpy as np
import matplotlib.pyplot as plt

# pyvale imports
import pyvale.mooseherder as mh
import pyvale.sensorsim as sens
import pyvale.dataset as dataset

1. Load physics simulation data

data_path: Path = dataset.element_case_output_path(dataset.EElemTest.HEX20)
sim_data: mh.SimData = mh.ExodusLoader(data_path).load_all_sim_data()

disp_keys = ("disp_x","disp_y","disp_z")
sim_data: mh.SimData = sens.scale_length_units(scale=1000.0,
                                               sim_data=sim_data,
                                               disp_keys=disp_keys)

2. Build virtual sensor arrays

# Simulation is a 10mm cube
sensor_positions = np.array(((5.0,0.0,5.0),     # cube x-z face
                             (5.0,10.0,5.0),    # cube x-z face
                             (5.0,5.0,0.0),     # cube x-y face
                             (5.0,5.0,10.0),    # cube x-y face
                             (0.0,5.0,5.0),     # cube y-z face
                             (10.0,5.0,5.0),))  # cube y-z face

sample_times = np.linspace(0.0,np.max(sim_data.time),50)

sens_data = sens.SensorData(positions=sensor_positions,
                            sample_times=sample_times)


sens_array: sens.SensorsPoint = sens.SensorFactory.vector_point(
    sim_data,
    sens_data,
    comp_keys=disp_keys,
    spatial_dims=sens.EDim.THREED,
    descriptor=sens.DescriptorFactory.displacement(),
)

2.1. Add simulated measurement errors

Our simulation is a cube and we have sensors on the midpoint of each face of the cube. We want to allow our sensors positions and orientation to be perturbed as part of our uncertainty simulation but we want the sensors to stay attached to the faces of the cube they are located on

pos_uncert = 1.0 # units = mm
pos_rand = (sens.GenUniform(low=-pos_uncert,high=pos_uncert),
            sens.GenUniform(low=-pos_uncert,high=pos_uncert),
            sens.GenUniform(low=-pos_uncert,high=pos_uncert))

pos_lock = np.full(sensor_positions.shape,False,dtype=bool)
pos_lock[0,1] = True # cube x-z face, lock y
pos_lock[1,1] = True # cube x-z face, lock y
pos_lock[2,2] = True # cube x-y face, lock z
pos_lock[3,2] = True # cube x-y face, lock z
pos_lock[4,0] = True # cube y-z face, lock x
pos_lock[5,0] = True # cube y-z face, lock x

ang_uncert = 2.0 # units = degrees
ang_rand = (sens.GenUniform(low=-ang_uncert,high=ang_uncert),
            sens.GenUniform(low=-ang_uncert,high=ang_uncert),
            sens.GenUniform(low=-ang_uncert,high=ang_uncert))

ang_lock = np.full(sensor_positions.shape,True,dtype=bool)
ang_lock[0,1] = False # cube x-z face, unlock/rotate about y
ang_lock[1,1] = False # cube x-z face, unlock/rotate about y
ang_lock[2,2] = False # cube x-y face, unlock/rotate about z
ang_lock[3,2] = False # cube x-y face, unlock/rotate about z
ang_lock[4,0] = False # cube y-z face, unlock/rotate about x
ang_lock[5,0] = False # cube y-z face, unlock/rotate about x

field_err_data = sens.ErrFieldData(pos_rand_xyz=pos_rand,
                                   pos_lock_xyz=pos_lock)

error_chain: list[sens.IErrSimulator] = [
    sens.ErrSysField(sens_array.get_field(),field_err_data),
    sens.ErrSysGenPercent(sens.GenUniform(low=-1.0,high=1.0),
                          sens.EErrDep.DEPENDENT),
    sens.ErrRandGenPercent(sens.GenNormal(std=1.0),sens.EErrDep.DEPENDENT),
]

sens_array.set_error_chain(error_chain)

3. Run a simulated experiment

measurements: np.ndarray = sens_array.sim_measurements()

truth: np.ndarray = sens_array.get_truth()
sys_errs: np.ndarray = sens_array.get_errors_systematic()
rand_errs: np.ndarray = sens_array.get_errors_random()

print(80*"-")
print("measurement = truth + sysematic error + random error")

print(f"measurements.shape = {measurements.shape} = "
        + "(n_sensors,n_field_components,n_timesteps)")
print(f"truth.shape     = {truth.shape}")
print(f"sys_errs.shape  = {sys_errs.shape}")
print(f"rand_errs.shape = {rand_errs.shape}")

sens_print: int = 2 # x-y face sensor
comp_print: int = 1 # component=disp_y
time_last: int = 5
time_print = slice(measurements.shape[2]-time_last,measurements.shape[2])

print(f"\nThese are the last {time_last} virtual measurements of sensor "
        + f"{sens_print}:\n")

sens.print_measurements(sens_array,sens_print,comp_print,time_print)
print("\n"+80*"-")

4. Analyse & visualise the results

output_path = Path.cwd() / "pyvale-output"
if not output_path.is_dir():
    output_path.mkdir(parents=True, exist_ok=True)

for kk in disp_keys:
    pv_plot = sens.plot_point_sensors_on_sim(sens_array,kk)
    pv_plot.camera_position = "yz"
    pv_plot.camera.azimuth = 45
    pv_plot.camera.elevation = 45

    # Set to False to show an interactive plot instead of saving the figure
    pv_plot.off_screen = True
    if pv_plot.off_screen:
        pv_plot.screenshot(output_path/f"ext_ex4d_locs_{kk}.png")
    else:
        pv_plot.show()

for kk in disp_keys:
    (fig,ax) = sens.plot_time_traces(sens_array,comp_key=kk)
    fig.savefig(output_path/f"ext_ex4d_traces_{kk}.png",
                dpi=300,
                bbox_inches="tight")

# Uncomment this to display the sensor trace plot
# plt.show()