Add support for plastic models

eg0_affine_thermoelastic_solver.py → eg00_affine_thermoelastic_solver.py

@@ -7,7 +7,7 @@
 from microstructures import *
 from utilities import verify_data, read_h5
 
-for microstructure in microstructures:
+for microstructure in microstructures[-9:]:
 
     file_name, data_path, temp1, temp2 = itemgetter('file_name', 'data_path', 'temp1', 'temp2')(microstructure)
 
@@ -17,6 +17,12 @@
 
     mesh, samples = read_h5(file_name, data_path, temperatures)
 
+    # print(mesh.keys())
+    # print(samples[0].keys())
+    print('strain localication shape:', samples[0]['strain_localization'].shape)
+    print('material stiffness shape:', samples[0]['mat_stiffness'].shape)
+    print('plastic modes shape:', samples[0]['plastic_modes'].shape)
+
     for sample in samples:
         verify_data(mesh, sample)
 

eg1_approximation_of_mat_properties.py → eg01_approximation_of_mat_properties.py

@@ -6,6 +6,8 @@
 from material_parameters import *
 from optimize_alpha import naive, opt1, opt2, opt4
 from utilities import plot_and_save, cm
+from matplotlib import rc
+rc('text', usetex=True)
 
 temp1 = 300
 temp2 = 1300

eg2_compare_approximations.py → eg02_compare_approximations.py

@@ -7,12 +7,15 @@
 
 from interpolate_fluctuation_modes import interpolate_fluctuation_modes
 from microstructures import *
+from material_parameters import *
 from optimize_alpha import opt1, opt2, opt4, naive
 from utilities import read_h5, construct_stress_localization, volume_average, compute_residual_efficient
+from matplotlib import rc
+rc('text', usetex=True) 
 
 np.random.seed(0)
 file_name, data_path, temp1, temp2, n_tests, sampling_alphas = itemgetter('file_name', 'data_path', 'temp1', 'temp2', 'n_tests',
-                                                                          'sampling_alphas')(microstructures[0])
+                                                                          'sampling_alphas')(microstructures[-8])
 print(file_name, '\t', data_path)
 
 n_loading_directions = 10
@@ -26,11 +29,11 @@
 strain_dof = mesh['strain_dof']
 global_gradient = mesh['global_gradient']
 n_gp = mesh['n_integration_points']
-n_modes = ref[0]['strain_localization'].shape[-1]
+n_modes = ref[0]['plastic_modes'].shape[-1]
 
 _, samples = read_h5(file_name, data_path, [temp1, temp2], get_mesh=False)
 
-strains = np.random.normal(size=(n_loading_directions, mesh['strain_dof']))
+strains = np.random.normal(size=(n_loading_directions, strain_dof))
 strains /= la.norm(strains, axis=1)[:, None]
 
 n_approaches = 5
@@ -54,38 +57,46 @@
     Eref = ref[idx]['strain_localization']
     ref_C = ref[idx]['mat_stiffness']
     ref_eps = ref[idx]['mat_thermal_strain']
+    ref_C_ = np.stack(([stiffness_cu(temperature), stiffness_wsc(temperature)]))
+    ref_eps_ = np.expand_dims(np.stack(([thermal_strain_cu(temperature), thermal_strain_wsc(temperature)])), axis=2)
+    print(np.linalg.norm(ref_C - ref_C_), np.linalg.norm(ref_eps - ref_eps_))
+    plastic_modes = ref[idx]['plastic_modes']
     normalization_factor_mech = ref[idx]['normalization_factor_mech']
 
-    Sref = construct_stress_localization(Eref, ref_C, ref_eps, mat_id, n_gauss, strain_dof)
+    Sref = construct_stress_localization(Eref, ref_C, ref_eps, plastic_modes, mat_id, n_gauss, strain_dof)
     effSref = volume_average(Sref)
 
     # interpolated quantities using an explicit interpolation scheme with one DOF
     approx_C, approx_eps = naive(alpha, sampling_C, sampling_eps, ref_C, ref_eps)
     Enaive = interpolate_temp(E0, E1)
-    Snaive = construct_stress_localization(Enaive, ref_C, ref_eps, mat_id, n_gauss, strain_dof)
+    Snaive = construct_stress_localization(Enaive, ref_C, ref_eps, plastic_modes, mat_id, n_gauss, strain_dof)
     effSnaive = volume_average(Snaive)
 
     # interpolated quantities using an explicit interpolation scheme with one DOF
