# Module pylabfea.basic
'''Module pylabfea.basic introduces basic methods and attributes like calculation of
equivalent stresses and strain, conversions of
cylindrical to principal stresses and vice versa, and the classes ``Stress`` and ``Strain``
for efficient operations with these quantities.
uses NumPy, pickle
Version: 4.0 (2021-11-27)
Author: Alexander Hartmaier, ICAMS/Ruhr University Bochum, Germany
Email: alexander.hartmaier@rub.de
distributed under GNU General Public License (GPLv3)'''
import numpy as np
import sys
import pickle
#===================================
#define global methods and variables
#===================================
a_vec = np.array([1., -0.5, -0.5])/np.sqrt(1.5)
'''First unit vector spanning deviatoric stress plane (real axis)'''
b_vec = np.array([0., 0.5, -0.5])*np.sqrt(2)
'''Second unit vector spanning deviatoric stress plane (imaginary axis)'''
ptol = 5.e-3
'''Tolerance: Plastic yielding if yield function > ptol'''
[docs]def sig_eq_j2(sig):
'''Calculate J2 equivalent stress from any stress tensor
Parameters
----------
sig : (3,), (6,) (N,3) or (N,6) array
(3,), (N,3): Principal stress or list of principal stresses;
(6,), (N,6): Voigt stress
Returns
-------
seq : float or (N,) array
J2 equivalent stresses
'''
N = len(sig)
sh = np.shape(sig)
if sh==(3,):
N = 1 # sig is single principle stress vector
sp = np.array([sig])
elif sh==(6,):
N = 1 # sig is single Voigt stress
sp = np.array([sig_princ(sig)[0]])
elif sh==(N,6):
sp = sig_princ(sig)[0]
elif sh==(N,3):
sp = np.array(sig)
else:
print('*** sig_eq_j2: N, sh', N, sh, sys._getframe().f_back.f_code.co_name)
raise TypeError('Error: Unknown format of stress in sig_eq_j2')
d12 = sp[:,0] - sp[:,1]
d23 = sp[:,1] - sp[:,2]
d31 = sp[:,2] - sp[:,0]
J2 = 0.5*(np.square(d12) + np.square(d23) + np.square(d31))
seq = np.sqrt(J2) # J2 eqiv. stress
if sh==(3,) or sh==(6,):
seq = seq[0]
return seq
[docs]def sig_polar_ang(sig):
'''Transform stresses into polar angle on deviatoric plane spanned by a_vec and b_vec
Parameters
----------
sig : (3,), (6,) (N,3) or (N,6) array
(3,), (N,3): Principal stresses;
(6,), (N,6): Voigt stress
Returns
-------
theta : float or (N,) array
polar angles in deviatoric plane as positive angle between sig_princ and a_vec in range [-pi,+p]
'''
sh = np.shape(sig)
N = len(sig)
if sh==(3,):
N = 1
sp = np.array([sig])
elif sh==(6,):
sp = [sig_princ(sig)[0]]
elif sh==(N,6):
sp = sig_princ(sig)[0]
elif sh==(N,3):
sp = np.array(sig)
else:
print('*** polar_angle: N, sh', N, sh, sys._getframe().f_back.f_code.co_name)
raise TypeError('Error: Unknown format of stress in polar_angle')
hyd = np.sum(sp,axis=1)/3. # hydrostatic component
dev = sp - hyd[:,None] # deviatoric princ. stress
vn = np.linalg.norm(dev,axis=1) #norms of princ. stress vectors
ind = np.nonzero(vn<1.e-4)[0]
vn[ind] = 1.
dsa = np.dot(dev/vn[:,None],a_vec)
dsb = np.dot(dev/vn[:,None],b_vec)
theta = np.angle(dsa + 1j*dsb)
if sh==(3,) or sh==(6,):
theta=theta[0]
return theta
[docs]def sig_princ(sig):
'''Convert Voigt stress tensors into principal stresses and eigenvectors.
Parameters
----------
sig : (6,), (N,6), (3,3), or (N,3,3) array
Voigt stress tensor (dim=6) or Cartesian stress tensor (dim=3x3)
Returns
-------
spa : (3,) or (N,3) array
Principal stresses
eva : (3,3) or (N,3,3) array
Eigenvectors/rotation matrices of stress tensor
'''
N = len(sig)
sh = np.shape(sig)
if sh==(3,3):
N = 1 # sig is Cartesian single stress tensor
st=np.array([sig])
elif sh==(N,3,3):
st=np.array(sig) # sig is array of Cart. stress tensors
elif sh==(6,):
N = 1 # sig is single Voigt stress
st = np.zeros((1,3,3))
st[0,0,0]=sig[0]
st[0,1,1]=sig[1]
st[0,2,2]=sig[2]
st[0,2,1]=st[0,1,2]=sig[3]
st[0,2,0]=st[0,0,2]=sig[4]
st[0,1,0]=st[0,0,1]=sig[5]
elif sh==(N,6):
st=np.zeros((N,3,3))
for i in range(N):
st[i,0,0]=sig[i,0]
st[i,1,1]=sig[i,1]
st[i,2,2]=sig[i,2]
st[i,2,1]=st[i,1,2]=sig[i,3]
st[i,2,0]=st[i,0,2]=sig[i,4]
st[i,1,0]=st[i,0,1]=sig[i,5]
else:
print('*** sig_princ: N, sh', N, sh, sys._getframe().f_back.f_code.co_name)
raise TypeError('Error: Unknown format of stress in sig_princ')
#calculate principal stresses and eigen vectors
spa = np.zeros((N,3))
eva = np.zeros((N,3,3))
for n in range(N):
sp, ev = np.linalg.eig(st[n]) # solve eigenvalue problem
#arrange principal stress components according to major force axes
iev = np.argmax(np.abs(ev),axis=1)
j = np.zeros(3, dtype=int)
i0 = [i for i,x in enumerate(iev) if x == 0] # positions of indices 0
i1 = [i for i,x in enumerate(iev) if x == 1]
i2 = [i for i,x in enumerate(iev) if x == 2]
k0 = len(i0)
for i in range(k0):
j[i] = i0[i]
for i in range(len(i1)):
j[k0+i] = i1[i]
k0 += len(i1)
for i in range(len(i2)):
j[k0+i] = i2[i]
ev = np.array([ev[j[0],:], ev[j[1],:], ev[j[2],:]])
sp = np.array([sp[j[0]], sp[j[1]], sp[j[2]]])
#ensure positive determinant
if np.linalg.det(ev)<0:
ev *= -1
spa[n,:] = sp
eva[n,:,:] = ev
if sh==(3,3) or sh==(6,):
spa = spa[0]
eva = eva[0,:,:]
return spa, eva
[docs]def sig_cyl2princ(scyl):
'''Convert cylindrical stress into 3D Cartesian principle stress
Parameters
----------
scyl : (2,), (3,), (N,2) or (N,3) array
Cylindrical stress in form (seq, theta, (optional: p))
Returns
-------
sig_princ : (3,) or (N,3) array
principle deviatoric stresses
'''
sh = np.shape(scyl)
if sh==(2,) or sh==(3,):
scyl = np.array([scyl])
seq = scyl[:,0]
theta = scyl[:,1]
sig_princ = (np.tensordot(np.cos(theta), a_vec, axes=0) +
np.tensordot(np.sin(theta), b_vec, axes=0)) \
*np.sqrt(2./3.)*np.array([seq,seq,seq]).T
if sh[0]==3:
p = scyl[:,2]
sig_princ += np.array([p,p,p]).T/3.
if sh==(2,) or sh==(3,):
sig_princ=sig_princ[0]
return sig_princ
[docs]def sig_cyl2voigt(scyl, evec):
'''Convert cylindrical stress and eigenvectors into Voigt stress tensor
Parameters
----------
scyl : (3,) or (N,3) array
Cylindrical stress in form (seq, theta, p)
evec : (3,3) or (N,3,3) array
Eigenvectors of stress tensor
Returns
-------
sig_cyl2voigt : (6,) or (N,6) array
Voigt stress tensor
'''
#sh = np.shape(scyl)
#if sh==(3,):
# scyl = np.array([scyl])
sp = sig_cyl2princ(scyl)
if np.linalg.det(evec)<0:
evec *= -1 # enforce right-handed system of eigenvectors
st = np.diag(sp) # diag. matrix of princ. stresses
hh = evec@st@evec.T # rotate back into original stress frame
sig_cyl2voigt = np.array([hh[0,0], hh[1,1], hh[2,2], hh[1,2], hh[0,2], hh[0,1]])
#if sh==(3,):
# sig_cyl2voigt=sig_cyl2voigt[0]
return sig_cyl2voigt
[docs]def sig_princ2cyl(sig, mat=None):
'''convert principal stress into cylindrical stress vector
Parameters
----------
sig : (3,), (6,), (N,3) or (N,6) array
stress to be converted, if (3,) or (N,3) principal stress is assumed
mat : object of class ``Material``
Material for Hill-type principal stress (optional)
Returns
-------
sc : (3,) or (N,3) array
stress in cylindrical coordinates (seq, theta, p)
'''
sh = np.shape(sig)
N = len(sig)
if sh==(3,):
N = 1 # sig is single principle stress vector
sp = np.array([sig])
sig = np.array([[sig[0], sig[1], sig[2], 0, 0,0]])
elif sh==(N,3):
sp = np.array(sig)
sig = np.append(sig,np.zeros((N,3)), axis=1)
elif sh==(6,):
N = 1
sp = np.array([Stress(sig).princ])
sig = np.array([sig])
elif sh==(N,6):
sp = sig_princ(sig)[0]
else:
print('*** sig_princ2cyl (N,sh): ',N,sh,sys._getframe().f_back.f_code.co_name)
raise TypeError('Error in sig_princ2cyl: Format not supported')
sc = np.zeros((N,3))
if mat is None:
sc[:,0] = sig_eq_j2(sp)
else:
sc[:,0] = mat.calc_seq(sig)
sc[:,1] = sig_polar_ang(sp)
sc[:,2] = np.sum(sp, axis=1)/3.
if sh==(3,) or sh==(6,):
sc=sc[0]
return sc
[docs]def sig_dev(sig):
'''Calculate deviatoric stress component from given stress tensor
Parameters
----------
sig : (3,), (6,) (N,3) or (N,6) array
Returns
-------
sd : float or (N,) array
deviatoric stresses
'''
sh = np.shape(sig)
hyd = np.zeros(sh)
if sh==(3,) or sh==(6,):
p = np.sum(sig[0:3])/3.
hyd[0:3] = p[None]
else:
p = np.sum(sig[:,0:3],axis=1)/3.
hyd[:,0:3] = p[:,None]
sd = sig - hyd
return sd
[docs]def eps_eq(eps):
'''Calculate equivalent strain
Parameters
----------
eps : (3,), (6,), (N,3) or (N,6) array
(3,) or (N,3): Principal strains;
(6,) or (N,6): Voigt strains
Returns
-------
eeq : float or (N,) array
equivalent strains
'''
sh = np.shape(eps)
if sh==(6,) or sh==(3,):
eps = np.array([eps])
N = 1
else:
N = len(eps)
#ev = np.sum(eps[:,0:3],axis=1)
#ed = eps[:,0:3] - np.array([ev,ev,ev]).T
if sh==(6,) or sh==(N,6):
eeq = np.sqrt(2.*(np.sum(eps[:,0:3]*eps[:,0:3],axis=1)+0.5*np.sum(eps[:,3:6]*eps[:,3:6],axis=1))/3.)
elif sh==(3,) or sh==(N,3):
eeq = np.sqrt(2.*np.sum(eps[:,0:3]*eps[:,0:3],axis=1)/3.)
else:
print('*** eps_eq (N,sh): ',N,sh,sys._getframe().f_back.f_code.co_name)
raise ValueError('Error in eps_eq: Format not supported')
if sh==(6,) or sh==(3,):
eeq = eeq[0]
return eeq
#=========================
#define class for stresses
#=========================
[docs]class Stress(object):
'''Stores and converts Voigt stress tensors into different formats,
calculates principle stresses, equivalent stresses and transforms into cylindrical coordinates.
Parameters
----------
sv : list-like object, must be 1D with length 6
Voigt-stress components
Attributes
----------
voigt, v : 1d-array (size 6)
Stress tensor in Voigt notation
tens, t : 3x3 array
Stress tensor in matrix notation
princ, p : 1d-array (size 3)
Principal stresses
hydrostatic, h : float
Hydrostatic stress component
'''
def __init__(self, sv):
self.v=self.voigt = np.array(sv)
#calculate (3x3)-tensorial representation
self.t=self.tens = np.zeros((3,3))
self.tens[0,0]=sv[0]
self.tens[1,1]=sv[1]
self.tens[2,2]=sv[2]
self.tens[2,1]=self.tens[1,2]=sv[3]
self.tens[2,0]=self.tens[0,2]=sv[4]
self.tens[1,0]=self.tens[0,1]=sv[5]
#calcualte principal stresses and eigen vectors
self.princ, self.evec = sig_princ(self.tens)
self.p = self.princ
self.h = self.hydrostatic = np.sum(self.p)/3.
self.d = self.dev = sv - np.array([self.h, self.h, self.h, 0., 0., 0.])
[docs] def seq(self, mat=None):
'''calculate Hill-type equivalent stress, invokes corresponding method of class ``Material``
Parameters
----------
mat: object of class ``Material``
containes Hill parameters and method needed for Hill-type equivalent stress (optional, defaul=None)
Returns
-------
seq : float
equivalent stress; material specific (J2, Hill, Tresca, Barlat) if mat is provided,
J2 equivalent stress otherwise
'''
if mat is None:
seq = sig_eq_j2(self.p) # give princ. stress as parameter to avoid re-calculation
else:
seq = mat.calc_seq(self.v)
return seq
[docs] def theta(self):
'''calculate polar angle in deviatoric plane
Returns
-------
ang : float
polar angle of stress in devitoric plane
'''
ang = sig_polar_ang(self.p)
return ang
[docs] def seq_j2(self):
'''calculate J2 principal stress
Returns
-------
seq_j2 : float
equivalent stress
'''
seq_j2 = sig_eq_j2(self.p)
return seq_j2
[docs] def cyl(self):
'''Calculate cylindrical stress tensor
Returns
-------
cyl : (3,) array
stress in cylindrical form: (J2 eqiv. stress, polar angle, hydrostatic)
'''
cyl = np.array([sig_eq_j2(self.p), sig_polar_ang(self.p), self.h])
return cyl
[docs] def lode_ang(self, X):
'''Calculate Lode angle:
Transforms principal stress space into hydrostatic stress, eqiv. stress, and Lode angle;
definition of positive cosine for Lode angle is applied
Parameters
----------
X : either float or object of class ``Material``
if float: interpreted as equivalent stress
if ``Material``: used to invoke method of class ``Material`` to calculate equivalent stress
Returns
-------
la : float
Lode angle
'''
if type(X) is float:
seq = X # float-type parameters are interpreted as equiv. stress
else:
seq = self.seq(X) # otherwise parameter is Material
J3 = np.linalg.det(self.tens - self.h*np.diag(np.ones(3)))
hh = 0.5*J3*(3./seq)**3
la = np.arccos(hh)/3.
return la
#=======================
#define class for strain
#=======================
[docs]class Strain(object):
'''Stores and converts Voigt strain tensors into different formats,
calculates principle strain and equivalent strain.
Parameters
----------
sv : list-like object, must be 1D with length 6
Voigt-strain components
Attributes
----------
voigt, v : 1d-array (size 6)
Strain tensor in Voigt notation
tens, t : 3x3 array
Strain tensor in matrix notation
princ, p : 1d-array (size 3)
Principal strains
'''
def __init__(self, sv):
self.v=self.voigt = np.array(sv)
#calculate (3x3)-tensorial representation
self.t=self.tens = np.zeros((3,3))
self.tens[0,0]=sv[0]
self.tens[1,1]=sv[1]
self.tens[2,2]=sv[2]
self.tens[2,1]=self.tens[1,2]=sv[3]
self.tens[2,0]=self.tens[0,2]=sv[4]
self.tens[1,0]=self.tens[0,1]=sv[5]
#calcualte principal stresses and eigen vectors
self.princ, self.evec = np.linalg.eig(self.tens)
self.p=self.princ
[docs] def eeq(self):
'''Calculate equivalent strain
Returns
-------
eeq : float
Equivalent strain
'''
eeq=eps_eq(self.v)
return eeq
[docs] def inv(self):
'''Calculate inverse of strain tensor ignoring zeros.
Returns
-------
inv : (6,) array
'''
inv = np.zeros(6)
for i in range(6):
if np.abs(self.voigt[i])>1.e-9 : inv[i] = 1./self.voigt[i]
return inv
#=================================================
#define subroutines for reading objects from files
#=================================================
[docs]def pickle2mat(name, path='./'):
'''Read pickled material file.
Parameters
----------
name : string
File name of pickled material to be read.
path : string
Path under which pckl-files are stored (optional, default: './')
Returns
-------
pckl : Material object
Unpickled material
'''
if name is None:
raise ValueError('Name for pickled material must be given.')
if path[-1] != '/':
path += '/'
with open(path+name, 'rb') as input:
pckl = pickle.load(input)
return pckl