223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

mechanics done

This commit is contained in:
Samad Vakili 2020-06-09 21:27:08 +02:00
parent d365cc9e12
commit 694b7ec3c5
1 changed files with 29 additions and 45 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

74
python/damask/mechanics.py
View File

@ -8,17 +8,13 @@ def Cauchy(P,F):
Parameters Parameters
---------- ----------
F : numpy.ndarray of shape (:,3,3) or (3,3) F : numpy.ndarray of shape (...,3,3)
Deformation gradient. Deformation gradient.
P : numpy.ndarray of shape (:,3,3) or (3,3) P : numpy.ndarray of shape (...,3,3)
First Piola-Kirchhoff stress. First Piola-Kirchhoff stress.
""" """
if _np.shape(F) == _np.shape(P) == (3,3): sigma = _np.einsum('...,...ij,...kj->...ik',1.0/_np.linalg.det(F),P,F)
sigma = 1.0/_np.linalg.det(F) * _np.dot(P,F.T)
else:
#sigma = _np.einsum('i,ijk,ilk->ijl',1.0/_np.linalg.det(F),P,F)
sigma = _np.einsum('...,...ij,...kj->...ik',1.0/_np.linalg.det(F),P,F)
return symmetric(sigma) return symmetric(sigma)
@ -44,7 +40,7 @@ def eigenvalues(T_sym):
Parameters Parameters
---------- ----------
T_sym : numpy.ndarray of shape (:,3,3) or (3,3) T_sym : numpy.ndarray of shape (...,3,3)
Symmetric tensor of which the eigenvalues are computed. Symmetric tensor of which the eigenvalues are computed.
""" """
@ -59,7 +55,7 @@ def eigenvectors(T_sym,RHS=False):
Parameters Parameters
---------- ----------
T_sym : numpy.ndarray of shape (:,3,3) or (3,3) T_sym : numpy.ndarray of shape (...,3,3)
Symmetric tensor of which the eigenvectors are computed. Symmetric tensor of which the eigenvectors are computed.
RHS: bool, optional RHS: bool, optional
Enforce right-handed coordinate system. Default is False. Enforce right-handed coordinate system. Default is False.
@ -68,10 +64,7 @@ def eigenvectors(T_sym,RHS=False):
(u,v) = _np.linalg.eigh(symmetric(T_sym)) (u,v) = _np.linalg.eigh(symmetric(T_sym))
if RHS: if RHS:
if _np.shape(T_sym) == (3,3): v[_np.linalg.det(v) < 0.0,:,2] *= -1.0
if _np.linalg.det(v) < 0.0: v[:,2] *= -1.0
else:
v[_np.linalg.det(v) < 0.0,:,2] *= -1.0
return v return v
@ -81,7 +74,7 @@ def left_stretch(T):
Parameters Parameters
---------- ----------
T : numpy.ndarray of shape (:,3,3) or (3,3) T : numpy.ndarray of shape (...,3,3)
Tensor of which the left stretch is computed. Tensor of which the left stretch is computed.
""" """
@ -94,13 +87,12 @@ def maximum_shear(T_sym):
Parameters Parameters
---------- ----------
T_sym : numpy.ndarray of shape (:,3,3) or (3,3) T_sym : numpy.ndarray of shape (...,3,3)
Symmetric tensor of which the maximum shear is computed. Symmetric tensor of which the maximum shear is computed.
""" """
w = eigenvalues(T_sym) w = eigenvalues(T_sym)
return (w[0] - w[2])*0.5 if _np.shape(T_sym) == (3,3) else \ return (w[...,0] - w[...,2])*0.5
(w[...,0] - w[...,2])*0.5
def Mises_strain(epsilon): def Mises_strain(epsilon):
@ -109,7 +101,7 @@ def Mises_strain(epsilon):
Parameters Parameters
---------- ----------
epsilon : numpy.ndarray of shape (:,3,3) or (3,3) epsilon : numpy.ndarray of shape (...,3,3)
Symmetric strain tensor of which the von Mises equivalent is computed. Symmetric strain tensor of which the von Mises equivalent is computed.
""" """
@ -122,7 +114,7 @@ def Mises_stress(sigma):
Parameters Parameters
---------- ----------
sigma : numpy.ndarray of shape (:,3,3) or (3,3) sigma : numpy.ndarray of shape (...,3,3)
Symmetric stress tensor of which the von Mises equivalent is computed. Symmetric stress tensor of which the von Mises equivalent is computed.
""" """
@ -135,16 +127,13 @@ def PK2(P,F):
Parameters Parameters
---------- ----------
P : numpy.ndarray of shape (...,3,3) or (3,3) P : numpy.ndarray of shape (...,3,3)
First Piola-Kirchhoff stress. First Piola-Kirchhoff stress.
F : numpy.ndarray of shape (...,3,3) or (3,3) F : numpy.ndarray of shape (...,3,3)
Deformation gradient. Deformation gradient.
""" """
if _np.shape(F) == _np.shape(P) == (3,3): S = _np.einsum('...jk,...kl->...jl',_np.linalg.inv(F),P)
S = _np.dot(_np.linalg.inv(F),P)
else:
S = _np.einsum('...jk,...kl->...jl',_np.linalg.inv(F),P)
return symmetric(S) return symmetric(S)
@ -154,7 +143,7 @@ def right_stretch(T):
Parameters Parameters
---------- ----------
T : numpy.ndarray of shape (:,3,3) or (3,3) T : numpy.ndarray of shape (...,3,3)
Tensor of which the right stretch is computed. Tensor of which the right stretch is computed.
""" """
@ -167,7 +156,7 @@ def rotational_part(T):
Parameters Parameters
---------- ----------
T : numpy.ndarray of shape (:,3,3) or (3,3) T : numpy.ndarray of shape (...,3,3)
Tensor of which the rotational part is computed. Tensor of which the rotational part is computed.
""" """
@ -199,7 +188,7 @@ def strain_tensor(F,t,m):
Parameters Parameters
---------- ----------
F : numpy.ndarray of shape (:,3,3) or (3,3) F : numpy.ndarray of shape (...,3,3)
Deformation gradient. Deformation gradient.
t : {‘V’, ‘U’} t : {‘V’, ‘U’}
Type of the polar decomposition, ‘V’ for left stretch tensor and ‘U’ for right stretch tensor. Type of the polar decomposition, ‘V’ for left stretch tensor and ‘U’ for right stretch tensor.
@ -207,12 +196,11 @@ def strain_tensor(F,t,m):
Order of the strain. Order of the strain.
""" """
F_ = F.reshape(1,3,3) if F.shape == (3,3) else F
if t == 'V': if t == 'V':
B = _np.matmul(F_,transpose(F_)) B = _np.matmul(F,transpose(F))
w,n = _np.linalg.eigh(B) w,n = _np.linalg.eigh(B)
elif t == 'U': elif t == 'U':
C = _np.matmul(transpose(F_),F_) C = _np.matmul(transpose(F),F)
w,n = _np.linalg.eigh(C) w,n = _np.linalg.eigh(C)
if m > 0.0: if m > 0.0:
@ -225,8 +213,7 @@ def strain_tensor(F,t,m):
else: else:
eps = _np.matmul(n,_np.einsum('...j,...kj->...jk',0.5*_np.log(w),n)) eps = _np.matmul(n,_np.einsum('...j,...kj->...jk',0.5*_np.log(w),n))
return eps.reshape(3,3) if _np.shape(F) == (3,3) else \ return eps
eps
def symmetric(T): def symmetric(T):
@ -235,7 +222,7 @@ def symmetric(T):
Parameters Parameters
---------- ----------
T : numpy.ndarray of shape (...,3,3) or (3,3) T : numpy.ndarray of shape (...,3,3)
Tensor of which the symmetrized values are computed. Tensor of which the symmetrized values are computed.
""" """
@ -248,12 +235,11 @@ def transpose(T):
Parameters Parameters
---------- ----------
T : numpy.ndarray of shape (...,3,3) or (3,3) T : numpy.ndarray of shape (...,3,3)
Tensor of which the transpose is computed. Tensor of which the transpose is computed.
""" """
return T.T if _np.shape(T) == (3,3) else \ return _np.swapaxes(T,axis2=-2,axis1=-1)
_np.swapaxes(T,axis2=-2,axis1=-1)
def _polar_decomposition(T,requested): def _polar_decomposition(T,requested):
@ -262,7 +248,7 @@ def _polar_decomposition(T,requested):
Parameters Parameters
---------- ----------
T : numpy.ndarray of shape (:,3,3) or (3,3) T : numpy.ndarray of shape (...,3,3)
Tensor of which the singular values are computed. Tensor of which the singular values are computed.
requested : iterable of str requested : iterable of str
Requested outputs: ‘R’ for the rotation tensor, Requested outputs: ‘R’ for the rotation tensor,
@ -270,16 +256,15 @@ def _polar_decomposition(T,requested):
""" """
u, s, vh = _np.linalg.svd(T) u, s, vh = _np.linalg.svd(T)
R = _np.dot(u,vh) if _np.shape(T) == (3,3) else \ R = _np.einsum('...ij,...jk->...ik',u,vh)
_np.einsum('...ij,...jk->...ik',u,vh)
output = [] output = []
if 'R' in requested: if 'R' in requested:
output.append(R) output.append(R)
if 'V' in requested: if 'V' in requested:
output.append(_np.dot(T,R.T) if _np.shape(T) == (3,3) else _np.einsum('...ij,...kj->...ik',T,R)) output.append(_np.einsum('...ij,...kj->...ik',T,R))
if 'U' in requested: if 'U' in requested:
output.append(_np.dot(R.T,T) if _np.shape(T) == (3,3) else _np.einsum('...ji,...jk->...ik',R,T)) output.append(_np.einsum('...ji,...jk->...ik',R,T))
return tuple(output) return tuple(output)
@ -290,12 +275,11 @@ def _Mises(T_sym,s):
Parameters Parameters
---------- ----------
T_sym : numpy.ndarray of shape (:,3,3) or (3,3) T_sym : numpy.ndarray of shape (...,3,3)
Symmetric tensor of which the von Mises equivalent is computed. Symmetric tensor of which the von Mises equivalent is computed.
s : float s : float
Scaling factor (2/3 for strain, 3/2 for stress). Scaling factor (2/3 for strain, 3/2 for stress).
""" """
d = deviatoric_part(T_sym) d = deviatoric_part(T_sym)
return _np.sqrt(s*(_np.sum(d**2.0))) if _np.shape(T_sym) == (3,3) else \ return _np.sqrt(s*_np.einsum('...jk->...',d**2.0))
_np.sqrt(s*_np.einsum('...jk->...',d**2.0))