-    Eopt0, _ = interpolate_fluctuation_modes(E01, approx_C, approx_eps, mat_id, n_gauss, strain_dof, n_modes, n_gp)
-    Sopt0 = construct_stress_localization(Eopt0, ref_C, ref_eps, mat_id, n_gauss, strain_dof)
+    Eopt0, _ = interpolate_fluctuation_modes(E01, approx_C, approx_eps, plastic_modes, mat_id, n_gauss, strain_dof, n_modes,
+                                             n_gp)
+    Sopt0 = construct_stress_localization(Eopt0, ref_C, ref_eps, plastic_modes, mat_id, n_gauss, strain_dof)
     effSopt0 = volume_average(Sopt0)
 
     # interpolated quantities using an implicit interpolation scheme with one DOF
     approx_C, approx_eps = opt1(sampling_C, sampling_eps, ref_C, ref_eps)
-    Eopt1, _ = interpolate_fluctuation_modes(E01, approx_C, approx_eps, mat_id, n_gauss, strain_dof, n_modes, n_gp)
-    Sopt1 = construct_stress_localization(Eopt1, ref_C, ref_eps, mat_id, n_gauss, strain_dof)
+    Eopt1, _ = interpolate_fluctuation_modes(E01, approx_C, approx_eps, plastic_modes, mat_id, n_gauss, strain_dof, n_modes,
+                                             n_gp)
+    Sopt1 = construct_stress_localization(Eopt1, ref_C, ref_eps, plastic_modes, mat_id, n_gauss, strain_dof)
     effSopt1 = volume_average(Sopt1)
 
     # interpolated quantities using an implicit interpolation scheme with two DOF
     approx_C, approx_eps = opt2(sampling_C, sampling_eps, ref_C, ref_eps)
-    Eopt2, _ = interpolate_fluctuation_modes(E01, approx_C, approx_eps, mat_id, n_gauss, strain_dof, n_modes, n_gp)
-    Sopt2 = construct_stress_localization(Eopt2, ref_C, ref_eps, mat_id, n_gauss, strain_dof)
+    Eopt2, _ = interpolate_fluctuation_modes(E01, approx_C, approx_eps, plastic_modes, mat_id, n_gauss, strain_dof, n_modes,
+                                             n_gp)
+    Sopt2 = construct_stress_localization(Eopt2, ref_C, ref_eps, plastic_modes, mat_id, n_gauss, strain_dof)
     effSopt2 = volume_average(Sopt2)
 
     # interpolated quantities using an implicit interpolation scheme with four DOF
     approx_C, approx_eps = opt4(sampling_C, sampling_eps, ref_C, ref_eps)
-    Eopt4, _ = interpolate_fluctuation_modes(E01, approx_C, approx_eps, mat_id, n_gauss, strain_dof, n_modes, n_gp)
-    Sopt4 = construct_stress_localization(Eopt4, ref_C, ref_eps, mat_id, n_gauss, strain_dof)
+    Eopt4, _ = interpolate_fluctuation_modes(E01, approx_C, approx_eps, plastic_modes, mat_id, n_gauss, strain_dof, n_modes,
+                                             n_gp)
+    Sopt4 = construct_stress_localization(Eopt4, ref_C, ref_eps, plastic_modes, mat_id, n_gauss, strain_dof)
     effSopt4 = volume_average(Sopt4)
 
     err = lambda x, y: np.mean(la.norm(x - y, axis=(-1, -2)) / la.norm(y, axis=(-1, -2))) * 100
@@ -97,7 +108,7 @@
                          err(effSopt2, effSref), err(effSopt4, effSref)]
 
     for strain_idx, strain in enumerate(strains):
-        zeta = np.hstack((strain, 1))
+        zeta = np.hstack((strain, 1, np.ones(plastic_modes.shape[-1])))
 
         eff_stress_ref = effSref @ zeta
         err_eff_stress[:, idx * n_loading_directions + strain_idx] = \
@@ -112,7 +123,7 @@
         stress_opt4 = np.einsum('ijk,k', Sopt4, zeta, optimize='optimal')
 
         residuals = compute_residual_efficient([stress_naive, stress_opt0, stress_opt1, stress_opt2, stress_opt4],
-                                               mesh['global_gradient'])
+                                               global_gradient)
 
         err_f[:, idx * n_loading_directions + strain_idx] = la.norm(residuals, np.inf, axis=0) / normalization_factor_mech * 100
 

eg3_hierarchical_sampling.py → eg03_hierarchical_sampling.py

@@ -10,14 +10,16 @@
 from optimize_alpha import opt4
 from utilities import read_h5, construct_stress_localization, volume_average, compute_residual_efficient, \
     compute_err_indicator_efficient
