4 space indent is python standard
This commit is contained in:
Martin Diehl
parent
ffb112b0d8
commit
8a85123abc
|
|
@ -1,247 +1,247 @@
|
||||||
import numpy as np
|
import numpy as np
|
||||||
|
|
||||||
def Cauchy(F,P):
|
def Cauchy(F,P):
|
||||||
"""
|
"""
|
||||||
Return Cauchy stress calculated from 1. Piola-Kirchhoff stress and deformation gradient.
|
Return Cauchy stress calculated from 1. Piola-Kirchhoff stress and deformation gradient.
|
||||||
|
|
||||||
Resulting tensor is symmetrized as the Cauchy stress needs to be symmetric.
|
Resulting tensor is symmetrized as the Cauchy stress needs to be symmetric.
|
||||||
|
|
||||||
Parameters
|
Parameters
|
||||||
----------
|
----------
|
||||||
F : numpy.array of shape (:,3,3) or (3,3)
|
F : numpy.array of shape (:,3,3) or (3,3)
|
||||||
Deformation gradient.
|
Deformation gradient.
|
||||||
P : numpy.array of shape (:,3,3) or (3,3)
|
P : numpy.array of shape (:,3,3) or (3,3)
|
||||||
1. Piola-Kirchhoff stress.
|
1. Piola-Kirchhoff stress.
|
||||||
|
|
||||||
"""
|
"""
|
||||||
if np.shape(F) == np.shape(P) == (3,3):
|
if np.shape(F) == np.shape(P) == (3,3):
|
||||||
sigma = 1.0/np.linalg.det(F) * np.dot(P,F.T)
|
sigma = 1.0/np.linalg.det(F) * np.dot(P,F.T)
|
||||||
else:
|
else:
|
||||||
sigma = np.einsum('i,ijk,ilk->ijl',1.0/np.linalg.det(F),P,F)
|
sigma = np.einsum('i,ijk,ilk->ijl',1.0/np.linalg.det(F),P,F)
|
||||||
return symmetric(sigma)
|
return symmetric(sigma)
|
||||||
|
|
||||||
|
|
||||||
def strain_tensor(F,t,m):
|
def strain_tensor(F,t,m):
|
||||||
"""
|
"""
|
||||||
Return strain tensor calculated from deformation gradient.
|
Return strain tensor calculated from deformation gradient.
|
||||||
|
|
||||||
For details refer to https://en.wikipedia.org/wiki/Finite_strain_theory and
|
For details refer to https://en.wikipedia.org/wiki/Finite_strain_theory and
|
||||||
https://de.wikipedia.org/wiki/Verzerrungstensor
|
https://de.wikipedia.org/wiki/Verzerrungstensor
|
||||||
|
|
||||||
Parameters
|
Parameters
|
||||||
----------
|
----------
|
||||||
F : numpy.array of shape (:,3,3) or (3,3)
|
F : numpy.array of shape (:,3,3) or (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.
|
||||||
m : float
|
m : float
|
||||||
Order of the strain.
|
Order of the strain.
|
||||||
|
|
||||||
"""
|
"""
|
||||||
F_ = F.reshape((1,3,3)) if F.shape == (3,3) else F
|
F_ = F.reshape((1,3,3)) if F.shape == (3,3) else F
|
||||||
if t == 'U':
|
if t == 'U':
|
||||||
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 == 'V':
|
elif t == 'V':
|
||||||
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:
|
||||||
eps = 1.0/(2.0*abs(m)) * (+ np.matmul(n,np.einsum('ij,ikj->ijk',w**m,n))
|
eps = 1.0/(2.0*abs(m)) * (+ np.matmul(n,np.einsum('ij,ikj->ijk',w**m,n))
|
||||||
- np.broadcast_to(np.ones(3),[F_.shape[0],3]))
|
- np.broadcast_to(np.ones(3),[F_.shape[0],3]))
|
||||||
elif m < 0.0:
|
elif m < 0.0:
|
||||||
eps = 1.0/(2.0*abs(m)) * (- np.matmul(n,np.einsum('ij,ikj->ijk',w**m,n))
|
eps = 1.0/(2.0*abs(m)) * (- np.matmul(n,np.einsum('ij,ikj->ijk',w**m,n))
|
||||||
+ np.broadcast_to(np.ones(3),[F_.shape[0],3]))
|
+ np.broadcast_to(np.ones(3),[F_.shape[0],3]))
|
||||||
else:
|
else:
|
||||||
eps = np.matmul(n,np.einsum('ij,ikj->ijk',0.5*np.log(w),n))
|
eps = np.matmul(n,np.einsum('ij,ikj->ijk',0.5*np.log(w),n))
|
||||||
|
|
||||||
return eps.reshape((3,3)) if np.shape(F) == (3,3) else \
|
return eps.reshape((3,3)) if np.shape(F) == (3,3) else \
|
||||||
eps
|
eps
|
||||||
|
|
||||||
|
|
||||||
def deviatoric_part(x):
|
def deviatoric_part(x):
|
||||||
"""
|
"""
|
||||||
Return deviatoric part of a tensor.
|
Return deviatoric part of a tensor.
|
||||||
|
|
||||||
Parameters
|
Parameters
|
||||||
----------
|
----------
|
||||||
x : numpy.array of shape (:,3,3) or (3,3)
|
x : numpy.array of shape (:,3,3) or (3,3)
|
||||||
Tensor of which the deviatoric part is computed.
|
Tensor of which the deviatoric part is computed.
|
||||||
|
|
||||||
"""
|
"""
|
||||||
return x - np.eye(3)*spherical_part(x) if np.shape(x) == (3,3) else \
|
return x - np.eye(3)*spherical_part(x) if np.shape(x) == (3,3) else \
|
||||||
x - np.einsum('ijk,i->ijk',np.broadcast_to(np.eye(3),[x.shape[0],3,3]),spherical_part(x))
|
x - np.einsum('ijk,i->ijk',np.broadcast_to(np.eye(3),[x.shape[0],3,3]),spherical_part(x))
|
||||||
|
|
||||||
|
|
||||||
def spherical_part(x):
|
def spherical_part(x):
|
||||||
"""
|
"""
|
||||||
Return spherical (hydrostatic) part of a tensor.
|
Return spherical (hydrostatic) part of a tensor.
|
||||||
|
|
||||||
A single scalar is returned, i.e. the hydrostatic part is not mapped on the 3rd order identity
|
A single scalar is returned, i.e. the hydrostatic part is not mapped on the 3rd order identity
|
||||||
matrix.
|
matrix.
|
||||||
|
|
||||||
Parameters
|
Parameters
|
||||||
----------
|
----------
|
||||||
x : numpy.array of shape (:,3,3) or (3,3)
|
x : numpy.array of shape (:,3,3) or (3,3)
|
||||||
Tensor of which the hydrostatic part is computed.
|
Tensor of which the hydrostatic part is computed.
|
||||||
|
|
||||||
"""
|
"""
|
||||||
return np.trace(x)/3.0 if np.shape(x) == (3,3) else \
|
return np.trace(x)/3.0 if np.shape(x) == (3,3) else \
|
||||||
np.trace(x,axis1=1,axis2=2)/3.0
|
np.trace(x,axis1=1,axis2=2)/3.0
|
||||||
|
|
||||||
|
|
||||||
def Mises_stress(sigma):
|
def Mises_stress(sigma):
|
||||||
"""
|
"""
|
||||||
Return the Mises equivalent of a stress tensor.
|
Return the Mises equivalent of a stress tensor.
|
||||||
|
|
||||||
Parameters
|
Parameters
|
||||||
----------
|
----------
|
||||||
sigma : numpy.array of shape (:,3,3) or (3,3)
|
sigma : numpy.array of shape (:,3,3) or (3,3)
|
||||||
Symmetric stress tensor of which the von Mises equivalent is computed.
|
Symmetric stress tensor of which the von Mises equivalent is computed.
|
||||||
|
|
||||||
"""
|
"""
|
||||||
s = deviatoric_part(sigma)
|
s = deviatoric_part(sigma)
|
||||||
return np.sqrt(3.0/2.0*(np.sum(s**2.0))) if np.shape(sigma) == (3,3) else \
|
return np.sqrt(3.0/2.0*(np.sum(s**2.0))) if np.shape(sigma) == (3,3) else \
|
||||||
np.sqrt(3.0/2.0*np.einsum('ijk->i',s**2.0))
|
np.sqrt(3.0/2.0*np.einsum('ijk->i',s**2.0))
|
||||||
|
|
||||||
|
|
||||||
def Mises_strain(epsilon):
|
def Mises_strain(epsilon):
|
||||||
"""
|
"""
|
||||||
Return the Mises equivalent of a strain tensor.
|
Return the Mises equivalent of a strain tensor.
|
||||||
|
|
||||||
Parameters
|
Parameters
|
||||||
----------
|
----------
|
||||||
epsilon : numpy.array of shape (:,3,3) or (3,3)
|
epsilon : numpy.array of shape (:,3,3) or (3,3)
|
||||||
Symmetric strain tensor of which the von Mises equivalent is computed.
|
Symmetric strain tensor of which the von Mises equivalent is computed.
|
||||||
|
|
||||||
"""
|
"""
|
||||||
s = deviatoric_part(epsilon)
|
s = deviatoric_part(epsilon)
|
||||||
return np.sqrt(2.0/3.0*(np.sum(s**2.0))) if np.shape(epsilon) == (3,3) else \
|
return np.sqrt(2.0/3.0*(np.sum(s**2.0))) if np.shape(epsilon) == (3,3) else \
|
||||||
np.sqrt(2.0/3.0*np.einsum('ijk->i',s**2.0))
|
np.sqrt(2.0/3.0*np.einsum('ijk->i',s**2.0))
|
||||||
|
|
||||||
|
|
||||||
def symmetric(x):
|
def symmetric(x):
|
||||||
"""
|
"""
|
||||||
Return the symmetrized tensor.
|
Return the symmetrized tensor.
|
||||||
|
|
||||||
Parameters
|
Parameters
|
||||||
----------
|
----------
|
||||||
x : numpy.array of shape (:,3,3) or (3,3)
|
x : numpy.array of shape (:,3,3) or (3,3)
|
||||||
Tensor of which the symmetrized values are computed.
|
Tensor of which the symmetrized values are computed.
|
||||||
|
|
||||||
"""
|
"""
|
||||||
return (x+transpose(x))*0.5
|
return (x+transpose(x))*0.5
|
||||||
|
|
||||||
|
|
||||||
def maximum_shear(x):
|
def maximum_shear(x):
|
||||||
"""
|
"""
|
||||||
Return the maximum shear component of a symmetric tensor.
|
Return the maximum shear component of a symmetric tensor.
|
||||||
|
|
||||||
Parameters
|
Parameters
|
||||||
----------
|
----------
|
||||||
x : numpy.array of shape (:,3,3) or (3,3)
|
x : numpy.array of shape (:,3,3) or (3,3)
|
||||||
Symmetric tensor of which the maximum shear is computed.
|
Symmetric tensor of which the maximum shear is computed.
|
||||||
|
|
||||||
"""
|
"""
|
||||||
w = np.linalg.eigvalsh(symmetric(x)) # eigenvalues in ascending order
|
w = np.linalg.eigvalsh(symmetric(x)) # eigenvalues in ascending order
|
||||||
return (w[2] - w[0])*0.5 if np.shape(x) == (3,3) else \
|
return (w[2] - w[0])*0.5 if np.shape(x) == (3,3) else \
|
||||||
(w[:,2] - w[:,0])*0.5
|
(w[:,2] - w[:,0])*0.5
|
||||||
|
|
||||||
|
|
||||||
def principal_components(x):
|
def principal_components(x):
|
||||||
"""
|
"""
|
||||||
Return the principal components of a symmetric tensor.
|
Return the principal components of a symmetric tensor.
|
||||||
|
|
||||||
The principal components (eigenvalues) are sorted in descending order, each repeated according to
|
The principal components (eigenvalues) are sorted in descending order, each repeated according to
|
||||||
its multiplicity.
|
its multiplicity.
|
||||||
|
|
||||||
Parameters
|
Parameters
|
||||||
----------
|
----------
|
||||||
x : numpy.array of shape (:,3,3) or (3,3)
|
x : numpy.array of shape (:,3,3) or (3,3)
|
||||||
Symmetric tensor of which the principal compontents are computed.
|
Symmetric tensor of which the principal compontents are computed.
|
||||||
|
|
||||||
"""
|
"""
|
||||||
w = np.linalg.eigvalsh(symmetric(x)) # eigenvalues in ascending order
|
w = np.linalg.eigvalsh(symmetric(x)) # eigenvalues in ascending order
|
||||||
return w[::-1] if np.shape(x) == (3,3) else \
|
return w[::-1] if np.shape(x) == (3,3) else \
|
||||||
w[:,::-1]
|
w[:,::-1]
|
||||||
|
|
||||||
|
|
||||||
def transpose(x):
|
def transpose(x):
|
||||||
"""
|
"""
|
||||||
Return the transpose of a tensor.
|
Return the transpose of a tensor.
|
||||||
|
|
||||||
Parameters
|
Parameters
|
||||||
----------
|
----------
|
||||||
x : numpy.array of shape (:,3,3) or (3,3)
|
x : numpy.array of shape (:,3,3) or (3,3)
|
||||||
Tensor of which the transpose is computed.
|
Tensor of which the transpose is computed.
|
||||||
|
|
||||||
"""
|
"""
|
||||||
return x.T if np.shape(x) == (3,3) else \
|
return x.T if np.shape(x) == (3,3) else \
|
||||||
np.transpose(x,(0,2,1))
|
np.transpose(x,(0,2,1))
|
||||||
|
|
||||||
|
|
||||||
def rotational_part(x):
|
def rotational_part(x):
|
||||||
"""
|
"""
|
||||||
Return the rotational part of a tensor.
|
Return the rotational part of a tensor.
|
||||||
|
|
||||||
Parameters
|
Parameters
|
||||||
----------
|
----------
|
||||||
x : numpy.array of shape (:,3,3) or (3,3)
|
x : numpy.array of shape (:,3,3) or (3,3)
|
||||||
Tensor of which the rotational part is computed.
|
Tensor of which the rotational part is computed.
|
||||||
|
|
||||||
"""
|
"""
|
||||||
return __polar_decomposition(x,'R')
|
return __polar_decomposition(x,'R')
|
||||||
|
|
||||||
|
|
||||||
def left_stretch(x):
|
def left_stretch(x):
|
||||||
"""
|
"""
|
||||||
Return the left stretch of a tensor.
|
Return the left stretch of a tensor.
|
||||||
|
|
||||||
Parameters
|
Parameters
|
||||||
----------
|
----------
|
||||||
x : numpy.array of shape (:,3,3) or (3,3)
|
x : numpy.array of shape (:,3,3) or (3,3)
|
||||||
Tensor of which the left stretch is computed.
|
Tensor of which the left stretch is computed.
|
||||||
|
|
||||||
"""
|
"""
|
||||||
return __polar_decomposition(x,'V')
|
return __polar_decomposition(x,'V')
|
||||||
|
|
||||||
|
|
||||||
def right_stretch(x):
|
def right_stretch(x):
|
||||||
"""
|
"""
|
||||||
Return the right stretch of a tensor.
|
Return the right stretch of a tensor.
|
||||||
|
|
||||||
Parameters
|
Parameters
|
||||||
----------
|
----------
|
||||||
x : numpy.array of shape (:,3,3) or (3,3)
|
x : numpy.array of shape (:,3,3) or (3,3)
|
||||||
Tensor of which the right stretch is computed.
|
Tensor of which the right stretch is computed.
|
||||||
|
|
||||||
"""
|
"""
|
||||||
return __polar_decomposition(x,'U')
|
return __polar_decomposition(x,'U')
|
||||||
|
|
||||||
|
|
||||||
def __polar_decomposition(x,requested):
|
def __polar_decomposition(x,requested):
|
||||||
"""
|
"""
|
||||||
Singular value decomposition.
|
Singular value decomposition.
|
||||||
|
|
||||||
Parameters
|
Parameters
|
||||||
----------
|
----------
|
||||||
x : numpy.array of shape (:,3,3) or (3,3)
|
x : numpy.array of shape (:,3,3) or (3,3)
|
||||||
Tensor of which the singular values are computed.
|
Tensor of which the singular values are computed.
|
||||||
requested : list of str
|
requested : list of str
|
||||||
Requested outputs: ‘R’ for the rotation tensor,
|
Requested outputs: ‘R’ for the rotation tensor,
|
||||||
‘V’ for left stretch tensor and ‘U’ for right stretch tensor.
|
‘V’ for left stretch tensor and ‘U’ for right stretch tensor.
|
||||||
|
|
||||||
"""
|
"""
|
||||||
u, s, vh = np.linalg.svd(x)
|
u, s, vh = np.linalg.svd(x)
|
||||||
R = np.dot(u,vh) if np.shape(x) == (3,3) else \
|
R = np.dot(u,vh) if np.shape(x) == (3,3) else \
|
||||||
np.einsum('ijk,ikl->ijl',u,vh)
|
np.einsum('ijk,ikl->ijl',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(x,R.T) if np.shape(x) == (3,3) else np.einsum('ijk,ilk->ijl',x,R))
|
output.append(np.dot(x,R.T) if np.shape(x) == (3,3) else np.einsum('ijk,ilk->ijl',x,R))
|
||||||
if 'U' in requested:
|
if 'U' in requested:
|
||||||
output.append(np.dot(R.T,x) if np.shape(x) == (3,3) else np.einsum('ikj,ikl->ijl',R,x))
|
output.append(np.dot(R.T,x) if np.shape(x) == (3,3) else np.einsum('ikj,ikl->ijl',R,x))
|
||||||
|
|
||||||
return tuple(output)
|
return tuple(output)
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue