Revert "add initial support for plastic models"

This reverts commit 2aa7c46.

Author: juliusgh

Diff for: eg2_compare_approximations.py

@@ -30,7 +30,7 @@
 
 _, samples = read_h5(file_name, data_path, [temp1, temp2], get_mesh=False)
 
-strains = np.random.normal(size=(n_loading_directions, strain_dof))
+strains = np.random.normal(size=(n_loading_directions, mesh['strain_dof']))
 strains /= la.norm(strains, axis=1)[:, None]
 
 n_approaches = 5
@@ -54,43 +54,38 @@
     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']
 
-    Sref = construct_stress_localization(Eref, ref_C, ref_eps, plastic_modes, mat_id, n_gauss, strain_dof)
+    Sref = construct_stress_localization(Eref, ref_C, ref_eps, 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, plastic_modes, mat_id, n_gauss, strain_dof)
+    Snaive = construct_stress_localization(Enaive, ref_C, ref_eps, 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, 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)
+    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)
     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, 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)
+    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)
     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, 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)
+    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)
     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, 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)
+    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)
     effSopt4 = volume_average(Sopt4)
 
     err = lambda x, y: np.mean(la.norm(x - y, axis=(-1, -2)) / la.norm(y, axis=(-1, -2))) * 100
@@ -102,7 +97,7 @@
                          err(effSopt2, effSref), err(effSopt4, effSref)]
 
     for strain_idx, strain in enumerate(strains):
-        zeta = np.hstack((strain, 1, np.ones(plastic_modes.shape[-1])))
+        zeta = np.hstack((strain, 1))
 
         eff_stress_ref = effSref @ zeta
         err_eff_stress[:, idx * n_loading_directions + strain_idx] = \
@@ -117,7 +112,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],
-                                               global_gradient)
+                                               mesh['global_gradient'])
 
         err_f[:, idx * n_loading_directions + strain_idx] = la.norm(residuals, np.inf, axis=0) / normalization_factor_mech * 100
 

Diff for: eg3_hierarchical_sampling.py

@@ -60,27 +60,22 @@
         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']
-        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
+        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, 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)
+        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)
         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, np.ones(plastic_modes.shape[-1])))
+            zeta = np.hstack((strain, 1))
             stress_opt4 = np.einsum('ijk,k', Sopt4, zeta, optimize='optimal')
             residual = compute_residual_efficient(stress_opt4, global_gradient)
 
@@ -94,11 +89,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[:, 7]), la.solve(Cref, effSref[:, 7]))
+        err_eff_eps[level, idx] = err(la.solve(Capprox, effSopt[:, -1]), la.solve(Cref, effSref[:, -1]))
 
     # 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.unique(np.sort(np.hstack((alphas, test_alphas[max_err_idx])))))
+    alpha_levels.append(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]) = }')
 

Diff for: eg4_hierarchical_sampling_efficient.py

@@ -6,8 +6,8 @@
 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, opt4
-from utilities import read_h5, construct_stress_localization, compute_err_indicator_efficient, volume_average
+from optimize_alpha import opt4_alphas
+from utilities import read_h5, construct_stress_localization, compute_err_indicator_efficient
 
 np.random.seed(0)
 # np.set_printoptions(precision=3)
@@ -40,20 +40,16 @@
 
     # 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))
     normalization_factor_mech = np.zeros((n_tests))
-    plastic_modes = refs[0]['plastic_modes']  # temperature independent
     for idx, alpha in enumerate(test_alphas):
-        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']
-        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)
+        effSref[idx] = np.hstack(
+            (refs[idx]['eff_stiffness'], -np.reshape(refs[idx]['eff_stiffness'] @ refs[idx]['eff_thermal_strain'], (-1, 1))))
 
     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)
@@ -109,26 +105,20 @@
                 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, plastic_modes, n_modes, n_phases,
-                    n_gp, strain_dof,
-                    mat_id, n_gauss)
+                    E01s[current_sampling_id], sampling_C, sampling_eps, 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, plastic_modes, mat_id, n_gauss,
-                                                           strain_dof,
+                _, effSopt = interpolate_fluctuation_modes(E01s[current_sampling_id], C, eps, mat_id, n_gauss, strain_dof,
                                                            n_modes, n_gp)
             elif speed == 1:
-                # TODO verify the result when plasticity is on
-                effSopt = effective_stress_localization(E01s[current_sampling_id], phi, ref_Cs[idx], ref_epss[idx], plastic_modes,
-                                                        mat_id,
+                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:
-                # 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
                 effSopt, phi = effective_S(phi, S001, S101, S103, S002, S102, S104, alpha_C, np.squeeze(alpha_eps, axis=-1),
@@ -138,9 +128,7 @@
 
             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], plastic_modes, mat_id, n_gauss,
-                                                      strain_dof)
-                # effSopt = volume_average(Sopt4)
+                Sopt4 = construct_stress_localization(Eopt4, ref_Cs[idx], ref_epss[idx], mat_id, n_gauss, strain_dof)
                 err_indicators[level,
                                idx] = np.mean(np.max(np.abs(compute_err_indicator_efficient(Sopt4, global_gradient)),
                                                      axis=0)) / normalization_factor_mech[idx] * 100
@@ -152,7 +140,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[:, 7]), la.solve(Cref, effSref[idx][:, 7]))
+            err_eff_eps[level, idx] = err(la.solve(Capprox, effSopt[:, -1]), la.solve(Cref, effSref[idx][:, -1]))
 
             # TODO remove dtype='f'
             group = file.require_group(f'{data_path}_level{level}')
@@ -161,7 +149,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[:, 7])
+            dset_thermal_strain[:] = la.solve(Capprox, effSopt[:, -1])
 
         if not given_alpha_levels:
             max_err_idx = np.argmax(err_indicators[level])

Diff for: interpolate_fluctuation_modes.py

@@ -2,8 +2,7 @@
 from numba import jit, prange
 
 @jit(nopython=True, cache=True, parallel=True, nogil=True)
-def interpolate_fluctuation_modes(E01, mat_stiffness, mat_thermal_strain, plastic_modes, mat_id, n_gauss, strain_dof, n_modes,
-                                  n_gp):
+def interpolate_fluctuation_modes(E01, mat_stiffness, mat_thermal_strain, mat_id, n_gauss, strain_dof, n_modes, n_gp):
     K = np.zeros((2 * n_modes, 2 * n_modes))
     F = np.zeros((2 * n_modes, n_modes))
     E_approx = np.zeros((n_gp, strain_dof, n_modes))
@@ -13,7 +12,7 @@ def interpolate_fluctuation_modes(E01, mat_stiffness, mat_thermal_strain, plasti
     for gp_id in prange(n_gp):
         phase_id = mat_id[gp_id // n_gauss]
 
-        P = np.hstack((-I, mat_thermal_strain[phase_id], plastic_modes[gp_id]))
+        P = np.hstack((-I, mat_thermal_strain[phase_id]))
 
         E01t_C = E01[gp_id].T @ mat_stiffness[phase_id]
 
@@ -30,8 +29,7 @@ def interpolate_fluctuation_modes(E01, mat_stiffness, mat_thermal_strain, plasti
 
     return E_approx, S_average[:6, :]
 
-def update_affine_decomposition(E01, sampling_C, sampling_eps, plastic_modes, n_modes, n_phases, n_gp, strain_dof, mat_id,
-                                n_gauss):
+def update_affine_decomposition(E01, sampling_C, sampling_eps, n_modes, n_phases, n_gp, strain_dof, mat_id, n_gauss):
     I = np.eye(strain_dof)
 
     K0 = np.zeros((2 * n_modes, 2 * n_modes))
@@ -61,7 +59,7 @@ def update_affine_decomposition(E01, sampling_C, sampling_eps, plastic_modes, n_
         eps0phase = eps0[phase_id]
         deltaCphase = deltaC[phase_id]
         deltaEPSphase = delta_eps[phase_id]
-        P0_phase = np.hstack((-I, eps0phase, plastic_modes[gp_id]))
+        P0_phase = np.hstack((-I, eps0phase))
         E01_transposed = E01[gp_id].T
 
         # essential terms that construct all other precomputed terms
@@ -118,12 +116,11 @@ def transform_strain_localization(E01, phi, n_gp, strain_dof, n_modes):
     return E_approx
 
 @jit(nopython=True, cache=True, parallel=True, nogil=True)
-def effective_stress_localization(E01, phi, mat_stiffness, mat_thermal_strain, plastic_modes, mat_id, n_gauss, n_gp, strain_dof,
-                                  n_modes):
+def effective_stress_localization(E01, phi, mat_stiffness, mat_thermal_strain, mat_id, n_gauss, n_gp, strain_dof, n_modes):
     S_average = np.zeros((strain_dof, n_modes))
     I = np.eye(strain_dof)
     for gp_id in prange(n_gp):
         phase_id = mat_id[gp_id // n_gauss]
-        P = np.hstack((-I, mat_thermal_strain[phase_id], plastic_modes[gp_id]))
+        P = np.hstack((-I, mat_thermal_strain[phase_id]))
         S_average += mat_stiffness[phase_id] @ (E01[gp_id] @ phi - P)
     return S_average / n_gp