+from matplotlib import rc
+rc('text', usetex=True) 
 
 np.random.seed(0)
 file_name, data_path, temp1, temp2, n_tests, sampling_alphas = itemgetter('file_name', 'data_path', 'temp1', 'temp2', 'n_tests',
-                                                                          'sampling_alphas')(microstructures[6])
+                                                                          'sampling_alphas')(microstructures[-8])
 print(file_name, '\t', data_path)
 
 n_loading_directions = 1
-n_hierarchical_levels = 5
+n_hierarchical_levels = 2
 test_temperatures = np.linspace(temp1, temp2, num=n_tests)
 test_alphas = np.linspace(0, 1, num=n_tests)
 
@@ -28,7 +30,7 @@
 global_gradient = mesh['global_gradient']
 n_gp = mesh['n_integration_points']
 n_phases = len(np.unique(mat_id))
-n_modes = ref[0]['strain_localization'].shape[-1]
+n_modes = ref[0]['plastic_modes'].shape[-1]
 
 strains = np.random.normal(size=(n_loading_directions, strain_dof))
 strains /= la.norm(strains, axis=1)[:, None]
@@ -60,22 +62,27 @@
         sampling_eps = np.stack((samples[id0]['mat_thermal_strain'], samples[id1]['mat_thermal_strain'])).transpose([1, 0, 2, 3])
 
         # reference values
+        Eref = ref[idx]['strain_localization']
         ref_C = ref[idx]['mat_stiffness']
         ref_eps = ref[idx]['mat_thermal_strain']
+        plastic_modes = ref[idx]['plastic_modes']
         normalization_factor_mech = ref[idx]['normalization_factor_mech']
-        effSref = np.vstack((ref[idx]['eff_stiffness'], -ref[idx]['eff_stiffness'] @ ref[idx]['eff_thermal_strain'])).T
+        Sref = construct_stress_localization(Eref, ref_C, ref_eps, plastic_modes, mat_id, n_gauss, strain_dof)
+        effSref = volume_average(Sref)
+        # effSref = np.vstack((ref[idx]['eff_stiffness'], -ref[idx]['eff_stiffness'] @ ref[idx]['eff_thermal_strain'])).T
 
         # interpolated quantities using an implicit interpolation scheme with four DOF
         approx_C, approx_eps = opt4(sampling_C, sampling_eps, ref_C, ref_eps)
-        Eopt4, _ = interpolate_fluctuation_modes(E01, approx_C, approx_eps, mat_id, n_gauss, strain_dof, n_modes, n_gp)
-        Sopt4 = construct_stress_localization(Eopt4, ref_C, ref_eps, mat_id, n_gauss, strain_dof)
+        Eopt4, _ = interpolate_fluctuation_modes(E01, approx_C, approx_eps, plastic_modes, mat_id, n_gauss, strain_dof, n_modes,
+                                                 n_gp)
+        Sopt4 = construct_stress_localization(Eopt4, ref_C, ref_eps, plastic_modes, mat_id, n_gauss, strain_dof)
         effSopt = volume_average(Sopt4)
 
         err_indicators[level, idx] = np.mean(np.max(np.abs(compute_err_indicator_efficient(Sopt4, global_gradient)),
                                                     axis=0)) / normalization_factor_mech * 100
 
         for strain_idx, strain in enumerate(strains):
-            zeta = np.hstack((strain, 1))
+            zeta = np.hstack((strain, 1, np.ones(plastic_modes.shape[-1])))
             stress_opt4 = np.einsum('ijk,k', Sopt4, zeta, optimize='optimal')
             residual = compute_residual_efficient(stress_opt4, global_gradient)
 
@@ -89,11 +96,11 @@
         invL = la.inv(la.cholesky(Cref))
         err_eff_C[level, idx] = la.norm(invL @ Capprox @ invL.T - np.eye(6)) / la.norm(np.eye(6)) * 100
 
-        err_eff_eps[level, idx] = err(la.solve(Capprox, effSopt[:, -1]), la.solve(Cref, effSref[:, -1]))
+        err_eff_eps[level, idx] = err(la.solve(Capprox, effSopt[:, 7]), la.solve(Cref, effSref[:, 7]))
 
     # max_err_idx = np.argmax(np.mean(err_nodal_force[level], axis=1))
     max_err_idx = np.argmax(err_indicators[level])
-    alpha_levels.append(np.sort(np.hstack((alphas, test_alphas[max_err_idx]))))
+    alpha_levels.append(np.unique(np.sort(np.hstack((alphas, test_alphas[max_err_idx])))))
     print(f'{np.max(np.mean(err_nodal_force[level], axis=1)) = }')
     print(f'{np.max(err_indicators[level]) = }')
 

