#!python
import argparse
import pathlib
import sys
import h5py
import matplotlib.pyplot
import numpy
import pandas
from scipy.linalg import norm
from scipy.linalg import polar
import xdmf_reader_tools as XRT
[docs]
def str2bool(v):
'''Function for converting string to Boolean. Borrowed from: https://stackoverflow.com/questions/15008758/parsing-boolean-values-with-argparse
:param str/bool v: A string or boolean indicating a True or False value
:returns: True or False
'''
if isinstance(v, bool):
return v
if v.lower() in ('yes', 'true', 't', 'y', '1'):
return True
elif v.lower() in ('no', 'false', 'f', 'n', '0'):
return False
else:
raise argparse.ArgumentTypeError('Boolean value expected.')
[docs]
def plot_reference_stresses(E, stress, average, output_name, dim=3):
''' Plot stress vs strain in the reference configuration
:param dict E: The quantities dict storing Green-Lagrance strain
:param dict stress: The quantities dict storing either second Piola Kirchhoff or Symmetric micro stress
:param bool average: Average over quadrature points if True
:param str output_name: The output plot name
:param int dim: The dimension of the arrays, '2' for 2d or '3' for 3d
:returns: ``output_name``
'''
name = output_name.replace('.PNG','')
fig1 = matplotlib.pyplot.figure(name, figsize=(11,9))
axes1 = [[fig1.add_subplot(dim,dim,dim * i + j + 1) for j in range(dim)] for i in range(dim)]
ybounds = [-1, 1]
colors = matplotlib.pyplot.rcParams['axes.prop_cycle'].by_key()['color']
for i in range(dim):
for j in range(dim):
ax1 = axes1[i][j]
if 'PK2' in output_name:
plot_label = r"$PK2_{" + str(i+1) + str(j+1) + "}$ (MPa)"
if 'SIGMA' in output_name:
plot_label = r"$\Sigma_{" + str(i+1) + str(j+1) + "}$ (MPa)"
if average:
ax1.plot(E[0][:,0,i,j], stress[0][:,0,i,j], '-o')
else:
for qp in stress.keys():
ax1.plot(E[qp][:,0,i,j], stress[qp][:,0,i,j], '-o', color=colors[qp], label=f'qp #{qp+1}')
if (i == 2) and (j == 2):
matplotlib.pyplot.xticks(rotation=45)
ax1.set_xlabel(r"$E_{" + str(i+1) + str(j+1) + "}$", fontsize=14)
ax1.set_ylabel(plot_label, fontsize=14)
matplotlib.pyplot.ticklabel_format(style='sci', axis='x')
matplotlib.pyplot.ticklabel_format(style='sci', axis='y')
ax1.tick_params('x', labelrotation=45)
ax1 = axes1[dim-2][dim-1]
if average == False:
ax1.legend(bbox_to_anchor=(1.2, 0.9))
matplotlib.pyplot.figure(name)
matplotlib.pyplot.tight_layout()
matplotlib.pyplot.savefig(f'{output_name}')
return 0
[docs]
def plot_current_stresses(e, stress, average, output_name, dim=3):
''' Plot stress vs strain in the current configuration
:param dict e: The quantities dict storing Euler-Almansi strain
:param dict stress: The quantities dict storing either Cauchy or symmetric micro stress
:param bool average: Average over quadrature points if True
:param str output_name: The output plot name
:param int dim: The dimension of the arrays, '2' for 2d or '3' for 3d
:returns: ``output_name``
'''
name = output_name.replace('.PNG','')
fig1 = matplotlib.pyplot.figure(name, figsize=(11,9))
axes1 = [[fig1.add_subplot(dim,dim,dim * i + j + 1) for j in range(dim)] for i in range(dim)]
ybounds = [-1, 1]
colors = matplotlib.pyplot.rcParams['axes.prop_cycle'].by_key()['color']
for i in range(dim):
for j in range(dim):
ax1 = axes1[i][j]
if 'cauchy' in output_name:
plot_label = r"$\sigma_{" + str(i+1) + str(j+1) + "}$ (MPa)"
if 'symm' in output_name:
plot_label = r"$s_{" + str(i+1) + str(j+1) + "}$ (MPa)"
if average:
ax1.plot(e[0][:,0,i,j], stress[0][:,0,i,j], '-o')
else:
for qp in stress.keys():
ax1.plot(e[qp][:,0,i,j], stress[qp][:,0,i,j], '-o', color=colors[qp], label=f'qp #{qp+1}')
if (i == 2) and (j ==2):
matplotlib.pyplot.xticks(rotation=45)
ax1.set_xlabel(r"$e_{" + str(i+1) + str(j+1) + "}$", fontsize=14)
ax1.set_ylabel(plot_label, fontsize=14)
matplotlib.pyplot.ticklabel_format(style='sci', axis='x')
matplotlib.pyplot.ticklabel_format(style='sci', axis='y')
ax1.tick_params('x', labelrotation=45)
ax1 = axes1[dim-2][dim-1]
if average == False:
ax1.legend(bbox_to_anchor=(1.2, 0.9))
matplotlib.pyplot.figure(name)
matplotlib.pyplot.tight_layout()
matplotlib.pyplot.savefig(f'{output_name}')
return 0
[docs]
def plot_axial_vector_mag(cauchy, times, average, output_name, dim=3):
''' Plot the magnitude of the Cauchy couple
:param dict cauchy: The quantities dict storing Cauchy stress
:param array-like times: The time increments
:param bool average: Average over quadrature points if True
:param str output_name: The output plot name
:param int dim: The dimension of the arrays, '2' for 2d or '3' for 3d
:returns: ``output_name``
'''
colors = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd', '#8c564b', '#e377c2', '#7f7f7f', '#bcbd22', '#17becf']
fig1 = matplotlib.pyplot.figure('axial', figsize=(5,5),dpi=300)
nel = numpy.shape(cauchy[0])[1]
for el in range(nel):
for qp in cauchy.keys():
axial = []
for i,t in enumerate(times):
# skew symmetric
sig = cauchy[qp][i,el,:,:]
sig_t = numpy.transpose(sig)
skew = 0.5*(sig - sig_t)
# off diagonals
if dim == 3:
ax = numpy.sqrt((skew[0][1]**2)+(skew[0][2]**2)+(skew[1][2]**2))
elif dim == 2:
ax = skew[0][1]
# Append results
axial.append(ax)
# Plot magnitude of axial vector
if average:
matplotlib.pyplot.plot(times, axial)
else:
matplotlib.pyplot.plot(times, axial, color=colors[qp])
matplotlib.pyplot.figure('axial')
matplotlib.pyplot.xlabel('Simulation time (s)', fontsize=14)
matplotlib.pyplot.ylabel('|$\omega^{\sigma}$| (MPa)', fontsize=14)
matplotlib.pyplot.tight_layout()
matplotlib.pyplot.savefig(f'{output_name}')
return 0
[docs]
def plot_stress_diffs(cauchy, symm, times, average, output_name, dim=3):
''' Plot differences between Cauchy and micro symmetric stresses
:param dict cauchy: The quantities dict storing Cauchy stress
:param dict symm: The quantities dict storing symmetric micro stress
:param array-like times: The time increments
:param bool average: Average over quadrature points if True
:param str output_name: The output plot name
:param int dim: The dimension of the arrays, '2' for 2d or '3' for 3d
:returns: ``output_name``
'''
fig1 = matplotlib.pyplot.figure('stress_diff', figsize=(11,9),dpi=300)
axes1 = [[fig1.add_subplot(dim,dim,dim * i + j + 1) for j in range(dim)] for i in range(dim)]
colors = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd', '#8c564b', '#e377c2', '#7f7f7f', '#bcbd22', '#17becf']
nel = numpy.shape(cauchy[0])[1]
for el in range(nel):
for i in range(dim):
for j in range(dim):
ax1 = axes1[i][j]
for qp in cauchy.keys():
sig = cauchy[qp][:,el,i,j]
sym = symm[qp][:,el,i,j]
diff = sig-sym
if average:
ax1.plot(times, diff)
else:
ax1.plot(times, diff, color=colors[qp])
ax1.set_xlabel(r'Simulation time (s)', fontsize=14)
ax1.set_ylabel(r"($\sigma_{" + str(i+1) + str(j+1) + "}$ - $s_{" + str(i+1) + str(j+1) + "}$)", fontsize=14)
ax1 = axes1[dim-2][dim-1]
matplotlib.pyplot.figure('stress_diff')
matplotlib.pyplot.tight_layout()
matplotlib.pyplot.savefig(f'{output_name}')
return 0
[docs]
def get_R_and_U(PK2, F, chi, times, dim=3):
''' Calculate macro and micro stretches and deformations from deformation gradients
:param dict PK2: The quantities dict storing second Piola Kirchhoff stress
:param dict F: The quantities dict storing the macro deformation gradient
:param dict chi: The quantities dict storing the micro deformation
:param array-like times: The time increments
:param int dim: The dimension of the arrays, '2' for 2d or '3' for 3d
:returns: macro rotation (R), macro stretch (U), micro rotation (Rchi), and micro stretch (Uchi)
'''
R, U = {}, {}
Rchi, Uchi = {}, {}
nel = numpy.shape(PK2[0])[1]
ntimes = len(times)
for qp in PK2.keys():
Rs, Us = numpy.zeros([ntimes,nel,dim,dim]), numpy.zeros([ntimes,nel,dim,dim])
Rcs, Ucs = numpy.zeros([ntimes,nel,dim,dim]), numpy.zeros([ntimes,nel,dim,dim])
for el in range(nel):
for i,t in enumerate(times):
#if i > 0:
Rr, Uu = polar(F[qp][i,el,:,:],side='left')
Rcr, Ucu = polar(chi[qp][i,el,:,:], side='left')
#check orthogonality
tol = 1.e-9
if (norm((numpy.dot(Rr,numpy.transpose(Rr)) - numpy.eye(dim)),ord=2)) > tol:
print('Error!!! R is not orthogonal!')
if (norm((numpy.dot(Rcr,numpy.transpose(Rcr)) - numpy.eye(dim)),ord=2)) > tol:
print('Error!!! Rc is not orthogonal!')
Rs[i,el,:,:] = Rr
Us[i,el,:,:] = Uu
Rcs[i,el,:,:] = Rcr
Ucs[i,el,:,:] = Ucu
R[qp] = Rs
U[qp] = Us
Rchi[qp] = Rcs
Uchi[qp] = Ucs
return(R,U,Rchi,Uchi)
[docs]
def plot_rot_diffs(PK2, times, R, Rchi, average, output_name, dim=3):
''' Plot differences between macro and micro rotations
:param dict PK2: The quantities dict storing second Piola Kirchhoff stress
:param array-like times: The time increments
:param dict R: The quantities dict storing macro rotations
:param dict Rchi: The quantities dict storing micro rotations
:param bool average: Average over quadrature points if True
:param str output_name: The output plot name
:param int dim: The dimension of the arrays, '2' for 2d or '3' for 3d
:returns: ``output_name``
'''
fig1 = matplotlib.pyplot.figure('Rot_diff', figsize=(12,9))
axes1 = [[fig1.add_subplot(dim,dim,dim * i + j + 1) for j in range(dim)] for i in range(dim)]
colors = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd', '#8c564b', '#e377c2', '#7f7f7f', '#bcbd22', '#17becf']
nel = numpy.shape(PK2[0])[1]
for el in range(nel):
for i in range(dim):
for j in range(dim):
ax1 = axes1[i][j]
for qp in PK2.keys():
diffs = []
for k,t in enumerate(times):
diff = R[qp][k,el,i,j] - Rchi[qp][k,el,i,j]
diffs.append(diff)
if average:
ax1.plot(times[:], diffs)
else:
ax1.plot(times[:], diffs, color=colors[qp])
ax1.set_xlabel(r'Simulation time (s)', fontsize=14)
ax1.set_ylabel(r"$R_{" + str(i+1) + str(j+1) + "}$ - $R^{\chi}_{" + str(i+1) + str(j+1) + "}$", fontsize=14)
ax1.ticklabel_format(axis='y',style='sci')
matplotlib.pyplot.ticklabel_format(axis='y',style='sci')
matplotlib.pyplot.figure('Rot_diff')
matplotlib.pyplot.tight_layout()
matplotlib.pyplot.savefig(f'{output_name}')
return 0
[docs]
def plot_stretch_diffs(PK2, times, U, Uchi, average, output_name, dim=3):
''' Plot differences between macro and micro stretches
:param dict PK2: The quantities dict storing second Piola Kirchhoff stress
:param array-like times: The time increments
:param dict U: The quantities dict storing macro stretches
:param dict Uchi: The quantities dict storing micro stretches
:param bool average: Average over quadrature points if True
:param str output_name: The output plot name
:param int dim: The dimension of the arrays, '2' for 2d or '3' for 3d
:returns: ``output_name``
'''
fig1 = matplotlib.pyplot.figure('U', figsize=(11,9))
axes1 = [[fig1.add_subplot(dim,dim,dim * i + j + 1) for j in range(dim)] for i in range(dim)]
colors = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd', '#8c564b', '#e377c2', '#7f7f7f', '#bcbd22', '#17becf']
nel = numpy.shape(PK2[0])[1]
for el in range(nel):
for i in range(dim):
for j in range(dim):
ax1 = axes1[i][j]
for qp in PK2.keys():
diffs = []
for k,t in enumerate(times):
diff = U[qp][k,el,i,j] - Uchi[qp][k,el,i,j]
diffs.append(diff)
if average:
ax1.plot(times[:], diffs)
else:
ax1.plot(times[:],diffs, color=colors[qp])
ax1.set_xlabel(r'Simulation time (s)', fontsize=14)
ax1.set_ylabel(r"$U_{" + str(i+1) + str(j+1) + "}$ - $U^{\chi}_{" + str(i+1) + str(j+1) + "}$", fontsize=14)
ax1.ticklabel_format(axis='y',style='sci')
matplotlib.pyplot.figure('U')
matplotlib.pyplot.tight_layout()
matplotlib.pyplot.savefig(f'{output_name}')
return 0
[docs]
def plot_first_moment_of_momentum_measures(coup, spins, times, average,
plot_body_couples=None, plot_spin_inertias=None, plot_spin_diff=None,
dim=3):
'''Plot the first moment of momentum measures
:param dict coup: The quantities dict storing body couples
:param dict spins: THe quantities dict stroring micro spin inertias
:param array-like times: The time increments
:param bool average: Average over quadrature points if True
:param str plot_body_couples: Optional filename to plot body couples vs. simulation time
:param str plot_spin_inertias: Optional filename to plot micro spin inertias vs. simulation time
:param str plot_spin_diff: Optional filename to plot difference between body couples and micro spin inertias
:param int dim: The dimension of the arrays, '2' for 2d or '3' for 3d
:returns: ``plot_body_couples``, ``plot_spin_inertias``, and ``plot_spin_diff``
'''
name1 = 'body_couples'
name2 = 'micro_spin_inertias'
name3 = 'couple_spin_diff'
fig1 = matplotlib.pyplot.figure(name1, figsize=(11,9))
axes1 = [[fig1.add_subplot(dim,dim,dim * i + j + 1) for j in range(dim)] for i in range(dim)]
fig2 = matplotlib.pyplot.figure(name2, figsize=(11,9))
axes2 = [[fig2.add_subplot(dim,dim,dim * i + j + 1) for j in range(dim)] for i in range(dim)]
fig3 = matplotlib.pyplot.figure(name3, figsize=(11,9))
axes3 = [[fig3.add_subplot(dim,dim,dim * i + j + 1) for j in range(dim)] for i in range(dim)]
colors = matplotlib.pyplot.rcParams['axes.prop_cycle'].by_key()['color']
for i in range(dim):
for j in range(dim):
ax1 = axes1[i][j]
ax2 = axes2[i][j]
ax3 = axes3[i][j]
if average:
ax1.plot(times[:], coup[0][:,0,i,j], '-o')
ax2.plot(times[:], spins[0][:,0,i,j], '-o')
ax3.plot(times[:], coup[0][:,0,i,j] - spins[0][:,0,i,j], '-o')
else:
for qp in coup.keys():
ax1.plot(times[:], coup[qp][:,0,i,j], '-o', color=colors[qp], label=f'qp #{qp+1}')
ax2.plot(times[:], spins[qp][:,0,i,j], '-o', color=colors[qp], label=f'qp #{qp+1}')
ax3.plot(times[:], coup[qp][:,0,i,j] - spins[qp][:,0,i,j], '-o', color=colors[qp], label=f'qp #{qp+1}')
ax1.set_xlabel(r'Simulation time (s)', fontsize=14)
ax1.set_ylabel(r'$l_{' + str(i+1) + str(j+1) + '}$', fontsize=14)
ax2.set_xlabel(r'Simulation time (s)', fontsize=14)
ax2.set_ylabel(r'$\omega_{' + str(i+1) + str(j+1) + '}$', fontsize=14)
ax3.set_xlabel(r'Simulation time (s)', fontsize=14)
ax3.set_ylabel(r'$l_{' + str(i+1) + str(j+1) + '}$ - $\omega_{' + str(i+1) + str(j+1) + '}$', fontsize=14)
ax1, ax2, ax3 = axes1[dim-2][dim-1], axes2[dim-2][dim-1], axes3[dim-2][dim-1]
if average == False:
ax1.legend(bbox_to_anchor=(1.2,0.9))
ax2.legend(bbox_to_anchor=(1.2,0.9))
ax3.legend(bbox_to_anchor=(1.2,0.9))
matplotlib.pyplot.figure(name1)
matplotlib.pyplot.tight_layout()
matplotlib.pyplot.figure(name2)
matplotlib.pyplot.tight_layout()
matplotlib.pyplot.figure(name3)
matplotlib.pyplot.tight_layout()
if plot_body_couples:
fig1.savefig(plot_body_couples)
if plot_spin_inertias:
fig2.savefig(plot_spin_inertias)
if plot_spin_diff:
fig3.savefig(plot_spin_diff)
return 0
[docs]
def deviatoric(stress, dim=3):
'''Calculate the deviatoric component of a stress quantity
:param array-like stress: A second order tensor or 9-component slice of a third order tensor
:param int dim: The dimension of the arrays, '2' for 2d or '3' for 3d
:returns: deviatoric of ``stress``
'''
return(stress - (1/dim)*numpy.eye(dim)*numpy.trace(stress))
[docs]
def plot_stress_norm(cauchy, symm, m, nqp, nel, ninc, times, output_name, dim=3):
'''Plot the infinity norms of deviatoric Cauchy, symmetric micro, and higher stresses
:param dict cauchy: The quantities dict storing Cauchy stress
:param dict symm: The quantities dict storing symmetric micro stress
:param dict m: THe quantities dict storing higher order stress
:param int nqp: The number of quadrature points
:param int nel: The number of elements
:param int ninc: The number of time increments
:param array-like times: The time increments
:param str output_name: Output filename
:param int dim: The dimension of the arrays, '2' for 2d or '3' for 3d
:returns: ``output_name`` plot
'''
fig = matplotlib.pyplot.figure('stress_norms', figsize=(9,9), dpi=300)
ax1 = fig.add_subplot(2,2,1)
ax2 = fig.add_subplot(2,2,2)
ax3 = fig.add_subplot(2,2,3)
ax4 = fig.add_subplot(2,2,4)
for el in range(nel):
for qp in range(nqp):
# take norms
cauchy_norm, symm_norm, diff_norm = [], [], []
m1_norm, m2_norm, m3_norm = [], [], []
for t in range(ninc):
cauchy_norm.append(numpy.linalg.norm(deviatoric(cauchy[qp][t,el,:,:], dim=dim), ord='fro'))
symm_norm.append(numpy.linalg.norm(deviatoric(symm[qp][t,el,:,:], dim=dim), ord='fro'))
diff_norm.append(numpy.linalg.norm(deviatoric(cauchy[qp][t,el,:,:]-symm[qp][t,el,:,:], dim=dim), ord='fro'))
m1_norm.append(numpy.linalg.norm(deviatoric(m[qp][t,el,:,:,0], dim=dim), ord='fro'))
m2_norm.append(numpy.linalg.norm(deviatoric(m[qp][t,el,:,:,1], dim=dim), ord='fro'))
if dim == 3:
m3_norm.append(numpy.linalg.norm(deviatoric(m[qp][t,el,:,:,2], dim=dim), ord='fro'))
label = f"qp #{(qp+1)+(nel*nqp)}"
ax1.plot(times, cauchy_norm, label=label)
ax1.set_xlabel(r'Simulation time (s)', fontsize=14)
ax1.set_ylabel(r'$||dev\left(\sigma_{ij}\right)||$', fontsize=14)
ax2.plot(times, symm_norm, label=label)
ax2.set_ylabel(r'$||dev\left(s_{ij}\right)||$', fontsize=14)
ax2.set_xlabel(r'Simulation time (s)', fontsize=14)
ax3.plot(times, diff_norm, label=label)
ax3.set_xlabel(r'Simulation time (s)', fontsize=14)
ax3.set_ylabel(r'$||dev\left(\sigma_{ij} - s_{ij}\right)||$', fontsize=14)
ax4.plot(times, m1_norm, '-o', label=f"index k = 1, {label}")
ax4.plot(times, m2_norm, '-^', label=f"index k = 2, {label}")
if dim == 3:
ax4.plot(times, m3_norm, '-v', label=f"index k = 3, {label}")
ax4.set_xlabel(r'Simulation time (s)', fontsize=14)
ax4.set_ylabel(r'$||dev\left(m_{ijk}\right)||$', fontsize=14)
# ax4.legend([(label='index 1', color='b'),
# (label='index 2', color='b'),
# (label='index 3', color='b')])
if dim == 3:
ax4.legend(['index k = 1', 'index k = 2', 'index k = 3'])
elif dim == 2:
ax4.legend(['index k = 1', 'index k = 2', 'index k = 3'])
ax1.set_xlabel(r'Simulation time (s)', fontsize=14)
fig.savefig(output_name)
return 0
[docs]
def triple_deviatoric(stress, dim=3):
'''Calculate the deviatoric component of a stress quantity
:param array-like stress: A second order tensor or (dim*dim)-component slice of a third order tensor
:param int dim: The dimension of the arrays, '2' for 2d or '3' for 3d
:returns: deviatoric of ``stress``
'''
stress_1 = stress[:,:,0] - (1/dim)*numpy.eye(dim)*numpy.trace(stress[:,:,0])
stress_2 = stress[:,:,1] - (1/dim)*numpy.eye(dim)*numpy.trace(stress[:,:,1])
if dim == 3:
stress_3 = stress[:,:,2] - (1/dim)*numpy.eye(dim)*numpy.trace(stress[:,:,2])
return(stress_1 + stress_2 + stress_3)
elif dim == 2:
return(stress_1 + stress_2)
[docs]
def plot_better_stress_norm(PK2, SIGMA, M, E, Ecal, Gamma, nqp, nel, ninc, output_name, dim=3):
'''Plot the infinity norms of deviatoric Cauchy, symmetric micro, and higher stresses
:param dict PK2: The quantities dict storing second Piola-Kirchhoff stress
:param dict SIGMA: The quantities dict storing symmetric micro stress
:param dict M: THe quantities dict storing higher order stress
:param dict E: The quantities dict storing Green-Lagrange strain
:param dict Ecal: The quantities dict storing micro-strain
:param dict Gamma: The quantities dict storing micro deformation gradient
:param int nqp: The number of quadrature points
:param int nel: The number of elements
:param int ninc: The number of time increments
:param str output_name: Output filename
:param int dim: The dimension of the arrays, '2' for 2d or '3' for 3d
:returns: ``output_name`` plot
'''
fig, axes = matplotlib.pyplot.subplots(4, 3)
for el in range(nel):
for qp in range(nqp):
# take norms
PK2_norm, SIGMA_norm, diff_norm = [], [], []
M1_norm, M2_norm, M3_norm = [], [], []
E_norm, Ecal_norm, Gamma_norm = [], [], []
for t in range(ninc):
PK2_norm.append(numpy.linalg.norm(deviatoric(PK2[qp][t,el,:,:], dim=dim), ord='fro'))
SIGMA_norm.append(numpy.linalg.norm(deviatoric(SIGMA[qp][t,el,:,:], dim=dim), ord='fro'))
diff_norm.append(numpy.linalg.norm(deviatoric(PK2[qp][t,el,:,:]-SIGMA[qp][t,el,:,:], dim=dim), ord='fro'))
M1_norm.append(numpy.linalg.norm(deviatoric(M[qp][t,el,:,:,0], dim=dim), ord='fro'))
M2_norm.append(numpy.linalg.norm(deviatoric(M[qp][t,el,:,:,1], dim=dim), ord='fro'))
if dim == 3:
M3_norm.append(numpy.linalg.norm(deviatoric(M[qp][t,el,:,:,2], dim=dim), ord='fro'))
E_norm.append(numpy.linalg.norm(deviatoric(E[qp][t,el,:,:], dim=dim), ord='fro'))
Ecal_norm.append(numpy.linalg.norm(deviatoric(Ecal[qp][t,el,:,:], dim=dim), ord='fro'))
Gamma_norm.append(numpy.linalg.norm(triple_deviatoric(Gamma[qp][t,el,:,:,:], dim=dim), ord='fro'))
label = f"qp #{(qp+1)+(nel*nqp)}"
axes[0][0].plot(E_norm, PK2_norm)
axes[1][0].plot(E_norm, SIGMA_norm)
axes[2][0].plot(E_norm, diff_norm)
axes[3][0].plot(E_norm, M1_norm, '-o', label=f"index K = 1, {label}")
axes[3][0].plot(E_norm, M2_norm, '-^', label=f"index K = 2, {label}")
if dim == 3:
axes[3][0].plot(E_norm, M3_norm, '-v', label=f"index K = 3, {label}")
axes[0][1].plot(Ecal_norm, PK2_norm)
axes[1][1].plot(Ecal_norm, SIGMA_norm)
axes[2][1].plot(Ecal_norm, diff_norm)
axes[3][1].plot(Ecal_norm, M1_norm, '-o', label=f"index K = 1, {label}")
axes[3][1].plot(Ecal_norm, M2_norm, '-^', label=f"index K = 2, {label}")
if dim == 3:
axes[3][1].plot(Ecal_norm, M3_norm, '-v', label=f"index K = 3, {label}")
axes[0][2].plot(Gamma_norm, PK2_norm)
axes[1][2].plot(Gamma_norm, SIGMA_norm)
axes[2][2].plot(Gamma_norm, diff_norm)
axes[3][2].plot(Gamma_norm, M1_norm, '-o', label=f"index K = 1, {label}")
axes[3][2].plot(Gamma_norm, M2_norm, '-^', label=f"index K = 2, {label}")
if dim == 3:
axes[3][2].plot(Gamma_norm, M3_norm, '-v', label=f"index K = 3, {label}")
axes[0][0].set_ylabel(r'$||dev\left(S_{IJ}\right)|| \left( MPa \right)$', fontsize=14)
axes[1][0].set_ylabel(r'$||dev\left(\Sigma_{IJ}\right)|| \left( MPa \right)$', fontsize=14)
axes[2][0].set_ylabel(r'$||dev\left(S_{IJ} - \Sigma_{IJ}\right)|| \left( MPa \right)$', fontsize=14)
axes[3][0].set_ylabel(r'$||dev\left(M_{IJK}\right)|| \left( MPa \cdot mm^2 \right)$', fontsize=14)
axes[3][0].set_xlabel(r'$||dev\left(E_{IJ}\right)||$', fontsize=14)
axes[3][1].set_xlabel(r'$||dev\left(\mathcal{E}_{IJ}\right)||$', fontsize=14)
axes[3][2].set_xlabel(r'$||dev\left(\Gamma_{IJK}\right)||$', fontsize=14)
fig.set_figheight(16)
fig.set_figwidth(12)
fig.tight_layout()
fig.savefig(output_name)
return 0
[docs]
def filter_stress_above_threshold(stress_lists, strain_lists, threshold, elem_qp_string):
'''Filter out stresses above some threshold for a given element and quadrature point
:param list stress_lists: List of arrays/lists containing stress histories
:param list strain_lists: List of array/lists containing strain histories
:param float threshold: Upper limit threshold for filtering stresses
:param str elem_qp_string: The element id and quadrature point id for the provided stress and strain lists
:returns: Stress lists filtered below threshold, strain lists masked by filtered stress lists
'''
stress_norms, strain_norms = [], []
stress_names = ['PK2', 'SIGMA', 'M1', 'M2', 'M3']
for stress, strain, stress_name in zip(stress_lists, strain_lists, stress_names):
mask = numpy.array(stress) < threshold
if False in mask:
print(f'Oh no!!! {stress_name} stress above threshold for element-qp pair {elem_qp_string}')
print(f'\t {stress_name} = {stress}')
stress_out = numpy.array(stress)[mask]
strain_out = numpy.array(strain)[mask]
else:
stress_out = stress
strain_out = strain
stress_norms.append(stress_out)
strain_norms.append(strain_out)
return stress_norms, strain_norms
[docs]
def plot_best_stress_norm(PK2, SIGMA, M, E, Ecal, Gamma, nqp, nel, ninc, output_name,
transparency=None, dim=3, upper_limit=None):
'''Plot the infinity norms of deviatoric second Piola-Kirchhoff, symmetric micro, and higher order stresses versus relevant deformation measures
:param dict cauchy: The quantities dict storing Cauchy stress
:param dict symm: The quantities dict storing symmetric micro stress
:param dict m: THe quantities dict storing higher order stress
:param int nqp: The number of quadrature points
:param int nel: The number of elements
:param array-like times: The time increments
:param str output_name: Output filename
:param int dim: The dimension of the arrays, '2' for 2d or '3' for 3d
:param float upper_limit: Optional value above which stresses will be ignored
:returns: ``output_name`` plot
'''
fig, axes = matplotlib.pyplot.subplots(1, 3)
for el in range(nel):
for qp in range(nqp):
# take norms
PK2_norm, SIGMA_norm, diff_norm = [], [], []
M1_norm, M2_norm, M3_norm = [], [], []
E_norm, Ecal_norm, Gamma_norm = [], [], []
for t in range(ninc):
PK2_norm.append(numpy.linalg.norm(deviatoric(PK2[qp][t,el,:,:], dim=dim), ord='fro'))
SIGMA_norm.append(numpy.linalg.norm(deviatoric(SIGMA[qp][t,el,:,:], dim=dim), ord='fro'))
diff_norm.append(numpy.linalg.norm(deviatoric(PK2[qp][t,el,:,:]-SIGMA[qp][t,el,:,:], dim=dim), ord='fro'))
M1_norm.append(numpy.linalg.norm(deviatoric(M[qp][t,el,:,:,0], dim=dim), ord='fro'))
M2_norm.append(numpy.linalg.norm(deviatoric(M[qp][t,el,:,:,1], dim=dim), ord='fro'))
if dim == 3:
M3_norm.append(numpy.linalg.norm(deviatoric(M[qp][t,el,:,:,2], dim=dim), ord='fro'))
E_norm.append(numpy.linalg.norm(deviatoric(E[qp][t,el,:,:], dim=dim), ord='fro'))
Ecal_norm.append(numpy.linalg.norm(deviatoric(Ecal[qp][t,el,:,:], dim=dim), ord='fro'))
Gamma_norm.append(numpy.linalg.norm(triple_deviatoric(Gamma[qp][t,el,:,:,:], dim=dim), ord='fro'))
label = f"qp #{(qp+1)+(nel*nqp)}"
# Filter out bad points if requested
if upper_limit is not None:
PK2_mask = numpy.array(PK2_norm) < upper_limit
SIGMA_mask = numpy.array(SIGMA_norm) < upper_limit
M1_mask = numpy.array(M1_norm) < upper_limit
M2_mask = numpy.array(M2_norm) < upper_limit
M3_mask = numpy.array(M3_norm) < upper_limit
stress_norms, strain_norms = filter_stress_above_threshold(
[PK2_norm, SIGMA_norm, M1_norm, M2_norm, M3_norm],
[E_norm, Ecal_norm, Gamma_norm, Gamma_norm, Gamma_norm],
upper_limit, f'el {el} qp {qp}')
PK2_norm, SIGMA_norm, M1_norm, M2_norm, M3_norm = stress_norms
E_norm, Ecal_norm, Gamma_norm_1, Gamma_norm_2, Gamma_norm_3 = strain_norms
else:
Gamma_norm_1, Gamma_norm_2, Gamma_norm_3 = Gamma_norm, Gamma_norm, Gamma_norm
if transparency:
axes[0].plot(E_norm, PK2_norm, color='r', alpha=transparency)
axes[1].plot(Ecal_norm, SIGMA_norm, color='r', alpha=transparency)
axes[2].plot(Gamma_norm, M1_norm, '-o', label=f"index K = 1, {label}", color='r', alpha=transparency)
axes[2].plot(Gamma_norm, M2_norm, '-^', label=f"index K = 2, {label}", color='r', alpha=transparency)
if dim == 3:
axes[2].plot(Gamma_norm, M3_norm, '-v', label=f"index K = 3, {label}", color='r', alpha=transparency)
else:
axes[0].plot(E_norm, PK2_norm)
axes[1].plot(Ecal_norm, SIGMA_norm)
axes[2].plot(Gamma_norm, M1_norm, '-o', label=f"index K = 1, {label}")
axes[2].plot(Gamma_norm, M2_norm, '-^', label=f"index K = 2, {label}")
if dim == 3:
axes[2].plot(Gamma_norm, M3_norm, '-v', label=f"index K = 3, {label}")
axes[0].set_ylabel(r'$||dev\left(S_{IJ}\right)|| \left( MPa \right)$', fontsize=14)
axes[0].set_xlabel(r'$||dev\left(E_{IJ}\right)||$', fontsize=14)
axes[1].set_ylabel(r'$||dev\left(\Sigma_{IJ}\right)|| \left( MPa \right)$', fontsize=14)
axes[1].set_xlabel(r'$||dev\left(\mathcal{E}_{IJ}\right)|| $', fontsize=14)
axes[2].set_ylabel(r'$||dev\left(M_{IJK}\right)|| \left( MPa \cdot mm^2 \right)$', fontsize=14)
axes[2].set_xlabel(r'$||dev\left(\Gamma_{IJK}\right)||$', fontsize=14)
fig.set_figheight(5)
fig.set_figwidth(12)
fig.tight_layout()
fig.savefig(output_name)
return 0
[docs]
def plot_norm_history(PK2, SIGMA, M, E, Ecal, Gamma, nqp, nel, ninc, times, output_name, dim=3):
'''Plot the infinity norms of deviatoric second Piola-Kirchhoff, symmetric micro, and higher order stresses versus time
:param dict PK2: The quantities dict storing second Piola-Kirchhoff stress
:param dict SIGMA: The quantities dict storing symmetric micro stress
:param dict M: THe quantities dict storing higher order stress
:param dict E: The quantities dict storing Green-Lagrange strain
:param dict Ecal: The quantities dict storing micro-strain
:param dict Gamma: The quantities dict storing micro deformation gradient
:param int nqp: The number of quadrature points
:param int nel: The number of elements
:param int ninc: The number of time increments
:param array-like times: The time increments
:param str output_name: Output filename
:param int dim: The dimension of the arrays, '2' for 2d or '3' for 3d
:returns: ``output_name`` plot
'''
fig, axes = matplotlib.pyplot.subplots(3, 3)
for el in range(nel):
for qp in range(nqp):
# take norms
E_norm, Ecal_norm, Gamma_norm = [], [], []
PK2_p, PK2_q, SIGMA_p, SIGMA_q =[], [], [], []
diff_p, diff_q = [], []
M1_p, M1_q, M2_p, M2_q, M3_p, M3_q = [], [], [], [], [], []
for t in range(ninc):
E_norm.append(numpy.linalg.norm(deviatoric(E[qp][t,el,:,:], dim=dim), ord='fro'))
Ecal_norm.append(numpy.linalg.norm(deviatoric(Ecal[qp][t,el,:,:], dim=dim), ord='fro'))
Gamma_norm.append(numpy.linalg.norm(triple_deviatoric(Gamma[qp][t,el,:,:,:], dim=dim), ord='fro'))
PK2_p.append((-1/dim)*numpy.trace(PK2[qp][t,el,:,:]))
PK2_q.append(numpy.sqrt(numpy.tensordot(deviatoric(PK2[qp][t,el,:,:], dim=dim), deviatoric(PK2[qp][t,el,:,:], dim=dim))))
SIGMA_p.append((-1/dim)*numpy.trace(SIGMA[qp][t,el,:,:]))
SIGMA_q.append(numpy.sqrt(numpy.tensordot(deviatoric(SIGMA[qp][t,el,:,:], dim=dim), deviatoric(SIGMA[qp][t,el,:,:], dim=dim))))
diff_p.append((-1/dim)*numpy.trace(PK2[qp][t,el,:,:]-SIGMA[qp][t,el,:,:]))
diff_q.append(numpy.sqrt(numpy.tensordot(deviatoric(PK2[qp][t,el,:,:]-SIGMA[qp][t,el,:,:], dim=dim),
deviatoric(PK2[qp][t,el,:,:]-SIGMA[qp][t,el,:,:], dim=dim))))
M1_p.append((-1/dim)*numpy.trace(M[qp][t,el,:,:,0]))
M1_q.append(numpy.sqrt(numpy.tensordot(deviatoric(M[qp][t,el,:,:,0], dim=dim), deviatoric(M[qp][t,el,:,:,0], dim=dim))))
M2_p.append((-1/dim)*numpy.trace(M[qp][t,el,:,:,1]))
M2_q.append(numpy.sqrt(numpy.tensordot(deviatoric(M[qp][t,el,:,:,1], dim=dim), deviatoric(M[qp][t,el,:,:,1], dim=dim))))
if dim == 3:
M3_p.append((-1/dim)*numpy.trace(M[qp][t,el,:,:,2]))
M3_q.append(numpy.sqrt(numpy.tensordot(deviatoric(M[qp][t,el,:,:,2], dim=dim), deviatoric(M[qp][t,el,:,:,1], dim=dim))))
label = f"qp #{(qp+1)+(nel*nqp)}"
axes[0][0].plot(times, PK2_p)
axes[0][1].plot(times, PK2_q)
axes[0][2].plot(times, E_norm)
axes[1][0].plot(times, SIGMA_p)
axes[1][1].plot(times, SIGMA_q)
axes[1][2].plot(times, Ecal_norm)
axes[2][0].plot(times, M1_p, '-o', label=f"index K = 1, {label}")
axes[2][0].plot(times, M2_p, '-^', label=f"index K = 2, {label}")
axes[2][1].plot(times, M1_q, '-o', label=f"index K = 1, {label}")
axes[2][1].plot(times, M2_q, '-^', label=f"index K = 2, {label}")
if dim == 3:
axes[2][0].plot(times, M3_p, '-v', label=f"index K = 3, {label}")
axes[2][1].plot(times, M3_q, '-v', label=f"index K = 3, {label}")
axes[2][2].plot(times, Gamma_norm)
axes[0][0].set_ylabel(r'$p\left(\mathbf{S}\right) (MPa)$', fontsize=14)
axes[0][1].set_ylabel(r'$q\left(\mathbf{S}\right) (MPa)$', fontsize=14)
axes[0][2].set_ylabel(r'$||dev\left(E_{IJ}\right)||$', fontsize=14)
axes[1][0].set_ylabel(r'$p\left(\mathbf{\Sigma}\right) (MPa)$', fontsize=14)
axes[1][1].set_ylabel(r'$q\left(\mathbf{\Sigma}\right) (MPa)$', fontsize=14)
axes[1][2].set_ylabel(r'$||dev\left(\mathcal{E}_{IJ}\right)||$', fontsize=14)
axes[2][0].set_ylabel(r'$p\left(\mathbf{M}\right) (MPa \cdot mm^2)$', fontsize=14)
axes[2][0].legend(['index K = 1', 'index K = 2', 'index K = 3'])
axes[2][1].set_ylabel(r'$q\left(\mathbf{M}\right) (MPa \cdot mm^2)$', fontsize=14)
axes[2][1].legend(['index K = 1', 'index K = 2', 'index K = 3'])
axes[2][2].set_ylabel(r'$||dev\left(\Gamma_{IJK}\right)||$', fontsize=14)
axes[2][0].set_xlabel(r'Simulation time (s)', fontsize=14)
axes[2][1].set_xlabel(r'Simulation time (s)', fontsize=14)
axes[2][2].set_xlabel(r'Simulation time (s)', fontsize=14)
fig.set_figheight(12)
fig.set_figwidth(12)
fig.tight_layout()
fig.savefig(output_name)
return 0
[docs]
def p_q_plot(PK2, SIGMA, M, E, nqp, nel, ninc, output_name, dim=3):
'''Plot the pressure versus deviatoric norm (P-Q) plot for second Piola-Kirchhoff, symmetric micro, and higher order stresses
:param dict PK2: The quantities dict storing second Piola-Kirchhoff stress
:param dict SIGMA: The quantities dict storing symmetric micro stress
:param dict M: THe quantities dict storing higher order stress
:param dict E: The quantities dict storing Green-Lagrange strain
:param dict Ecal: The quantities dict storing micro-strain
:param dict Gamma: The quantities dict storing micro deformation gradient
:param int nqp: The number of quadrature points
:param int nel: The number of elements
:param int ninc: The number of time increments
:param str output_name: Output filename
:param int dim: The dimension of the arrays, '2' for 2d or '3' for 3d
:returns: ``output_name`` plot
'''
fig = matplotlib.pyplot.figure('p_q', figsize=(9,9), dpi=300)
ax1 = fig.add_subplot(2,2,1)
ax2 = fig.add_subplot(2,2,2)
ax3 = fig.add_subplot(2,2,3)
ax4 = fig.add_subplot(2,2,4)
for el in range(nel):
for qp in range(nqp):
# take norms
PK2_p, PK2_q, SIGMA_p, SIGMA_q =[], [], [], []
diff_p, diff_q = [], []
M1_p, M1_q, M2_p, M2_q, M3_p, M3_q = [], [], [], [], [], []
for t in range(ninc):
PK2_p.append((-1/dim)*numpy.trace(PK2[qp][t,el,:,:]))
PK2_q.append(numpy.sqrt(0.5*numpy.tensordot(deviatoric(PK2[qp][t,el,:,:], dim=dim), deviatoric(PK2[qp][t,el,:,:], dim=dim))))
SIGMA_p.append((-1/dim)*numpy.trace(SIGMA[qp][t,el,:,:]))
SIGMA_q.append(numpy.sqrt(0.5*numpy.tensordot(deviatoric(SIGMA[qp][t,el,:,:], dim=dim), deviatoric(SIGMA[qp][t,el,:,:], dim=dim))))
diff_p.append((-1/dim)*numpy.trace(PK2[qp][t,el,:,:]-SIGMA[qp][t,el,:,:]))
diff_q.append(numpy.sqrt(0.5*numpy.tensordot(deviatoric(PK2[qp][t,el,:,:]-SIGMA[qp][t,el,:,:], dim=dim),
deviatoric(PK2[qp][t,el,:,:]-SIGMA[qp][t,el,:,:], dim=dim))))
M1_p.append((-1/dim)*numpy.trace(M[qp][t,el,:,:,0]))
M1_q.append(numpy.sqrt(0.5*numpy.tensordot(deviatoric(M[qp][t,el,:,:,0], dim=dim), deviatoric(M[qp][t,el,:,:,0], dim=dim))))
M2_p.append((-1/dim)*numpy.trace(M[qp][t,el,:,:,1]))
M2_q.append(numpy.sqrt(0.5*numpy.tensordot(deviatoric(M[qp][t,el,:,:,1], dim=dim), deviatoric(M[qp][t,el,:,:,1], dim=dim))))
if dim == 3:
M3_p.append((-1/dim)*numpy.trace(M[qp][t,el,:,:,2]))
M3_q.append(numpy.sqrt(0.5*numpy.tensordot(deviatoric(M[qp][t,el,:,:,2], dim=dim), deviatoric(M[qp][t,el,:,:,1], dim=dim))))
label = f"qp #{(qp+1)+(nel*nqp)}"
ax1.plot(PK2_p, PK2_q, label=label)
ax1.set_xlabel(r'$p\left(\mathbf{S}\right) (MPa)$', fontsize=14)
ax1.set_ylabel(r'$q\left(\mathbf{S}\right) (MPa)$', fontsize=14)
ax2.plot(SIGMA_p, SIGMA_q, label=label)
ax2.set_xlabel(r'$p\left(\mathbf{\Sigma}\right) (MPa)$', fontsize=14)
ax2.set_ylabel(r'$q\left(\mathbf{\Sigma}\right) (MPa)$', fontsize=14)
ax3.plot(diff_p, diff_q, label=label)
ax3.set_xlabel(r'$p\left(\mathbf{S} - \mathbf{\Sigma}\right) (MPa)$', fontsize=14)
ax3.set_ylabel(r'$q\left(\mathbf{S} - \mathbf{\Sigma}\right) (MPa$)', fontsize=14)
ax4.plot(M1_p, M1_q, '-o', label=f"index K = 1, {label}")
ax4.plot(M2_p, M2_q, '-^', label=f"index K = 2, {label}")
if dim == 3:
ax4.plot(M3_p, M3_q, '-v', label=f"index K = 3, {label}")
ax4.set_xlabel(r'$p\left(\mathbf{M}\right) (MPa \cdot mm^2)$', fontsize=14)
ax4.set_ylabel(r'$q\left(\mathbf{M}\right) (MPa \cdot mm^2)$', fontsize=14)
ax4.legend(['index K = 1', 'index K = 2', 'index K = 3'])
fig.tight_layout()
fig.savefig(output_name)
return 0
[docs]
def csv_machine(quantity, output_file, nqp, nel, time, three_point=False, dim=3):
'''Write a text dump summarizing a variety of statistics
:param dict quantities: A 2nd or third order tensor dictionary with keys for quadrature points and values storing a 4d array where indices correspond to time, element number, component i, and component j
:param str output_file: Output filename
:param int nqp: The dictionary key of the quantityt to output
:param int nel: The desired element of the quantity to output
:param int time: The desired time increment to output the quantity
:param bool three_point: Parameter indicating if tensor is 2nd or 3rd order. "False" by default
:param int dim: The dimension of the arrays, '2' for 2d or '3' for 3d
:returns: ``output_file``
'''
comps, means, mins, maxs, devs, covs = [], [], [], [], [], []
if three_point == False:
for i in range(dim):
for j in range(dim):
comps.append(f"{i+1}{j+1}")
everything = numpy.array([quantity[qp][time,el,i,j] for qp in range(nqp) for el in range(nel)]).flatten()
means.append(numpy.mean(everything))
mins.append(numpy.min(everything))
maxs.append(numpy.max(everything))
devs.append(numpy.std(everything))
covs.append(numpy.std(everything)/numpy.mean(everything))
else:
for i in range(dim):
for j in range(dim):
for k in range(dim):
comps.append(f"{i+1}{j+1}{k+1}")
everything = numpy.array([quantity[qp][time,el,i,j,k] for qp in range(nqp) for el in range(nel)]).flatten()
means.append(numpy.mean(everything))
mins.append(numpy.min(everything))
maxs.append(numpy.max(everything))
devs.append(numpy.std(everything))
covs.append(numpy.std(everything)/numpy.mean(everything))
df = pandas.DataFrame({'component': comps,
'mean': means,
'min': mins,
'max': maxs,
'dev': devs,
'cov': covs})
print(f'building csv file for {output_file}')
df.to_csv(output_file, header=True, sep=',', index=False)
return 0
[docs]
def dump_all(quantity, output_file, icomp, jcomp, nqp, nel, time):
'''Dump all values for particular components of a 2nd order tensor for a given time value and element number
:param dict quantities: A 2nd order tensor dictionary with keys for quadrature points
and values storing a 4d array where indices correspond to time, element number,
component i, and component j
:param str output_file: Output filename
:param int icomp: the desired i component of the quantity to output
:param int jcomp: the desired j component of the quantity to output
:param int nqp: The number of quadrature points
:param int nel: The desired element of the quantity to output
:param int time: The desired time increment to output the quantity
:returns: ``output_file``
'''
everything = numpy.array([quantity[qp][time,el,icomp,jcomp] for qp in range(nqp) for el in range(nel)]).flatten()
df = pandas.DataFrame({f'quantity': everything})
df.to_csv(output_file, header=True, sep=',', index=False)
return 0
[docs]
def average_quantities(quantities, dim=3):
'''Average 2nd order tensor quantites over quadrature points
:param dict quantities: A 2nd order tensor dictionary with keys for quadrature points
and values storing a 4d array where indices correspond to time, element number,
component i, and component j
:param int dim: The dimension of the arrays, '2' for 2d or '3' for 3d
:returns: ``output`` dict with same indices as ``quantities`` and a single key
'''
output = {}
output[0] = numpy.zeros_like(quantities[0])
for i in range(dim):
for j in range(dim):
mean_field = []
for qp in quantities.keys():
mean_field.append(quantities[qp][:,0,i,j])
means = numpy.mean(mean_field, axis=0)
output[0][:,0,i,j] = means
return(output)
[docs]
def csv_all_quantities(output_file, times, nqp, quantities, quantity_names, dim=3):
'''Create csv file of a list of quantities for a single domain
:param str output_file: Output filename
:param array-like times: Array of times
:param int nqp: The number of quadrature points
:param list quantities: A list containing Micromorphic Filter output quantities
:param list quantity_names: A list containing the names of Micromorphic Filter output quantities
:param int dim: The dimension of the arrays, '2' for 2d or '3' for 3d
:returns: Write ``output_file``
'''
# Collect all results
output_dict = {}
output_dict['time'] = numpy.array([[time for qp in range(nqp) for time in times]]).flatten()
output_dict['qp'] = numpy.array([[qp for qp in range(nqp) for time in times]]).flatten()
for quantity, name in zip(quantities, quantity_names):
# volume fractions
if len(numpy.shape(quantity[0])) == 0:
key = f'{name}'
qp_data = []
for qp in range(nqp):
for t in range(len(times)):
qp_data.append(quantity[qp])
output_dict[key] = numpy.array(qp_data).flatten()
# vectors
elif len(numpy.shape(quantity[0])) == 3:
for i in range(dim):
key = f'{name}_{i+1}'
qp_data = []
for qp in range(nqp):
qp_data.append(quantity[qp][:,0,i])
output_dict[key] = numpy.array(qp_data).flatten()
# second order tensors
elif len(numpy.shape(quantity[0])) == 4:
for i in range(dim):
for j in range(dim):
key = f'{name}_{i+1}{j+1}'
qp_data = []
for qp in range(nqp):
qp_data.append(quantity[qp][:,0,i,j])
output_dict[key] = numpy.array(qp_data).flatten()
# third order tensors
elif len(numpy.shape(quantity[0])) == 5:
for i in range(dim):
for j in range(dim):
for k in range(dim):
key = f'{name}_{i+1}{j+1}{k+1}'
qp_data = []
for qp in range(nqp):
qp_data.append(quantity[qp][:,0,i,j,k])
output_dict[key] = numpy.array(qp_data).flatten()
else:
print(f'Data shape not supported! quantity = {quantity}')
# Output CSV
df = pandas.DataFrame(output_dict)
print(f'building csv file for {output_file}')
df.to_csv(output_file, header=True, sep=',', index=False)
return 0
[docs]
def calculate_initial_grain_volume_fraction(density, rho_binder, rho_grain, nqp):
'''Calculate the initial grain volume fraction given the homogenized density and two reference densities
:param dict density: Dictionary containing homogenized densities
:param str rho_binder: The density of the binder material, required if 'csv-all-quantities-single-domain' is specified
:param str rho_grain: The density of the grain material, required if 'csv-all-quantities-single-domain' is specified
:param int nqp: The number of quadrature points
:returns: Dictionary of grain volume fractions for each quadrature point
'''
fgs = {}
for qp in range(nqp):
fgs[qp] = (density[qp][0]-rho_binder)/(rho_grain-rho_binder)
return fgs
[docs]
def visualize_results(input_file, average, num_domains,
plot_cauchy_couple=None,
plot_cauchy_stress=None,
plot_PK2_stress=None,
plot_symm_stress=None,
plot_SIGMA_stress=None,
plot_stress_diff=None,
plot_body_couples=None,
plot_spin_inertias=None,
plot_spin_diff=None,
plot_rotation_diff=None,
plot_stretch_diff=None,
plot_stress_norms=None,
plot_better_stress_norms=None,
plot_best_stress_norms=None,
plot_norm_histories=None,
p_q_plots=None,
csv_cauchy=None,
csv_PK2=None,
csv_GLstrain=None,
csv_estrain=None,
csv_ref_mod=None,
csv_cur_mod=None,
csv_symm=None,
csv_stress_diff=None,
csv_m=None,
csv_M=None,
csv_stress33_all=None,
csv_all_quantities_single_domain=None,
rho_binder=None,
rho_grain=None,
dim=3):
''' Post-process Micromorphic Filter output
:param str input_file: The XDMF Micromorphic Filter results file
:param bool average: Average over quadrature points if True
:param int num_domains: The number of filter domains
:param str plot_cauchy_couple: Optional filename to plot Cauchy couple vs. simulation time
:param str plot_cauchy_stress: Optional filename to plot Cauchy stress vs. Eulerian strain
:param str plot_PK2_stress: Optional filename to plot PK2 stress vs. Green-Lagrange strain
:param str plot_symm_stress: Optional filename to plot symmetric micro stress vs. Eulerian strain
:param str plot_SIGMA_stress: Optional filename to plot Symmetric micro stress vs. Green-Lagrange strain
:param str plot_stress_diff: Optional filename to plot difference between Cauchy and symmetric micro stresses vs. simulation time
:param str plot_body_couples: Optional filename to plot body couples vs. simulation time
:param str plot_spin_inertias: Optional filename to plot micro spin inertias vs. simulation time
:param str plot_spin_diff: Optional filename to plot difference between body couples and micro spin inertias vs. simulation time
:param str plot_rotation_diff: Optional filename to plot difference between macro and micro rotations vs. simulation time
:param str plot_stress_diff: Optional filename to plot differences between macro and micro stretches vs. simulation time
:param str plot_stress_norms: Optional filename to plot norms of cauchy stress, symmetric micro stress, difference between Cauchy and symmetric micro stresses, and higher order stress
:param str plot_better_stress_norms: Optional filename to plot norms of PK2 stress, Symmetric micro stress, difference between PK2 and Symmetric micro stresses, and higher order stress, all against norms of Green-Lagrange strain, Micro strain, and Micro-deformation
:param str plot_best_stress_norms: Optional filename to plot norm of PK2 stress vs Green-Lagrange strain, Symmetric micro stress vs Micro strain, and higher order stress vs micro-deformation
:param str csv_cauchy: Optional filename for csv output of Cauchy stress summary statistics
:param str csv_PK2: Optional filename for csv output of PK2 stress summary statistics
:param str csv_GLstrain: Optional filename for csv output of Green-Lagrange strain summary statistics
:param str csv_ref_mod: Optional filename for csv output of 'moduli' calculation (S_{ij} / E_{ij}) in reference configuration summary statistics
:param str csv_cur_mod: Optional filename for csv output of 'moduli' calculation (\sigma_{ij} / e_{ij}) in the current configuration summary statistics
:param str csv_strain: Optional filename for csv output of Eulerian strain summary statistics
:param str csv_symm: Optional filename for csv output of symmetric micro stress summary statistics
:param str csv_stress_diff: Optional filename for csv output of difference between Cauchy and symmetric micro stresses summary statistics
:param str csv_m: Optional filename for csv output of couple stress (current configuration) summary statistics
:param str csv_M: Optional filename for csv output of couple stress (reference configuration) summary statistics
:param str csv_stress33_all: Optional filename for csv output of all Cauchy 33 values
:param str csv_all_quantities_single_domain: Optional filename for csv output of all quantities for a single domain
:param str rho_binder: The density of the binder material, required if 'csv-all-quantities-single-domain' is specified
:param str rho_grain: The density of the grain material, required if 'csv-all-quantities-single-domain' is specified
:param int dim: The dimension of the arrays, '2' for 2d or '3' for 3d
'''
# Read in data
data, geometry, topology = XRT.parse_xdmf_output(input_file)
## TODO: generalize for any element type
if dim == 3:
nqp = 8
elif dim == 2:
nqp = 4
ninc = numpy.shape(data['time'])[0]
nel = num_domains
# Read in the DOF information
displacement, gradu, phi, gradphi = XRT.construct_degrees_of_freedom(data, nqp, nel, dim=dim)
# Read in the stress information
cauchy, symm, m = XRT.collect_stresses(data, nqp, nel, dim=dim)
PK2, SIGMA, M = XRT.get_reference_configuration_stresses(data, nqp, nel, dim=dim)
# Read in the strain information
E, Ecal, Gamma, F, chi, grad_chi, e, h = XRT.compute_deformations(data, nqp, nel, dim=dim)
# Read in first moment of momentum measures
coup, spins = XRT.collect_first_moment_of_momentum_measures(data, nqp, nel, dim=dim)
# Get times
times = numpy.unique(data['time'])
# average fields if selected as an option
if average:
cauchy = average_quantities(cauchy, dim)
symm = average_quantities(symm, dim)
E = average_quantities(E, dim)
PK2 = average_quantities(PK2, dim)
SIGMA = average_quantities(SIGMA, dim)
F = average_quantities(F, dim)
chi = average_quantities(chi, dim)
e = average_quantities(e, dim)
h = average_quantities(h, dim)
coup = average_quantities(coup, dim)
spins = average_quantities(spins, dim)
# Plot stresses
## plot cauchy and PK2 stress
if plot_cauchy_stress:
plot_current_stresses(e, cauchy, average, output_name=plot_cauchy_stress, dim=dim)
if plot_PK2_stress:
plot_reference_stresses(E, PK2, average, output_name=plot_PK2_stress, dim=dim)
## plot symmetric micro-stress, s and SIGMA
if plot_symm_stress:
plot_current_stresses(e, symm, average, output_name=plot_symm_stress, dim=dim)
if plot_SIGMA_stress:
plot_reference_stresses(E, SIGMA, average, output_name=plot_SIGMA_stress, dim=dim)
## Check #1: Symmetry of Cauchy Stress using Axial Vector
if plot_cauchy_couple:
plot_axial_vector_mag(cauchy, times, average, output_name=plot_cauchy_couple, dim=dim)
## Check #2: Difference between Cauchy and Symmetric Micro-Stresses
if plot_stress_diff:
plot_stress_diffs(cauchy, symm, times, average, output_name=plot_stress_diff, dim=dim)
## Check #3: Micromorphic Effects - Difference between macro- and micro-deformation gradients
# Get R and U
if plot_rotation_diff or plot_stretch_diff:
R, U, Rchi, Uchi = get_R_and_U(PK2, F, chi, times, dim=dim)
# plot difference between rotations
if plot_rotation_diff:
plot_rot_diffs(PK2, times, R, Rchi, average, output_name=plot_rotation_diff, dim=dim)
# plot difference between stretches
if plot_stretch_diff:
plot_stretch_diffs(PK2, times, U, Uchi, average, output_name=plot_stretch_diff, dim=dim)
# Check #4: Different between body_couple and micro_spin_inertias, plot them too
if plot_body_couples or plot_spin_inertias or plot_spin_diff:
plot_first_moment_of_momentum_measures(coup, spins, times, average, plot_body_couples, plot_spin_inertias, plot_spin_diff, dim=dim)
# Plot stress norms
if plot_stress_norms:
plot_stress_norm(cauchy, symm, m, nqp, nel, ninc, times, plot_stress_norms, dim=dim)
# Plot better and best stress norms
if plot_better_stress_norms:
plot_better_stress_norm(PK2, SIGMA, M, E, Ecal, Gamma, nqp, nel, ninc, plot_better_stress_norms, dim=dim)
if plot_best_stress_norms:
plot_best_stress_norm(PK2, SIGMA, M, E, Ecal, Gamma, nqp, nel, ninc, plot_best_stress_norms, dim=dim)
transparent_name = f"{plot_best_stress_norms.split('.')[0]}_transparent.{plot_best_stress_norms.split('.')[1]}"
plot_best_stress_norm(PK2, SIGMA, M, E, Ecal, Gamma, nqp, nel, ninc, transparent_name, 0.05, dim=dim)
filtered_name = f"{plot_best_stress_norms.split('.')[0]}_filtered.{plot_best_stress_norms.split('.')[1]}"
plot_best_stress_norm(PK2, SIGMA, M, E, Ecal, Gamma, nqp, nel, ninc, filtered_name, upper_limit=5000., dim=dim)
if plot_norm_histories:
plot_norm_history(PK2, SIGMA, M, E, Ecal, Gamma, nqp, nel, ninc, times, plot_norm_histories, dim=dim)
# p_q_plots
if p_q_plots:
p_q_plot(PK2, SIGMA, M, E, nqp, nel, ninc, p_q_plots, dim=dim)
# Output csvs
#TODO: add an argument for specifying what time to gather statistics on
if csv_cauchy:
csv_machine(cauchy, csv_cauchy, nqp, nel, ninc-1, dim=dim)
if csv_PK2:
csv_machine(PK2, csv_PK2, nqp, nel, ninc-1, dim=dim)
if csv_GLstrain:
csv_machine(E, csv_GLstrain, nqp, nel, ninc-1, dim=dim)
if csv_estrain:
csv_machine(e, csv_estrain, nqp, nel, ninc-1, dim=dim)
if csv_ref_mod:
mod = {}
for qp in range(nqp):
values = numpy.zeros((ninc, nel, dim, dim))
for el in range(nel):
for t in range(ninc):
print(f'shape of PK2[qp][t,el,:,:] = {numpy.shape(PK2[qp][t,el,:,:])}')
print(f'shape of E[qp][t,el,:,:] = {numpy.shape(E[qp][t,el,:,:])}')
values[t,el,:,:] = PK2[qp][t,el,:,:] / E[qp][t,el,:,:]
mod.update({qp:numpy.copy(values)})
csv_machine(mod, csv_ref_mod, nqp, nel, ninc-1, dim=dim)
if csv_cur_mod:
mod = {}
for qp in range(nqp):
values = numpy.zeros((ninc, nel, dim, dim))
for el in range(nel):
for t in range(ninc):
values[t,el,:,:] = cauchy[qp][t,el,:,:] / e[qp][t,el,:,:]
mod.update({qp:numpy.copy(values)})
csv_machine(mod, csv_cur_mod, nqp, nel, ninc-1, dim=dim)
if csv_symm:
csv_machine(symm, csv_symm, nqp, nel, ninc-1, dim=dim)
# TODO: make this better lol
if csv_stress_diff:
diff = {}
for qp in range(nqp):
values = numpy.zeros((ninc, nel, dim, dim))
for el in range(nel):
for t in range(ninc):
values[t,el,:,:] = cauchy[qp][t,el,:,:] - symm[qp][t,el,:,:]
diff.update({qp:numpy.copy(values)})
csv_machine(diff, csv_stress_diff, nqp, nel, ninc-1, dim=dim)
# Higher order stresses
if csv_m:
csv_machine(m, csv_m, nqp, nel, ninc-1, three_point=True, dim=dim)
if csv_M:
csv_machine(M, csv_M, nqp, nel, ninc-1, three_point=True, dim=dim)
# Dump all cauchy
if csv_stress33_all:
dump_all(cauchy, csv_stress33_all, dim-1, dim-1, nqp, nel, ninc-1)
# Summary csv for quantities of a single single domains
if csv_all_quantities_single_domain:
# get initial volume fractions
assert rho_binder
assert rho_grain
fgs = calculate_initial_grain_volume_fraction(data['density_0'], rho_binder, rho_grain, nqp)
print(f'fgs = {fgs}')
quantities = [fgs, displacement, gradu, phi, gradphi, F, chi, E, Ecal, Gamma, e, h, cauchy, symm, m]
quantities_names = ['f_g', 'u', 'gradu', 'phi', 'gradphi', 'F', 'chi', 'E', 'Ecal', 'Gamma', 'e', 'h', 'cauchy', 'symm', 'm']
csv_all_quantities(csv_all_quantities_single_domain, times, nqp, quantities, quantities_names, dim=dim)
return 0
def get_parser():
script_name = pathlib.Path(__file__)
prog = f"python {script_name.name} "
cli_description = "Post-process Micromorphic Filter Output"
parser = argparse.ArgumentParser(description=cli_description, prog=prog)
parser.add_argument('-i', '--input-file', type=str, required=True,
help="The XDMF Micromorphic Filter results file")
parser.add_argument('--average', type=str, required=False, default=False,
help='Boolean whether or not homogenized DNS results will be averaged')
parser.add_argument('--num-domains', type=int, required=False, default=1,
help="Specify the number of filter domains")
parser.add_argument('--plot-cauchy-couple', type=str, required=False, default=None,
help="Optional filename to plot Cauchy couple vs. simulation time")
parser.add_argument('--plot-cauchy-stress', type=str, required=False, default=None,
help="Optional filename to plot Cauchy stress vs. Eulerian strain")
parser.add_argument('--plot-PK2-stress', type=str, required=False, default=None,
help="Optional filename to plot PK2 stress vs. Green-Lagrange strain")
parser.add_argument('--plot-symm-stress', type=str, required=False, default=None,
help="Optional filename to plot symmetric micro stress vs. Eulerian strain")
parser.add_argument('--plot-SIGMA-stress', type=str, required=False, default=None,
help="Optional filename to plot Symmetric micro stress vs. Green-Lagrange strain")
parser.add_argument('--plot-stress-diff', type=str, required=False, default=None,
help="Optional filename to plot difference between Cauchy and symmetric micro\
stresses vs. simulation time")
parser.add_argument('--plot-body-couples', type=str, required=False, default=None,
help="Optional filename to plot body couples vs. simulation time")
parser.add_argument('--plot-spin-inertias', type=str, required=False, default=None,
help="Optional filename to plot micro spin inertias vs. simulation time")
parser.add_argument('--plot-spin-diff', type=str, required=False, default=None,
help="Optional filename to plot difference between body couples and micro spin\
inertias vs. simulation time")
parser.add_argument('--plot-rotation-diff', type=str, required=False, default=None,
help="Optional filename to plot difference between macro and micro rotations\
vs. simulation time")
parser.add_argument('--plot-stretch-diff', type=str, required=False, default=None,
help="Optional filename to plot differences between macro and micro stretches\
vs. simulation time")
parser.add_argument('--plot-stress-norms', type=str, required=False, default=None,
help="Optional filename to plot norms of cauchy stress, symmetric micro stress,\
difference between Cauchy and symmetric micro stresses, and higher order\
stress.")
parser.add_argument('--plot-better-stress-norms', type=str, required=False, default=None,
help="Optional filename to plot norms of PK2 stress, Symmetric micro stress,\
difference between PK2 and Symmetric micro stresses, and higher order\
stress, all against norms of Green-Lagrange strain, Micro strain,\
and Micro-deformation.")
parser.add_argument('--plot-best-stress-norms', type=str, required=False, default=None,
help="Optional filename to plot norm of PK2 stress vs Green-Lagrange strain,\
Symmetric micro stress vs Micro strain,\
and higher order stress vs micro-deformation")
parser.add_argument('--plot-norm-histories', type=str, required=False, default=None,
help="Optional filename to plot norm of p and q invariants of Pk2 stress,\
Symmetric micro stress, and higher order stress, and the norms\
of Green-Lagrange strain, Micro strain, and Micro-deformation.\
All plots have an x-axis of pseudo-time")
parser.add_argument('--p-q-plots', type=str, required=False, default=None,
help="Optional filename to plot p-q invariants of PK2 stress, Symmetric\
micro stress, difference btween PK2 and Symmtric micro stresses,\
and higher order stress.")
parser.add_argument('--csv-cauchy', type=str, required=False, default=None,
help="Optional filename for csv output of Cauchy stress summary statistics")
parser.add_argument('--csv-PK2', type=str, required=False, default=None,
help="Optional filename for csv output of PK2 stress summary statistics")
parser.add_argument('--csv-GLstrain', type=str, required=False, default=None,
help="Optional filename for csv output of Green-Lagrange strain summary statistics")
parser.add_argument('--csv-ref-mod', type=str, required=False, default=None,
help="Optional filename for csv output of 'moduli' calculation (S_{ij} / E_{ij})\
in reference configuration summary statistics")
parser.add_argument('--csv-cur-mod', type=str, required=False, default=None,
help="Optional filename for csv output of 'moduli' calculation (\sigma_{ij} / e_{ij})\
in the current configuration summary statistics")
parser.add_argument('--csv-estrain', type=str, required=False, default=None,
help="Optional filename for csv output of Eulerian strain summary statistics")
parser.add_argument('--csv-symm', type=str, required=False, default=None,
help="Optional filename for csv output of symmetric micro stress summary statistics")
parser.add_argument('--csv-stress-diff', type=str, required=False, default=None,
help="Optional filename for csv output of difference between Cauchy and symmetric\
micro stresses summary statistics")
parser.add_argument('--csv-m', type=str, required=False, default=None,
help="Optional filename for csv output of couple stress (current configuration)\
summary statistics")
parser.add_argument('--csv-M', type=str, required=False, default=None,
help="Optional filename for csv output of couple stress (reference configuration)\
summary statistics")
parser.add_argument('--csv-stress33-all', type=str, required=False, default=None,
help="Optional filename for csv output of all Cauchy 33 values")
parser.add_argument('--csv-all-quantities-single-domain', type=str, required=False, default=None,
help="Optional filename for csv output of all quantities for a single domain")
parser.add_argument('--rho-binder', type=float, required=False, default=None,
help="The density of the binder material, required if '--csv-all-quantities-single-domain' is specified")
parser.add_argument('--rho-grain', type=float, required=False, default=None,
help="The density of the grain material, required if '--csv-all-quantities-single-domain' is specified")
parser.add_argument('--dim', type=int, required=False, default=3,
help="The dimension of the arrays, '2' for 2d or '3' for 3d")
return parser
if __name__ == '__main__':
parser = get_parser()
args, unknown = parser.parse_known_args()
sys.exit(visualize_results(
input_file=args.input_file,
average=str2bool(args.average),
num_domains=args.num_domains,
plot_cauchy_couple=args.plot_cauchy_couple,
plot_cauchy_stress=args.plot_cauchy_stress,
plot_PK2_stress=args.plot_PK2_stress,
plot_symm_stress=args.plot_symm_stress,
plot_SIGMA_stress=args.plot_SIGMA_stress,
plot_stress_diff=args.plot_stress_diff,
plot_body_couples=args.plot_body_couples,
plot_spin_inertias=args.plot_spin_inertias,
plot_spin_diff=args.plot_spin_diff,
plot_rotation_diff=args.plot_rotation_diff,
plot_stretch_diff=args.plot_stretch_diff,
plot_stress_norms=args.plot_stress_norms,
plot_better_stress_norms=args.plot_better_stress_norms,
plot_best_stress_norms=args.plot_best_stress_norms,
plot_norm_histories=args.plot_norm_histories,
p_q_plots=args.p_q_plots,
csv_cauchy=args.csv_cauchy,
csv_PK2=args.csv_PK2,
csv_GLstrain=args.csv_GLstrain,
csv_estrain=args.csv_estrain,
csv_ref_mod=args.csv_ref_mod,
csv_cur_mod=args.csv_cur_mod,
csv_symm=args.csv_symm,
csv_stress_diff=args.csv_stress_diff,
csv_m=args.csv_m,
csv_M=args.csv_M,
csv_stress33_all=args.csv_stress33_all,
csv_all_quantities_single_domain=args.csv_all_quantities_single_domain,
rho_binder=args.rho_binder,
rho_grain=args.rho_grain,
dim=args.dim,
))