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

einsum simplifications

This commit is contained in:
Martin Diehl 2020-11-16 07:12:37 +01:00
parent 9b9d83d93c
commit 5ebde607a2
4 changed files with 23 additions and 26 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
python/damask/_orientation.py
View File

@ -1103,7 +1103,7 @@ class Orientation(Rotation):
(np.array(hkil),np.array([[1,0,0,0], (np.array(hkil),np.array([[1,0,0,0],
[0,1,0,0], [0,1,0,0],
[0,0,0,1]])) [0,0,0,1]]))
return np.einsum('il,...l->...i',basis,axis) return np.einsum('il,...l',basis,axis)
@classmethod @classmethod
@ -1133,7 +1133,7 @@ class Orientation(Rotation):
[ 0, 1, 0], [ 0, 1, 0],
[-1,-1, 0], [-1,-1, 0],
[ 0, 0, 1]])) [ 0, 0, 1]]))
return np.einsum('il,...l->...i',basis,axis) return np.einsum('il,...l',basis,axis)
def to_lattice(self,direction=None,plane=None): def to_lattice(self,direction=None,plane=None):
@ -1157,7 +1157,7 @@ class Orientation(Rotation):
axis,basis = (np.array(direction),self.basis_reciprocal.T) \ axis,basis = (np.array(direction),self.basis_reciprocal.T) \
if plane is None else \ if plane is None else \
(np.array(plane),self.basis_real.T) (np.array(plane),self.basis_real.T)
return np.einsum('il,...l->...i',basis,axis) return np.einsum('il,...l',basis,axis)
def to_frame(self,uvw=None,hkl=None,with_symmetry=False): def to_frame(self,uvw=None,hkl=None,with_symmetry=False):
@ -1183,9 +1183,9 @@ class Orientation(Rotation):
if hkl is None else \ if hkl is None else \
(np.array(hkl),self.basis_reciprocal) (np.array(hkl),self.basis_reciprocal)
return (self.symmetry_operations.broadcast_to(self.symmetry_operations.shape+axis.shape[:-1],mode='right') return (self.symmetry_operations.broadcast_to(self.symmetry_operations.shape+axis.shape[:-1],mode='right')
@ np.broadcast_to(np.einsum('il,...l->...i',basis,axis),self.symmetry_operations.shape+axis.shape) @ np.broadcast_to(np.einsum('il,...l',basis,axis),self.symmetry_operations.shape+axis.shape)
if with_symmetry else if with_symmetry else
np.einsum('il,...l->...i',basis,axis)) np.einsum('il,...l',basis,axis))
def to_pole(self,uvw=None,hkl=None,with_symmetry=False): def to_pole(self,uvw=None,hkl=None,with_symmetry=False):
@ -1227,7 +1227,7 @@ class Orientation(Rotation):
""" """
d = self.to_frame(uvw=self.kinematics[mode]['direction'],with_symmetry=False) d = self.to_frame(uvw=self.kinematics[mode]['direction'],with_symmetry=False)
p = self.to_frame(hkl=self.kinematics[mode]['plane'] ,with_symmetry=False) p = self.to_frame(hkl=self.kinematics[mode]['plane'] ,with_symmetry=False)
P = np.einsum('...i,...j->...ij',d/np.linalg.norm(d,axis=-1,keepdims=True), P = np.einsum('...i,...j',d/np.linalg.norm(d,axis=-1,keepdims=True),
p/np.linalg.norm(p,axis=-1,keepdims=True)) p/np.linalg.norm(p,axis=-1,keepdims=True))
return ~self.broadcast_to( self.shape+P.shape[:-2],mode='right') \ return ~self.broadcast_to( self.shape+P.shape[:-2],mode='right') \

2
python/damask/_rotation.py
View File

@ -553,7 +553,7 @@ class Rotation:
orthonormal = False # contains stretch orthonormal = False # contains stretch
if not orthonormal: if not orthonormal:
(U,S,Vh) = np.linalg.svd(om) # singular value decomposition (U,S,Vh) = np.linalg.svd(om) # singular value decomposition
om = np.einsum('...ij,...jl->...il',U,Vh) om = np.einsum('...ij,...jl',U,Vh)
if not np.all(np.isclose(np.linalg.det(om),1.0)): if not np.all(np.isclose(np.linalg.det(om),1.0)):
raise ValueError('Orientation matrix has determinant ≠ 1.') raise ValueError('Orientation matrix has determinant ≠ 1.')
if not np.all(np.isclose(np.einsum('...i,...i',om[...,0],om[...,1]), 0.0)) \ if not np.all(np.isclose(np.einsum('...i,...i',om[...,0],om[...,1]), 0.0)) \

26
python/damask/mechanics.py
View File