eg4_hierarchical_sampling_efficient.py → eg04_hierarchical_sampling_efficient.py

@@ -6,13 +6,13 @@
 from interpolate_fluctuation_modes import update_affine_decomposition, effective_S, effective_stress_localization, \
     interpolate_fluctuation_modes, get_phi, transform_strain_localization
 from microstructures import *
-from optimize_alpha import opt4_alphas
-from utilities import read_h5, construct_stress_localization, compute_err_indicator_efficient
+from optimize_alpha import opt4_alphas, opt4
+from utilities import read_h5, construct_stress_localization, compute_err_indicator_efficient, volume_average
 
 np.random.seed(0)
 # np.set_printoptions(precision=3)
 
-for ms_id in [6, 7, 8, 9]:
+for ms_id in [0]:
     file_name, data_path, temp1, temp2, n_tests, sampling_alphas = itemgetter('file_name', 'data_path', 'temp1', 'temp2',
                                                                               'n_tests',
                                                                               'sampling_alphas')(microstructures[ms_id])
@@ -36,20 +36,25 @@
     global_gradient = mesh['global_gradient']
     n_gp = mesh['n_integration_points']
     n_phases = len(np.unique(mat_id))
-    n_modes = refs[0]['strain_localization'].shape[-1]
+    N_modes = refs[0]['strain_localization'].shape[2] - 7
 
     # extract temperature dependent data from the reference solutions
     # such as: material stiffness and thermal strain at each temperature and for all phases
+    Erefs = np.zeros((n_tests, *refs[0]['strain_localization'].shape))  # n_tests x n_phases x 6 x 6
     ref_Cs = np.zeros((n_tests, *refs[0]['mat_stiffness'].shape))  # n_tests x n_phases x 6 x 6
     ref_epss = np.zeros((n_tests, *refs[0]['mat_thermal_strain'].shape))  # n_tests x n_phases x 6 x 1
-    effSref = np.zeros((n_tests, strain_dof, n_modes))
+    effSref = np.zeros((n_tests, strain_dof, N_modes + 7))  # n_tests x 6 x (N + 7)
     normalization_factor_mech = np.zeros((n_tests))
+    plastic_modes = refs[0]['plastic_modes']  # temperature independent
     for idx, alpha in enumerate(test_alphas):
+        print(idx)
+        Erefs[idx] = refs[idx]['strain_localization']
         ref_Cs[idx] = refs[idx]['mat_stiffness']
         ref_epss[idx] = refs[idx]['mat_thermal_strain']
         normalization_factor_mech[idx] = refs[idx]['normalization_factor_mech']
-        effSref[idx] = np.hstack(
-            (refs[idx]['eff_stiffness'], -np.reshape(refs[idx]['eff_stiffness'] @ refs[idx]['eff_thermal_strain'], (-1, 1))))
+        Sref = construct_stress_localization(Erefs[idx], ref_Cs[idx], ref_epss[idx], plastic_modes, mat_id, n_gauss,
+                                             strain_dof)
+        effSref[idx] = volume_average(Sref)
 
     err_indicators, err_eff_S, err_eff_C, err_eff_eps = [np.zeros((n_hierarchical_levels, n_tests)) for _ in range(4)]
     interpolate_temp = lambda x1, x2, alpha: x1 + alpha * (x2 - x1)
@@ -105,30 +110,39 @@
                 current_sampling_id = alphas_indexing[idx]
 
                 K0, K1, F0, F1, F2, F3, S001, S101, S103, S002, S102, S104 = update_affine_decomposition(
-                    E01s[current_sampling_id], sampling_C, sampling_eps, n_modes, n_phases, n_gp, strain_dof, mat_id, n_gauss)
+                    E01s[current_sampling_id], sampling_C, sampling_eps, plastic_modes, N_modes, n_phases,
+                    n_gp, strain_dof,
+                    mat_id, n_gauss)
 
             phi = get_phi(K0, K1, F0, F1, F2, F3, alpha_C, alpha_eps, alpha_C_eps)
 
             speed = 1
