Skip to content

Commit

Permalink
Bugfix in ntfa matrices, code cleanup
Browse files Browse the repository at this point in the history
  • Loading branch information
juliusgh committed Jul 16, 2023
1 parent 36212b0 commit c96a415
Show file tree
Hide file tree
Showing 2 changed files with 6 additions and 12 deletions.
4 changes: 2 additions & 2 deletions eg09_compute_ntfa_matrices.py
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@
from ntfa import read_snapshots, mode_identification, compute_tabular_data, save_tabular_data

np.random.seed(0)
for microstructure in microstructures[:]:
for microstructure in microstructures:
file_name, data_path, temp1, temp2, n_tests, sampling_alphas = itemgetter(
"file_name", "data_path", "temp1", "temp2", "n_tests", "sampling_alphas"
)(microstructure)
Expand Down Expand Up @@ -46,7 +46,7 @@
# Compare computed plastic modes with plastic modes from h5 file
plastic_modes = samples[0]['plastic_modes']
plastic_modes = plastic_modes_svd
# assert np.allclose(plastic_modes / np.sign(np.expand_dims(np.expand_dims(plastic_modes[0,0,:], axis=0), axis=0)), plastic_modes_svd), 'Identified plastic modes do not match plastic modes in h5 file'
assert np.allclose(plastic_modes / np.sign(np.expand_dims(np.expand_dims(plastic_modes[0,0,:], axis=0), axis=0)), plastic_modes_svd), 'Identified plastic modes do not match plastic modes in h5 file'

# Mode processing to compute system matrices

Expand Down
14 changes: 4 additions & 10 deletions ntfa.py
Original file line number Diff line number Diff line change
Expand Up @@ -206,9 +206,7 @@ def compute_tabular_data(samples, mesh, temperatures):
mat_id = mesh['mat_id']
n_gauss = mesh['n_gauss']
strain_dof = mesh['strain_dof']
global_gradient = mesh['global_gradient']
n_gp = mesh['n_integration_points']
n_phases = len(np.unique(mat_id))
n_modes = samples[0]['plastic_modes'].shape[-1]
n_temps = len(temperatures)
C_bar = np.zeros((strain_dof, strain_dof, n_temps))
Expand All @@ -221,11 +219,7 @@ def compute_tabular_data(samples, mesh, temperatures):
C1 = np.zeros((strain_dof, strain_dof, n_temps))
A0 = np.zeros((strain_dof, 7 + n_modes, n_temps))
A1 = np.zeros((strain_dof, 7 + n_modes, n_temps))
vol_frac0 = mesh['volume_fraction'][0]
vol_frac1 = mesh['volume_fraction'][1]
plastic_modes = samples[0]['plastic_modes']
# interpolate_temp = lambda x1, x2, alpha: x1 + alpha * (x2 - x1)
# dns_temperatures = interpolate_temp(temp1, temp2, sample_alphas)
sample_temperatures = np.array([sample['temperature'] for sample in samples])
temp1, temp2 = min(sample_temperatures), max(sample_temperatures)
sample_alphas = (sample_temperatures - temp1) / (temp2 - temp1)
Expand All @@ -238,7 +232,7 @@ def compute_tabular_data(samples, mesh, temperatures):
if np.floor(alpha) == alpha:
# sample for given temperature exists, no need for interpolation
id = upper_bound
C, eps = ref_C, ref_eps
approx_C, approx_eps = ref_C, ref_eps
E = samples[id]['strain_localization']
S = construct_stress_localization(E, ref_C, ref_eps, plastic_modes, mat_id, n_gauss, strain_dof)
else:
Expand All @@ -253,13 +247,13 @@ def compute_tabular_data(samples, mesh, temperatures):
sampling_eps = np.stack((samples[id0]['mat_thermal_strain'], samples[id1]['mat_thermal_strain'])).transpose([1, 0, 2, 3])

# interpolated quantities using an implicit interpolation scheme with four DOF
C, eps = opt4(sampling_C, sampling_eps, ref_C, ref_eps)
E, _ = interpolate_fluctuation_modes(E01, C, eps, plastic_modes, mat_id, n_gauss, strain_dof, n_modes, n_gp)
approx_C, approx_eps = opt4(sampling_C, sampling_eps, ref_C, ref_eps)
E, _ = interpolate_fluctuation_modes(E01, approx_C, approx_eps, plastic_modes, mat_id, n_gauss, strain_dof, n_modes, n_gp)
S = construct_stress_localization(E, ref_C, ref_eps, plastic_modes, mat_id, n_gauss, strain_dof)

# Compute NTFA matrices
C_bar[:, :, idx], tau_theta[:, idx], A_bar[:, :, idx], tau_xi[:, idx], D_xi[:, :, idx], D_theta[idx] = \
compute_ntfa_matrices(E, S, plastic_modes, eps, mesh)
compute_ntfa_matrices(E, S, plastic_modes, ref_eps, mesh)

# Compute phase average stresses
A0_full, A1_full = compute_phase_average_stress_localizations(E, ref_C, ref_eps, plastic_modes, mesh)
Expand Down

0 comments on commit c96a415

Please sign in to comment.