@ -60,7 +60,7 @@ def Cauchy(P,F):
Cauchy stress. Cauchy stress.
""" """
sigma = _np.einsum('...,...ij,...kj->...ik',1.0/_np.linalg.det(F),P,F) sigma = _np.einsum('...,...ij,...kj',1.0/_np.linalg.det(F),P,F)
return tensor.symmetric(sigma) return tensor.symmetric(sigma)
@ -79,7 +79,7 @@ def deviatoric_part(T):
Deviatoric part of T. Deviatoric part of T.
""" """
return T - _np.einsum('...ij,...->...ij',_np.eye(3),spherical_part(T)) return T - _np.einsum('...ij,...',_np.eye(3),spherical_part(T))
def maximum_shear(T_sym): def maximum_shear(T_sym):
@ -157,7 +157,7 @@ def PK2(P,F):
Second Piola-Kirchhoff stress. Second Piola-Kirchhoff stress.
""" """
S = _np.einsum('...jk,...kl->...jl',_np.linalg.inv(F),P) S = _np.einsum('...jk,...kl',_np.linalg.inv(F),P)
return tensor.symmetric(S) return tensor.symmetric(S)
@ -199,7 +199,7 @@ def spherical_part(T,tensor=False):
""" """
sph = _np.trace(T,axis2=-2,axis1=-1)/3.0 sph = _np.trace(T,axis2=-2,axis1=-1)/3.0
return _np.einsum('...jk,...->...jk',_np.eye(3),sph) if tensor else sph return _np.einsum('...jk,...',_np.eye(3),sph) if tensor else sph
def strain(F,t,m): def strain(F,t,m):
@ -231,14 +231,12 @@ def strain(F,t,m):
w,n = _np.linalg.eigh(Cauchy_Green_deformation_right(F)) w,n = _np.linalg.eigh(Cauchy_Green_deformation_right(F))
if m > 0.0: if m > 0.0:
eps = 1.0/(2.0*abs(m)) * (+ _np.matmul(n,_np.einsum('...j,...kj->...jk',w**m,n)) eps = 1.0/(2.0*abs(m)) * (+ _np.einsum('...j,...kj,...lj',w**m,n,n) - _np.eye(3))
- _np.einsum('...jk->...jk',_np.eye(3)))
elif m < 0.0: elif m < 0.0:
eps = 1.0/(2.0*abs(m)) * (- _np.matmul(n,_np.einsum('...j,...kj->...jk',w**m,n)) eps = 1.0/(2.0*abs(m)) * (- _np.einsum('...j,...kj,...lj',w**m,n,n) + _np.eye(3))
+ _np.einsum('...jk->...jk',_np.eye(3)))
else: else:
eps = _np.matmul(n,_np.einsum('...j,...kj->...jk',0.5*_np.log(w),n)) eps = _np.einsum('...j,...kj,...lj',0.5*_np.log(w),n,n)
return eps return eps
@ -293,22 +291,22 @@ def _polar_decomposition(T,requested):
""" """
u, _, vh = _np.linalg.svd(T) u, _, vh = _np.linalg.svd(T)
R = _np.einsum('...ij,...jk->...ik',u,vh) R = _np.einsum('...ij,...jk',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.einsum('...ij,...kj->...ik',T,R)) output.append(_np.einsum('...ij,...kj',T,R))
if 'U' in requested: if 'U' in requested:
output.append(_np.einsum('...ji,...jk->...ik',R,T)) output.append(_np.einsum('...ji,...jk',R,T))
return tuple(output) return tuple(output)
def _Mises(T_sym,s): def _Mises(T_sym,s):
""" """
Base equation for Mises equivalent of a stres or strain tensor. Base equation for Mises equivalent of a stress or strain tensor.
Parameters Parameters
---------- ----------
@ -319,7 +317,7 @@ def _Mises(T_sym,s):
""" """
d = deviatoric_part(T_sym) d = deviatoric_part(T_sym)
return _np.sqrt(s*_np.einsum('...jk->...',d**2.0)) return _np.sqrt(s*_np.sum(d**2.0,axis=(-1,-2)))
# for compatibility # for compatibility

5
python/tests/test_mechanics.py
View File

@ -127,8 +127,7 @@ class TestMechanics:
for i,v in enumerate(np.reshape(vectorized(epsilon_vec),vectorized(epsilon).shape)): for i,v in enumerate(np.reshape(vectorized(epsilon_vec),vectorized(epsilon).shape)):
assert np.allclose(single(epsilon[i]),v) assert np.allclose(single(epsilon[i]),v)
@pytest.mark.parametrize('vectorized,single',[ @pytest.mark.parametrize('vectorized,single',[(mechanics.Cauchy,Cauchy),
(mechanics.Cauchy,Cauchy),
(mechanics.PK2 ,PK2 ) (mechanics.PK2 ,PK2 )
]) ])
def test_vectorize_2_arg(self,vectorized,single): def test_vectorize_2_arg(self,vectorized,single):
@ -227,7 +226,7 @@ class TestMechanics:
scalar**3.0) scalar**3.0)
def test_spherical_Mises(self): def test_spherical_Mises(self):
"""Ensure that Mises equivalent strrain of spherical strain is 0.""" """Ensure that Mises equivalent strain of spherical strain is 0."""
x = np.random.rand(self.n,3,3) x = np.random.rand(self.n,3,3)
sph = mechanics.spherical_part(x,True) sph = mechanics.spherical_part(x,True)
assert np.allclose(mechanics.Mises_strain(sph), assert np.allclose(mechanics.Mises_strain(sph),