-            if speed == 0:
-                C, eps = ref_Cs[idx], ref_epss[idx]
-                # C, eps = opt4(sampling_C, sampling_eps, ref_Cs[idx], ref_epss[idx])
-                _, effSopt = interpolate_fluctuation_modes(E01s[current_sampling_id], C, eps, mat_id, n_gauss, strain_dof,
-                                                           n_modes, n_gp)
-            elif speed == 1:
-                effSopt = effective_stress_localization(E01s[current_sampling_id], phi, ref_Cs[idx], ref_epss[idx], mat_id,
-                                                        n_gauss, n_gp, strain_dof, n_modes)
-            elif speed == 2:
-                # matches the result from interpolate_fluctuation_modes with a difference
-                # that depends on using ref_Cs[idx],ref_epss[idx] instead of alphas
-                effSopt, phi = effective_S(phi, S001, S101, S103, S002, S102, S104, alpha_C, np.squeeze(alpha_eps, axis=-1),
-                                           np.squeeze(alpha_C_eps, axis=-1))
-            else:
-                raise NotImplementedError()
+            # if speed == 0:
+            C, eps = ref_Cs[idx], ref_epss[idx]
+            # C, eps = opt4(sampling_C, sampling_eps, ref_Cs[idx], ref_epss[idx])
+            _, effSopt = interpolate_fluctuation_modes(E01s[current_sampling_id], C, eps, plastic_modes, mat_id, n_gauss,
+                                                        strain_dof,
+                                                        N_modes, n_gp)
+            #elif speed == 1:
+            # TODO verify the result when plasticity is on
+            effSopt1 = effective_stress_localization(E01s[current_sampling_id], phi, ref_Cs[idx], ref_epss[idx], plastic_modes,
+                                                    mat_id,
+                                                    n_gauss, n_gp, strain_dof, N_modes)
+            #elif speed == 2:
+            # TODO verify the result when plasticity is on
+            # matches the result from interpolate_fluctuatioN_modes with a difference
+            # that depends on using ref_Cs[idx],ref_epss[idx] instead of alphas
+            effSopt2, phi2 = effective_S(phi, S001, S101, S103, S002, S102, S104, alpha_C, np.squeeze(alpha_eps, axis=-1),
+                                        np.squeeze(alpha_C_eps, axis=-1))
+            #else:
+            #    raise NotImplementedError()
+            print(np.linalg.norm(effSopt - effSopt1))
 
             if not given_alpha_levels:
-                Eopt4 = transform_strain_localization(E01s[current_sampling_id], phi, n_gp, strain_dof, n_modes)
-                Sopt4 = construct_stress_localization(Eopt4, ref_Cs[idx], ref_epss[idx], mat_id, n_gauss, strain_dof)
+                Eopt4 = transform_strain_localization(E01s[current_sampling_id], phi, n_gp, strain_dof, N_modes)
+                Sopt4 = construct_stress_localization(Eopt4, ref_Cs[idx], ref_epss[idx], plastic_modes, mat_id, n_gauss,
+                                                      strain_dof)
+                # effSopt = volume_average(Sopt4)
                 err_indicators[level,
                                idx] = np.mean(np.max(np.abs(compute_err_indicator_efficient(Sopt4, global_gradient)),
                                                      axis=0)) / normalization_factor_mech[idx] * 100
@@ -140,7 +154,7 @@
             invL = la.inv(la.cholesky(Cref))
 
             err_eff_C[level, idx] = la.norm(invL @ Capprox @ invL.T - np.eye(6)) / la.norm(np.eye(6)) * 100
-            err_eff_eps[level, idx] = err(la.solve(Capprox, effSopt[:, -1]), la.solve(Cref, effSref[idx][:, -1]))
+            err_eff_eps[level, idx] = err(la.solve(Capprox, effSopt[:, 7]), la.solve(Cref, effSref[idx][:, 7]))
 
             # TODO remove dtype='f'
             group = file.require_group(f'{data_path}_level{level}')
@@ -149,7 +163,7 @@
             dset_stiffness = group.require_dataset(f'eff_stiffness_{temperature:07.2f}', (6, 6), dtype='f')
             dset_thermal_strain = group.require_dataset(f'eff_thermal_strain_{temperature:07.2f}', (6), dtype='f')
             dset_stiffness[:] = Capprox.T
-            dset_thermal_strain[:] = la.solve(Capprox, effSopt[:, -1])
+            dset_thermal_strain[:] = la.solve(Capprox, effSopt[:, 7])
 
         if not given_alpha_levels:
             max_err_idx = np.argmax(err_indicators[level])

eg8_plot_localization.py → eg08_plot_localization.py

@@ -3,12 +3,13 @@
 """
 #%%
 from operator import itemgetter
-
+import numpy as np
 import numpy.linalg as la
 import matplotlib.pyplot as plt
 from microstructures import *
 from utilities import read_h5, construct_stress_localization
 
+
 np.random.seed(0)
 file_name, data_path, temp1, temp2, n_tests, sampling_alphas = itemgetter(
     "file_name", "data_path", "temp1", "temp2", "n_tests", "sampling_alphas"