From 987c4a9e8d01856de7689cf46e8cb74a9f922e9a Mon Sep 17 00:00:00 2001 From: Samad Vakili Date: Sun, 24 May 2020 18:30:44 +0200 Subject: [PATCH 01/35] first commit mechanics --- python/damask/mechanics.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/python/damask/mechanics.py b/python/damask/mechanics.py index c81399d94..a144e40f5 100644 --- a/python/damask/mechanics.py +++ b/python/damask/mechanics.py @@ -17,7 +17,8 @@ def Cauchy(P,F): if _np.shape(F) == _np.shape(P) == (3,3): 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('i,ijk,ilk->ijl',1.0/_np.linalg.det(F),P,F) + sigma = _np.einsum('...,...ij,...kj->...ij',1.0/_np.linalg.det(F),P,F) return symmetric(sigma) @@ -190,7 +191,7 @@ def spherical_part(T,tensor=False): sph = _np.trace(T)/3.0 return sph if not tensor else _np.eye(3)*sph else: - sph = _np.trace(T,axis1=1,axis2=2)/3.0 + sph = _np.trace(T,axis2=-2,axis1=-1)/3.0 if not tensor: return sph else: From f9c33d9210ebd2638d996653e318a6a115ea80f1 Mon Sep 17 00:00:00 2001 From: Samad Vakili Date: Tue, 26 May 2020 16:27:27 +0200 Subject: [PATCH 02/35] mechanics checked for an array with arbitrary dimensions --- python/damask/mechanics.py | 31 +++++++++++++++++-------------- 1 file changed, 17 insertions(+), 14 deletions(-) diff --git a/python/damask/mechanics.py b/python/damask/mechanics.py index a144e40f5..e388b3569 100644 --- a/python/damask/mechanics.py +++ b/python/damask/mechanics.py @@ -18,7 +18,7 @@ def Cauchy(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->...ij',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) @@ -33,8 +33,7 @@ def deviatoric_part(T): """ return T - _np.eye(3)*spherical_part(T) if _np.shape(T) == (3,3) else \ - T - _np.einsum('ijk,i->ijk',_np.broadcast_to(_np.eye(3),[T.shape[0],3,3]),spherical_part(T)) - + T - _np.einsum('...ij,...->...ij',_np.eye(3),spherical_part(T)) def eigenvalues(T_sym): """ @@ -101,7 +100,7 @@ def maximum_shear(T_sym): """ w = eigenvalues(T_sym) return (w[0] - w[2])*0.5 if _np.shape(T_sym) == (3,3) else \ - (w[:,0] - w[:,2])*0.5 + (w[...,0] - w[...,2])*0.5 def Mises_strain(epsilon): @@ -195,7 +194,8 @@ def spherical_part(T,tensor=False): if not tensor: return sph else: - return _np.einsum('ijk,i->ijk',_np.broadcast_to(_np.eye(3),(T.shape[0],3,3)),sph) + #return _np.einsum('ijk,i->ijk',_np.broadcast_to(_np.eye(3),(T.shape[0],3,3)),sph) + return _np.einsum('...jk,...->...jk',_np.eye(3),sph) def strain_tensor(F,t,m): @@ -224,13 +224,16 @@ def strain_tensor(F,t,m): w,n = _np.linalg.eigh(C) if m > 0.0: - eps = 1.0/(2.0*abs(m)) * (+ _np.matmul(n,_np.einsum('ij,ikj->ijk',w**m,n)) - - _np.broadcast_to(_np.eye(3),[F_.shape[0],3,3])) + #eps = 1.0/(2.0*abs(m)) * (+ _np.matmul(n,_np.einsum('ij,ikj->ijk',w**m,n)) + # - _np.broadcast_to(_np.eye(3),[F_.shape[0],3,3])) + eps = 1.0/(2.0*abs(m)) * (+ _np.matmul(n,_np.einsum('...j,...kj->...jk',w**m,n)) + - _np.einsum('...jk->...jk',_np.eye(3))) + elif m < 0.0: - eps = 1.0/(2.0*abs(m)) * (- _np.matmul(n,_np.einsum('ij,ikj->ijk',w**m,n)) - + _np.broadcast_to(_np.eye(3),[F_.shape[0],3,3])) + eps = 1.0/(2.0*abs(m)) * (- _np.matmul(n,_np.einsum('...j,...kj->...jk',w**m,n)) + + _np.einsum('...jk->...jk',_np.eye(3))) else: - eps = _np.matmul(n,_np.einsum('ij,ikj->ijk',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 \ eps @@ -278,15 +281,15 @@ def _polar_decomposition(T,requested): """ u, s, vh = _np.linalg.svd(T) R = _np.dot(u,vh) if _np.shape(T) == (3,3) else \ - _np.einsum('ijk,ikl->ijl',u,vh) + _np.einsum('...ij,...jk->...ik',u,vh) output = [] if 'R' in requested: output.append(R) if 'V' in requested: - output.append(_np.dot(T,R.T) if _np.shape(T) == (3,3) else _np.einsum('ijk,ilk->ijl',T,R)) + output.append(_np.dot(T,R.T) if _np.shape(T) == (3,3) else _np.einsum('...ij,...kj->...ik',T,R)) if 'U' in requested: - output.append(_np.dot(R.T,T) if _np.shape(T) == (3,3) else _np.einsum('ikj,ikl->ijl',R,T)) + output.append(_np.dot(R.T,T) if _np.shape(T) == (3,3) else _np.einsum('...ji,...jk->...ik',R,T)) return tuple(output) @@ -305,4 +308,4 @@ def _Mises(T_sym,s): """ d = deviatoric_part(T_sym) return _np.sqrt(s*(_np.sum(d**2.0))) if _np.shape(T_sym) == (3,3) else \ - _np.sqrt(s*_np.einsum('ijk->i',d**2.0)) + _np.sqrt(s*_np.einsum('...jk->...',d**2.0)) From 56afc03f3cb108c63fc394e7ab2413864fac6a12 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 27 May 2020 18:05:08 +0200 Subject: [PATCH 03/35] only vectorized version needed use single point/simple versions only for testing --- python/damask/mechanics.py | 20 ++++++-------------- python/tests/test_mechanics.py | 15 +++++++++++++++ 2 files changed, 21 insertions(+), 14 deletions(-) diff --git a/python/damask/mechanics.py b/python/damask/mechanics.py index e388b3569..484341772 100644 --- a/python/damask/mechanics.py +++ b/python/damask/mechanics.py @@ -28,12 +28,12 @@ def deviatoric_part(T): Parameters ---------- - T : numpy.ndarray of shape (:,3,3) or (3,3) + T : numpy.ndarray of shape (...,3,3) Tensor of which the deviatoric part is computed. """ - return T - _np.eye(3)*spherical_part(T) if _np.shape(T) == (3,3) else \ - T - _np.einsum('...ij,...->...ij',_np.eye(3),spherical_part(T)) + return T - _np.einsum('...ij,...->...ij',_np.eye(3),spherical_part(T)) + def eigenvalues(T_sym): """ @@ -180,22 +180,14 @@ def spherical_part(T,tensor=False): Parameters ---------- - T : numpy.ndarray of shape (:,3,3) or (3,3) + T : numpy.ndarray of shape (...,3,3) Tensor of which the hydrostatic part is computed. tensor : bool, optional Map spherical part onto identity tensor. Default is false """ - if T.shape == (3,3): - sph = _np.trace(T)/3.0 - return sph if not tensor else _np.eye(3)*sph - else: - sph = _np.trace(T,axis2=-2,axis1=-1)/3.0 - if not tensor: - return sph - else: - #return _np.einsum('ijk,i->ijk',_np.broadcast_to(_np.eye(3),(T.shape[0],3,3)),sph) - return _np.einsum('...jk,...->...jk',_np.eye(3),sph) + sph = _np.trace(T,axis2=-2,axis1=-1)/3.0 + return _np.einsum('...jk,...->...jk',_np.eye(3),sph) if tensor else sph def strain_tensor(F,t,m): diff --git a/python/tests/test_mechanics.py b/python/tests/test_mechanics.py index f7128401f..da9e83f7f 100644 --- a/python/tests/test_mechanics.py +++ b/python/tests/test_mechanics.py @@ -3,12 +3,27 @@ import numpy as np from damask import mechanics +def deviatoric_part(T): + return T - np.eye(3)*spherical_part(T) + +def spherical_part(T,tensor=False): + sph = np.trace(T)/3.0 + return sph if not tensor else np.eye(3)*sph + class TestMechanics: n = 1000 c = np.random.randint(n) + @pytest.mark.parametrize('vectorized,single',[(mechanics.deviatoric_part, deviatoric_part), + (mechanics.spherical_part, spherical_part) + ]) + def test_vectorize_1_arg_(self,vectorized,single): + test_data = np.random.rand(self.n,3,3) + for i,v in enumerate(vectorized(test_data)): + assert np.allclose(single(test_data[i]),v) + @pytest.mark.parametrize('function',[mechanics.deviatoric_part, mechanics.eigenvalues, mechanics.eigenvectors, From 00168a5939b669deab27d1b1b30dcf7bb47cba40 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 27 May 2020 21:07:48 +0200 Subject: [PATCH 04/35] testing 2-dim array of tensors --- python/tests/test_mechanics.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/python/tests/test_mechanics.py b/python/tests/test_mechanics.py index da9e83f7f..5c757c71d 100644 --- a/python/tests/test_mechanics.py +++ b/python/tests/test_mechanics.py @@ -20,9 +20,10 @@ class TestMechanics: (mechanics.spherical_part, spherical_part) ]) def test_vectorize_1_arg_(self,vectorized,single): - test_data = np.random.rand(self.n,3,3) - for i,v in enumerate(vectorized(test_data)): - assert np.allclose(single(test_data[i]),v) + test_data_flat = np.random.rand(self.n,3,3) + test_data = np.reshape(test_data_flat,(self.n//10,10,3,3)) + for i,v in enumerate(np.reshape(vectorized(test_data),vectorized(test_data_flat).shape)): + assert np.allclose(single(test_data_flat[i]),v) @pytest.mark.parametrize('function',[mechanics.deviatoric_part, mechanics.eigenvalues, From 865a505186b6937480a45f55b2f63460612e4d9c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 27 May 2020 22:32:35 +0200 Subject: [PATCH 05/35] fix for logarithmic strain --- python/damask/mechanics.py | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/python/damask/mechanics.py b/python/damask/mechanics.py index 484341772..c7ba39bcd 100644 --- a/python/damask/mechanics.py +++ b/python/damask/mechanics.py @@ -216,8 +216,6 @@ def strain_tensor(F,t,m): w,n = _np.linalg.eigh(C) if m > 0.0: - #eps = 1.0/(2.0*abs(m)) * (+ _np.matmul(n,_np.einsum('ij,ikj->ijk',w**m,n)) - # - _np.broadcast_to(_np.eye(3),[F_.shape[0],3,3])) eps = 1.0/(2.0*abs(m)) * (+ _np.matmul(n,_np.einsum('...j,...kj->...jk',w**m,n)) - _np.einsum('...jk->...jk',_np.eye(3))) @@ -225,7 +223,7 @@ def strain_tensor(F,t,m): eps = 1.0/(2.0*abs(m)) * (- _np.matmul(n,_np.einsum('...j,...kj->...jk',w**m,n)) + _np.einsum('...jk->...jk',_np.eye(3))) 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 \ eps From 627e62439b1eb94a069b4b7ccc011f2c404778e3 Mon Sep 17 00:00:00 2001 From: Samad Vakili Date: Tue, 9 Jun 2020 14:17:31 +0200 Subject: [PATCH 06/35] the single case is included in test file --- python/tests/test_mechanics.py | 164 +++++++++++++++++++++++++++------ 1 file changed, 136 insertions(+), 28 deletions(-) mode change 100644 => 100755 python/tests/test_mechanics.py diff --git a/python/tests/test_mechanics.py b/python/tests/test_mechanics.py old mode 100644 new mode 100755 index 5c757c71d..0db6a1ba9 --- a/python/tests/test_mechanics.py +++ b/python/tests/test_mechanics.py @@ -3,13 +3,110 @@ import numpy as np from damask import mechanics +def Cauchy(P,F): + sigma = 1.0/np.linalg.det(F) * np.dot(P,F.T) + return mechanics.symmetric(sigma) + + def deviatoric_part(T): return T - np.eye(3)*spherical_part(T) + +def eigenvalues(T_sym): + return np.linalg.eigvalsh(symmetric(T_sym)) + + +def eigenvectors(T_sym,RHS=False): + (u,v) = np.linalg.eigh(symmetric(T_sym)) + + if RHS: + if np.linalg.det(v) < 0.0: v[:,2] *= -1.0 + return v + + +def left_stretch(T): + return polar_decomposition(T,'V')[0] + + +def maximum_shear(T_sym): + w = eigenvalues(T_sym) + return (w[0] - w[2])*0.5 + + +def Mises_strain(epsilon): + return Mises(epsilon,2.0/3.0) + + +def Mises_stress(sigma): + return Mises(sigma,3.0/2.0) + + +def PK2(P,F): + S = np.dot(_np.linalg.inv(F),P) + return symmetric(S) + + +def right_stretch(T): + return polar_decomposition(T,'U')[0] + + +def rotational_part(T): + return polar_decomposition(T,'R')[0] + def spherical_part(T,tensor=False): sph = np.trace(T)/3.0 return sph if not tensor else np.eye(3)*sph + +def strain_tensor(F,t,m): + F_ = F.reshape(1,3,3) + + if t == 'V': + B = np.matmul(F_,transpose(F_)) + w,n = np.linalg.eigh(B) + elif t == 'U': + C = np.matmul(transpose(F_),F_) + w,n = np.linalg.eigh(C) + + if m > 0.0: + eps = 1.0/(2.0*abs(m)) * (+ np.matmul(n,_np.einsum('ij,ikj->ijk',w**m,n)) + - np.einsum('ijk->ijk',np.eye(3))) + elif m < 0.0: + eps = 1.0/(2.0*abs(m)) * (- np.matmul(n,_np.einsum('ij,ikj->ijk',w**m,n)) + + np.einsum('ijk->ijk',np.eye(3))) + else: + eps = np.matmul(n,np.einsum('ij,ikj->ijk',0.5*np.log(w),n)) + + return eps.reshape(3,3) + + +def symmetric(T): + return (T+transpose(T))*0.5 + + +def transpose(T): + return T.T + + +def polar_decomposition(T,requested): + u, s, vh = np.linalg.svd(T) + R = np.dot(u,vh) + + output = [] + if 'R' in requested: + output.append(R) + if 'V' in requested: + output.append(np.dot(T,R.T)) + if 'U' in requested: + output.append(np.dot(R.T,T)) + + return tuple(output) + +def Mises(T_sym,s): + d = deviatoric_part(T_sym) + return np.sqrt(s*(np.sum(d**2.0))) + + class TestMechanics: n = 1000 @@ -20,42 +117,53 @@ class TestMechanics: (mechanics.spherical_part, spherical_part) ]) def test_vectorize_1_arg_(self,vectorized,single): + print("done") test_data_flat = np.random.rand(self.n,3,3) test_data = np.reshape(test_data_flat,(self.n//10,10,3,3)) for i,v in enumerate(np.reshape(vectorized(test_data),vectorized(test_data_flat).shape)): assert np.allclose(single(test_data_flat[i]),v) - @pytest.mark.parametrize('function',[mechanics.deviatoric_part, - mechanics.eigenvalues, - mechanics.eigenvectors, - mechanics.left_stretch, - mechanics.maximum_shear, - mechanics.Mises_strain, - mechanics.Mises_stress, - mechanics.right_stretch, - mechanics.rotational_part, - mechanics.spherical_part, - mechanics.symmetric, - mechanics.transpose, + @pytest.mark.parametrize('vectorized,single',[ + (mechanics.deviatoric_part, deviatoric_part), + (mechanics.eigenvalues , eigenvalues ), + (mechanics.eigenvectors , eigenvectors ), + (mechanics.left_stretch , left_stretch ), + (mechanics.maximum_shear , maximum_shear ), + (mechanics.Mises_strain , Mises_strain ), + (mechanics.Mises_stress , Mises_stress ), + (mechanics.right_stretch , right_stretch ), + (mechanics.rotational_part, rotational_part), + (mechanics.spherical_part , spherical_part ), + (mechanics.symmetric , symmetric ), + (mechanics.transpose , transpose ), ]) - def test_vectorize_1_arg(self,function): - epsilon = np.random.rand(self.n,3,3) - assert np.allclose(function(epsilon)[self.c],function(epsilon[self.c])) + def test_vectorize_1_arg(self,vectorized,single): + epsilon = np.random.rand(self.n,3,3) + epsilon_vec = np.reshape(epsilon,(self.n//10,10,3,3)) + for i,v in enumerate(np.reshape(vectorized(epsilon_vec),vectorized(epsilon).shape)): + assert np.allclose(single(epsilon[i]),v) - @pytest.mark.parametrize('function',[mechanics.Cauchy, - mechanics.PK2, - ]) - def test_vectorize_2_arg(self,function): - P = np.random.rand(self.n,3,3) - F = np.random.rand(self.n,3,3) - assert np.allclose(function(P,F)[self.c],function(P[self.c],F[self.c])) + @pytest.mark.parametrize('vectorized,single',[ + (mechanics.Cauchy,Cauchy), + (mechanics.PK2 ,PK2 ) + ]) + def test_vectorize_2_arg(self,vectorized,single): + P = np.random.rand(self.n,3,3) + F = np.random.rand(self.n,3,3) + P_vec = np.random.rand(self.n//10,10,3,3) + F_vec = np.random.rand(self.n//10,10,3,3) + for i,v in enumerate(np.reshape(vectorized(P_vec,F_vec),vectorized(P,F).shape)): + assert np.allclose(single(P[i],F[i]),v) - def test_vectorize_strain_tensor(self): - F = np.random.rand(self.n,3,3) - t = ['V','U'][np.random.randint(0,2)] - m = np.random.random()*10. -5.0 - assert np.allclose(mechanics.strain_tensor(F,t,m)[self.c], - mechanics.strain_tensor(F[self.c],t,m)) + + @pytest.mark.parametrize('vectorized,single',[(mechanics.strain_tensor,strain_tensor)]) + def test_vectorize_strain_tensor(self,vectroized,single): + F = np.random.rand(self.n,3,3) + F_vec = np.random.rand(self.n//10,10,3,3) + t = ['V','U'][np.random.randint(0,2)] + m = np.random.random()*10.0 -5.0 + for i,v in enumerate(np.reshape(vectorized(F_vec,t,m),vectorized(F,t,m).shape)): + assert np.allcloase(single(F[i],t,m),v) @pytest.mark.parametrize('function',[mechanics.Cauchy, mechanics.PK2, From 824682d236504485bd4bd06b207cd6cd4e4a4e84 Mon Sep 17 00:00:00 2001 From: Samad Vakili Date: Tue, 9 Jun 2020 14:26:22 +0200 Subject: [PATCH 07/35] checked for typo --- python/tests/test_mechanics.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/python/tests/test_mechanics.py b/python/tests/test_mechanics.py index 0db6a1ba9..ae3cb391e 100755 --- a/python/tests/test_mechanics.py +++ b/python/tests/test_mechanics.py @@ -42,7 +42,7 @@ def Mises_stress(sigma): def PK2(P,F): - S = np.dot(_np.linalg.inv(F),P) + S = np.dot(np.linalg.inv(F),P) return symmetric(S) @@ -69,10 +69,10 @@ def strain_tensor(F,t,m): w,n = np.linalg.eigh(C) 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.einsum('ijk->ijk',np.eye(3))) 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.einsum('ijk->ijk',np.eye(3))) else: eps = np.matmul(n,np.einsum('ij,ikj->ijk',0.5*np.log(w),n)) @@ -157,7 +157,7 @@ class TestMechanics: @pytest.mark.parametrize('vectorized,single',[(mechanics.strain_tensor,strain_tensor)]) - def test_vectorize_strain_tensor(self,vectroized,single): + def test_vectorize_strain_tensor(self,vectorized,single): F = np.random.rand(self.n,3,3) F_vec = np.random.rand(self.n//10,10,3,3) t = ['V','U'][np.random.randint(0,2)] From 2978cf72f7ac2dd96e0caf48ad36e2095199c232 Mon Sep 17 00:00:00 2001 From: Samad Vakili Date: Tue, 9 Jun 2020 14:29:58 +0200 Subject: [PATCH 08/35] check for test_mechanics --- python/tests/test_mechanics.py | 0 1 file changed, 0 insertions(+), 0 deletions(-) mode change 100755 => 100644 python/tests/test_mechanics.py diff --git a/python/tests/test_mechanics.py b/python/tests/test_mechanics.py old mode 100755 new mode 100644 From d365cc9e12b7515b956d81a954c4f6f400ada566 Mon Sep 17 00:00:00 2001 From: Samad Vakili Date: Tue, 9 Jun 2020 14:39:17 +0200 Subject: [PATCH 09/35] typo --- python/tests/test_mechanics.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/tests/test_mechanics.py b/python/tests/test_mechanics.py index ae3cb391e..10a36fbcc 100644 --- a/python/tests/test_mechanics.py +++ b/python/tests/test_mechanics.py @@ -163,7 +163,7 @@ class TestMechanics: t = ['V','U'][np.random.randint(0,2)] m = np.random.random()*10.0 -5.0 for i,v in enumerate(np.reshape(vectorized(F_vec,t,m),vectorized(F,t,m).shape)): - assert np.allcloase(single(F[i],t,m),v) + assert np.allclose(single(F[i],t,m),v) @pytest.mark.parametrize('function',[mechanics.Cauchy, mechanics.PK2, From 694b7ec3c52be9b7d5ace8d70c96adb9bb4a7ab2 Mon Sep 17 00:00:00 2001 From: Samad Vakili Date: Tue, 9 Jun 2020 21:27:08 +0200 Subject: [PATCH 10/35] mechanics done --- python/damask/mechanics.py | 74 +++++++++++++++----------------------- 1 file changed, 29 insertions(+), 45 deletions(-) diff --git a/python/damask/mechanics.py b/python/damask/mechanics.py index c7ba39bcd..e94c22705 100644 --- a/python/damask/mechanics.py +++ b/python/damask/mechanics.py @@ -8,17 +8,13 @@ def Cauchy(P,F): Parameters ---------- - F : numpy.ndarray of shape (:,3,3) or (3,3) + F : numpy.ndarray of shape (...,3,3) Deformation gradient. - P : numpy.ndarray of shape (:,3,3) or (3,3) + P : numpy.ndarray of shape (...,3,3) First Piola-Kirchhoff stress. """ - if _np.shape(F) == _np.shape(P) == (3,3): - 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) + sigma = _np.einsum('...,...ij,...kj->...ik',1.0/_np.linalg.det(F),P,F) return symmetric(sigma) @@ -44,7 +40,7 @@ def eigenvalues(T_sym): 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. """ @@ -59,7 +55,7 @@ def eigenvectors(T_sym,RHS=False): 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. RHS: bool, optional 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)) if RHS: - if _np.shape(T_sym) == (3,3): - if _np.linalg.det(v) < 0.0: v[:,2] *= -1.0 - else: - v[_np.linalg.det(v) < 0.0,:,2] *= -1.0 + v[_np.linalg.det(v) < 0.0,:,2] *= -1.0 return v @@ -81,7 +74,7 @@ def left_stretch(T): 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. """ @@ -94,13 +87,12 @@ def maximum_shear(T_sym): 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. """ w = eigenvalues(T_sym) - return (w[0] - w[2])*0.5 if _np.shape(T_sym) == (3,3) else \ - (w[...,0] - w[...,2])*0.5 + return (w[...,0] - w[...,2])*0.5 def Mises_strain(epsilon): @@ -109,7 +101,7 @@ def Mises_strain(epsilon): 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. """ @@ -122,7 +114,7 @@ def Mises_stress(sigma): 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. """ @@ -135,16 +127,13 @@ def PK2(P,F): Parameters ---------- - P : numpy.ndarray of shape (...,3,3) or (3,3) + P : numpy.ndarray of shape (...,3,3) First Piola-Kirchhoff stress. - F : numpy.ndarray of shape (...,3,3) or (3,3) + F : numpy.ndarray of shape (...,3,3) Deformation gradient. """ - if _np.shape(F) == _np.shape(P) == (3,3): - S = _np.dot(_np.linalg.inv(F),P) - else: - S = _np.einsum('...jk,...kl->...jl',_np.linalg.inv(F),P) + S = _np.einsum('...jk,...kl->...jl',_np.linalg.inv(F),P) return symmetric(S) @@ -154,7 +143,7 @@ def right_stretch(T): 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. """ @@ -167,7 +156,7 @@ def rotational_part(T): 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. """ @@ -199,7 +188,7 @@ def strain_tensor(F,t,m): Parameters ---------- - F : numpy.ndarray of shape (:,3,3) or (3,3) + F : numpy.ndarray of shape (...,3,3) Deformation gradient. t : {‘V’, ‘U’} 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. """ - F_ = F.reshape(1,3,3) if F.shape == (3,3) else F if t == 'V': - B = _np.matmul(F_,transpose(F_)) + B = _np.matmul(F,transpose(F)) w,n = _np.linalg.eigh(B) elif t == 'U': - C = _np.matmul(transpose(F_),F_) + C = _np.matmul(transpose(F),F) w,n = _np.linalg.eigh(C) if m > 0.0: @@ -225,8 +213,7 @@ def strain_tensor(F,t,m): else: 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 \ - eps + return eps def symmetric(T): @@ -235,7 +222,7 @@ def symmetric(T): 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. """ @@ -248,12 +235,11 @@ def transpose(T): 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. """ - return T.T if _np.shape(T) == (3,3) else \ - _np.swapaxes(T,axis2=-2,axis1=-1) + return _np.swapaxes(T,axis2=-2,axis1=-1) def _polar_decomposition(T,requested): @@ -262,7 +248,7 @@ def _polar_decomposition(T,requested): 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. requested : iterable of str Requested outputs: ‘R’ for the rotation tensor, @@ -270,16 +256,15 @@ def _polar_decomposition(T,requested): """ u, s, vh = _np.linalg.svd(T) - R = _np.dot(u,vh) if _np.shape(T) == (3,3) else \ - _np.einsum('...ij,...jk->...ik',u,vh) + R = _np.einsum('...ij,...jk->...ik',u,vh) output = [] if 'R' in requested: output.append(R) 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: - 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) @@ -290,12 +275,11 @@ def _Mises(T_sym,s): 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. s : float Scaling factor (2/3 for strain, 3/2 for stress). """ d = deviatoric_part(T_sym) - return _np.sqrt(s*(_np.sum(d**2.0))) if _np.shape(T_sym) == (3,3) else \ - _np.sqrt(s*_np.einsum('...jk->...',d**2.0)) + return _np.sqrt(s*_np.einsum('...jk->...',d**2.0)) From 36f04309b4bbd906f6c57108b920329566ec8e6a Mon Sep 17 00:00:00 2001 From: Samad Vakili Date: Tue, 9 Jun 2020 21:33:13 +0200 Subject: [PATCH 11/35] issue in test_mechanics --- python/tests/test_mechanics.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/python/tests/test_mechanics.py b/python/tests/test_mechanics.py index 10a36fbcc..9708d4ee8 100644 --- a/python/tests/test_mechanics.py +++ b/python/tests/test_mechanics.py @@ -150,8 +150,8 @@ class TestMechanics: def test_vectorize_2_arg(self,vectorized,single): P = np.random.rand(self.n,3,3) F = np.random.rand(self.n,3,3) - P_vec = np.random.rand(self.n//10,10,3,3) - F_vec = np.random.rand(self.n//10,10,3,3) + P_vec = np.reshape(P,(self.n//10,10,3,3)) + F_vec = np.reshape(F,(self.n//10,10,3,3)) for i,v in enumerate(np.reshape(vectorized(P_vec,F_vec),vectorized(P,F).shape)): assert np.allclose(single(P[i],F[i]),v) @@ -159,7 +159,7 @@ class TestMechanics: @pytest.mark.parametrize('vectorized,single',[(mechanics.strain_tensor,strain_tensor)]) def test_vectorize_strain_tensor(self,vectorized,single): F = np.random.rand(self.n,3,3) - F_vec = np.random.rand(self.n//10,10,3,3) + F_vec = np.reshape(F,(self.n//10,10,3,3)) t = ['V','U'][np.random.randint(0,2)] m = np.random.random()*10.0 -5.0 for i,v in enumerate(np.reshape(vectorized(F_vec,t,m),vectorized(F,t,m).shape)): From f3d59ddfe56c2d80d00d907761c69545dbda3850 Mon Sep 17 00:00:00 2001 From: Samad Vakili Date: Wed, 10 Jun 2020 12:00:33 +0200 Subject: [PATCH 12/35] F reshaped error --- python/tests/test_mechanics.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/tests/test_mechanics.py b/python/tests/test_mechanics.py index 9708d4ee8..4f12ecf9b 100644 --- a/python/tests/test_mechanics.py +++ b/python/tests/test_mechanics.py @@ -62,10 +62,10 @@ def strain_tensor(F,t,m): F_ = F.reshape(1,3,3) if t == 'V': - B = np.matmul(F_,transpose(F_)) + B = np.matmul(F_,mechanics.transpose(F_)) w,n = np.linalg.eigh(B) elif t == 'U': - C = np.matmul(transpose(F_),F_) + C = np.matmul(mechanics.transpose(F_),F_) w,n = np.linalg.eigh(C) if m > 0.0: From 7d7927e821362083a195b208d76657e62054e8ba Mon Sep 17 00:00:00 2001 From: Samad Vakili Date: Wed, 10 Jun 2020 13:22:57 +0200 Subject: [PATCH 13/35] reshaped F transpose error --- python/tests/test_mechanics.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/tests/test_mechanics.py b/python/tests/test_mechanics.py index 4f12ecf9b..0ae433254 100644 --- a/python/tests/test_mechanics.py +++ b/python/tests/test_mechanics.py @@ -62,10 +62,10 @@ def strain_tensor(F,t,m): F_ = F.reshape(1,3,3) if t == 'V': - B = np.matmul(F_,mechanics.transpose(F_)) + B = np.matmul(F_,F_[1].T) w,n = np.linalg.eigh(B) elif t == 'U': - C = np.matmul(mechanics.transpose(F_),F_) + C = np.matmul(F_[1].T,F_) w,n = np.linalg.eigh(C) if m > 0.0: From 4f4d8d8c926ab612b46bcadb735129961c51b305 Mon Sep 17 00:00:00 2001 From: Samad Vakili Date: Wed, 10 Jun 2020 13:27:36 +0200 Subject: [PATCH 14/35] reshaped F transpose error2 --- python/tests/test_mechanics.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/tests/test_mechanics.py b/python/tests/test_mechanics.py index 0ae433254..c740ed6f3 100644 --- a/python/tests/test_mechanics.py +++ b/python/tests/test_mechanics.py @@ -62,10 +62,10 @@ def strain_tensor(F,t,m): F_ = F.reshape(1,3,3) if t == 'V': - B = np.matmul(F_,F_[1].T) + B = np.matmul(F_,F_[0].T) w,n = np.linalg.eigh(B) elif t == 'U': - C = np.matmul(F_[1].T,F_) + C = np.matmul(F_[0].T,F_) w,n = np.linalg.eigh(C) if m > 0.0: From 0f42ff6b5dd7ac98b21f0533ebd1736fd6ca7b8b Mon Sep 17 00:00:00 2001 From: Samad Vakili Date: Wed, 10 Jun 2020 14:02:52 +0200 Subject: [PATCH 15/35] reshaped F transpose error3 --- python/tests/test_mechanics.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/python/tests/test_mechanics.py b/python/tests/test_mechanics.py index c740ed6f3..cba5a206c 100644 --- a/python/tests/test_mechanics.py +++ b/python/tests/test_mechanics.py @@ -62,18 +62,18 @@ def strain_tensor(F,t,m): F_ = F.reshape(1,3,3) if t == 'V': - B = np.matmul(F_,F_[0].T) + B = np.matmul(F_,mechanics.transpose(F_)) w,n = np.linalg.eigh(B) elif t == 'U': - C = np.matmul(F_[0].T,F_) + C = np.matmul(mechanics.transpose(F_),F_) w,n = np.linalg.eigh(C) if m > 0.0: eps = 1.0/(2.0*abs(m)) * (+ np.matmul(n,np.einsum('ij,ikj->ijk',w**m,n)) - - np.einsum('ijk->ijk',np.eye(3))) + - np.broadcast_to(np.eye(3),[F_.shape[0],3,3])) elif m < 0.0: eps = 1.0/(2.0*abs(m)) * (- np.matmul(n,np.einsum('ij,ikj->ijk',w**m,n)) - + np.einsum('ijk->ijk',np.eye(3))) + + np.broadcast_to(np.eye(3),[F_.shape[0],3,3])) else: eps = np.matmul(n,np.einsum('ij,ikj->ijk',0.5*np.log(w),n)) From 6f81f5278d04796b843cacdd22c39e499eeb09d4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 15 Nov 2020 23:14:46 +0100 Subject: [PATCH 16/35] separating general tensor math from mechanics operations --- python/damask/__init__.py | 3 +- python/damask/_orientation.py | 4 +- python/damask/_result.py | 5 +- python/damask/_rotation.py | 4 +- python/damask/mechanics.py | 80 ++++--------------------- python/damask/seeds.py | 5 ++ python/damask/tensor.py | 74 +++++++++++++++++++++++ python/tests/test_Result.py | 7 ++- python/tests/test_mechanics.py | 105 +++++++++------------------------ python/tests/test_tensor.py | 73 +++++++++++++++++++++++ 10 files changed, 203 insertions(+), 157 deletions(-) create mode 100644 python/damask/tensor.py create mode 100644 python/tests/test_tensor.py diff --git a/python/damask/__init__.py b/python/damask/__init__.py index 9d0350337..f3d5ad2fb 100644 --- a/python/damask/__init__.py +++ b/python/damask/__init__.py @@ -7,15 +7,16 @@ with open(_Path(__file__).parent/_Path('VERSION')) as _f: version = _re.sub(r'^v','',_f.readline().strip()) __version__ = version -# make classes directly accessible as damask.Class from ._environment import Environment as _ # noqa environment = _() from . import util # noqa from . import seeds # noqa +from . import tensor # noqa from . import mechanics # noqa from . import solver # noqa from . import grid_filters # noqa from . import lattice # noqa +# make classes directly accessible as damask.Class from ._rotation import Rotation # noqa from ._orientation import Orientation # noqa from ._table import Table # noqa diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index e0a00e0a0..9ae1d2ed3 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -2,7 +2,7 @@ import numpy as np from . import Rotation from . import util -from . import mechanics +from . import tensor __parameter_doc__ = \ """lattice : str @@ -362,7 +362,7 @@ class Orientation(Rotation): x = o.to_frame(uvw=uvw) z = o.to_frame(hkl=hkl) om = np.stack([x,np.cross(z,x),z],axis=-2) - return o.copy(rotation=Rotation.from_matrix(mechanics.transpose(om/np.linalg.norm(om,axis=-1,keepdims=True)))) + return o.copy(rotation=Rotation.from_matrix(tensor.transpose(om/np.linalg.norm(om,axis=-1,keepdims=True)))) @property diff --git a/python/damask/_result.py b/python/damask/_result.py index 5bed4347f..a5dcca8db 100644 --- a/python/damask/_result.py +++ b/python/damask/_result.py @@ -18,6 +18,7 @@ from . import Table from . import Orientation from . import grid_filters from . import mechanics +from . import tensor from . import util h5py3 = h5py.__version__[0] == '3' @@ -681,7 +682,7 @@ class Result: label,p = 'Minimum',0 return { - 'data': mechanics.eigenvalues(T_sym['data'])[:,p], + 'data': tensor.eigenvalues(T_sym['data'])[:,p], 'label': f"lambda_{eigenvalue}({T_sym['label']})", 'meta' : { 'Unit': T_sym['meta']['Unit'], @@ -713,7 +714,7 @@ class Result: elif eigenvalue == 'min': label,p = 'minimum',0 return { - 'data': mechanics.eigenvectors(T_sym['data'])[:,p], + 'data': tensor.eigenvectors(T_sym['data'])[:,p], 'label': f"v_{eigenvalue}({T_sym['label']})", 'meta' : { 'Unit': '1', diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index deae8b6ac..e60b14ecb 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -1,6 +1,6 @@ import numpy as np -from . import mechanics +from . import tensor from . import util from . import grid_filters @@ -549,7 +549,7 @@ class Rotation: raise ValueError('Invalid shape.') if reciprocal: - om = np.linalg.inv(mechanics.transpose(om)/np.pi) # transform reciprocal basis set + om = np.linalg.inv(tensor.transpose(om)/np.pi) # transform reciprocal basis set orthonormal = False # contains stretch if not orthonormal: (U,S,Vh) = np.linalg.svd(om) # singular value decomposition diff --git a/python/damask/mechanics.py b/python/damask/mechanics.py index e94c22705..dca26b9de 100644 --- a/python/damask/mechanics.py +++ b/python/damask/mechanics.py @@ -1,5 +1,10 @@ +"""Finite-strain continuum mechanics.""" + +from . import tensor + import numpy as _np + def Cauchy(P,F): """ Return Cauchy stress calculated from first Piola-Kirchhoff stress and deformation gradient. @@ -15,7 +20,7 @@ def Cauchy(P,F): """ sigma = _np.einsum('...,...ij,...kj->...ik',1.0/_np.linalg.det(F),P,F) - return symmetric(sigma) + return tensor.symmetric(sigma) def deviatoric_part(T): @@ -31,43 +36,6 @@ def deviatoric_part(T): return T - _np.einsum('...ij,...->...ij',_np.eye(3),spherical_part(T)) -def eigenvalues(T_sym): - """ - Return the eigenvalues, i.e. principal components, of a symmetric tensor. - - The eigenvalues are sorted in ascending order, each repeated according to - its multiplicity. - - Parameters - ---------- - T_sym : numpy.ndarray of shape (...,3,3) - Symmetric tensor of which the eigenvalues are computed. - - """ - return _np.linalg.eigvalsh(symmetric(T_sym)) - - -def eigenvectors(T_sym,RHS=False): - """ - Return eigenvectors of a symmetric tensor. - - The eigenvalues are sorted in ascending order of their associated eigenvalues. - - Parameters - ---------- - T_sym : numpy.ndarray of shape (...,3,3) - Symmetric tensor of which the eigenvectors are computed. - RHS: bool, optional - Enforce right-handed coordinate system. Default is False. - - """ - (u,v) = _np.linalg.eigh(symmetric(T_sym)) - - if RHS: - v[_np.linalg.det(v) < 0.0,:,2] *= -1.0 - return v - - def left_stretch(T): """ Return the left stretch of a tensor. @@ -91,7 +59,7 @@ def maximum_shear(T_sym): Symmetric tensor of which the maximum shear is computed. """ - w = eigenvalues(T_sym) + w = tensor.eigenvalues(T_sym) return (w[...,0] - w[...,2])*0.5 @@ -134,7 +102,7 @@ def PK2(P,F): """ S = _np.einsum('...jk,...kl->...jl',_np.linalg.inv(F),P) - return symmetric(S) + return tensor.symmetric(S) def right_stretch(T): @@ -197,16 +165,16 @@ def strain_tensor(F,t,m): """ if t == 'V': - B = _np.matmul(F,transpose(F)) + B = _np.matmul(F,tensor.transpose(F)) w,n = _np.linalg.eigh(B) elif t == 'U': - C = _np.matmul(transpose(F),F) + C = _np.matmul(tensor.transpose(F),F) w,n = _np.linalg.eigh(C) if m > 0.0: eps = 1.0/(2.0*abs(m)) * (+ _np.matmul(n,_np.einsum('...j,...kj->...jk',w**m,n)) - _np.einsum('...jk->...jk',_np.eye(3))) - + elif m < 0.0: eps = 1.0/(2.0*abs(m)) * (- _np.matmul(n,_np.einsum('...j,...kj->...jk',w**m,n)) + _np.einsum('...jk->...jk',_np.eye(3))) @@ -216,32 +184,6 @@ def strain_tensor(F,t,m): return eps -def symmetric(T): - """ - Return the symmetrized tensor. - - Parameters - ---------- - T : numpy.ndarray of shape (...,3,3) - Tensor of which the symmetrized values are computed. - - """ - return (T+transpose(T))*0.5 - - -def transpose(T): - """ - Return the transpose of a tensor. - - Parameters - ---------- - T : numpy.ndarray of shape (...,3,3) - Tensor of which the transpose is computed. - - """ - return _np.swapaxes(T,axis2=-2,axis1=-1) - - def _polar_decomposition(T,requested): """ Singular value decomposition. diff --git a/python/damask/seeds.py b/python/damask/seeds.py index 9aab953d0..c75dfb5e5 100644 --- a/python/damask/seeds.py +++ b/python/damask/seeds.py @@ -1,3 +1,8 @@ +""" +Functionality for generation of seed points for Voronoi +or Laguerre tessellation. +""" + from scipy import spatial as _spatial import numpy as _np diff --git a/python/damask/tensor.py b/python/damask/tensor.py new file mode 100644 index 000000000..71617e32c --- /dev/null +++ b/python/damask/tensor.py @@ -0,0 +1,74 @@ +""" +Tensor operations. + +Notes +----- +This is not a tensor class, but a collection of routines +to operate on numpy.ndarrays of shape (...,3,3). + +""" + +import numpy as _np + + +def symmetric(T): + """ + Return the symmetrized tensor. + + Parameters + ---------- + T : numpy.ndarray of shape (...,3,3) + Tensor of which the symmetrized values are computed. + + """ + return (T+transpose(T))*0.5 + + +def transpose(T): + """ + Return the transpose of a tensor. + + Parameters + ---------- + T : numpy.ndarray of shape (...,3,3) + Tensor of which the transpose is computed. + + """ + return _np.swapaxes(T,axis2=-2,axis1=-1) + + +def eigenvalues(T_sym): + """ + Return the eigenvalues, i.e. principal components, of a symmetric tensor. + + The eigenvalues are sorted in ascending order, each repeated according to + its multiplicity. + + Parameters + ---------- + T_sym : numpy.ndarray of shape (...,3,3) + Symmetric tensor of which the eigenvalues are computed. + + """ + return _np.linalg.eigvalsh(symmetric(T_sym)) + + +def eigenvectors(T_sym,RHS=False): + """ + Return eigenvectors of a symmetric tensor. + + The eigenvalues are sorted in ascending order of their associated eigenvalues. + + Parameters + ---------- + T_sym : numpy.ndarray of shape (...,3,3) + Symmetric tensor of which the eigenvectors are computed. + RHS: bool, optional + Enforce right-handed coordinate system. Default is False. + + """ + (u,v) = _np.linalg.eigh(symmetric(T_sym)) + + if RHS: + v[_np.linalg.det(v) < 0.0,:,2] *= -1.0 + return v diff --git a/python/tests/test_Result.py b/python/tests/test_Result.py index 929c0ef44..b739bb2ca 100644 --- a/python/tests/test_Result.py +++ b/python/tests/test_Result.py @@ -11,6 +11,7 @@ import h5py from damask import Result from damask import Rotation from damask import Orientation +from damask import tensor from damask import mechanics from damask import grid_filters @@ -152,7 +153,7 @@ class TestResult: default.add_eigenvalue('sigma',eigenvalue) loc = {'sigma' :default.get_dataset_location('sigma'), 'lambda':default.get_dataset_location(f'lambda_{eigenvalue}(sigma)')} - in_memory = function(mechanics.eigenvalues(default.read_dataset(loc['sigma'],0)),axis=1,keepdims=True) + in_memory = function(tensor.eigenvalues(default.read_dataset(loc['sigma'],0)),axis=1,keepdims=True) in_file = default.read_dataset(loc['lambda'],0) assert np.allclose(in_memory,in_file) @@ -162,7 +163,7 @@ class TestResult: default.add_eigenvector('sigma',eigenvalue) loc = {'sigma' :default.get_dataset_location('sigma'), 'v(sigma)':default.get_dataset_location(f'v_{eigenvalue}(sigma)')} - in_memory = mechanics.eigenvectors(default.read_dataset(loc['sigma'],0))[:,idx] + in_memory = tensor.eigenvectors(default.read_dataset(loc['sigma'],0))[:,idx] in_file = default.read_dataset(loc['v(sigma)'],0) assert np.allclose(in_memory,in_file) @@ -210,7 +211,7 @@ class TestResult: in_memory = mechanics.Mises_stress(default.read_dataset(loc['sigma'],0)).reshape(-1,1) in_file = default.read_dataset(loc['sigma_vM'],0) assert np.allclose(in_memory,in_file) - + def test_add_Mises_invalid(self,default): default.add_Cauchy('P','F') default.add_calculation('sigma_y','#sigma#',unit='y') diff --git a/python/tests/test_mechanics.py b/python/tests/test_mechanics.py index 9bcf6c6dc..4c0d68572 100644 --- a/python/tests/test_mechanics.py +++ b/python/tests/test_mechanics.py @@ -1,11 +1,12 @@ import pytest import numpy as np +from damask import tensor from damask import mechanics def Cauchy(P,F): sigma = 1.0/np.linalg.det(F) * np.dot(P,F.T) - return mechanics.symmetric(sigma) + return symmetric(sigma) def deviatoric_part(T): @@ -16,14 +17,6 @@ def eigenvalues(T_sym): return np.linalg.eigvalsh(symmetric(T_sym)) -def eigenvectors(T_sym,RHS=False): - (u,v) = np.linalg.eigh(symmetric(T_sym)) - - if RHS: - if np.linalg.det(v) < 0.0: v[:,2] *= -1.0 - return v - - def left_stretch(T): return polar_decomposition(T,'V')[0] @@ -53,39 +46,35 @@ def right_stretch(T): def rotational_part(T): return polar_decomposition(T,'R')[0] + def spherical_part(T,tensor=False): sph = np.trace(T)/3.0 return sph if not tensor else np.eye(3)*sph def strain_tensor(F,t,m): - F_ = F.reshape(1,3,3) - + if t == 'V': - B = np.matmul(F_,mechanics.transpose(F_)) + B = np.matmul(F,F.T) w,n = np.linalg.eigh(B) elif t == 'U': - C = np.matmul(mechanics.transpose(F_),F_) + C = np.matmul(F.T,F) w,n = np.linalg.eigh(C) if m > 0.0: - eps = 1.0/(2.0*abs(m)) * (+ np.matmul(n,np.einsum('ij,ikj->ijk',w**m,n)) - - np.broadcast_to(np.eye(3),[F_.shape[0],3,3])) + eps = 1.0/(2.0*abs(m)) * (+ np.matmul(n,np.einsum('j,kj->jk',w**m,n)) + - np.eye(3)) elif m < 0.0: - eps = 1.0/(2.0*abs(m)) * (- np.matmul(n,np.einsum('ij,ikj->ijk',w**m,n)) - + np.broadcast_to(np.eye(3),[F_.shape[0],3,3])) + eps = 1.0/(2.0*abs(m)) * (- np.matmul(n,np.einsum('j,kj->jk',w**m,n)) + + np.eye(3)) else: - eps = np.matmul(n,np.einsum('ij,ikj->ijk',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) + return eps def symmetric(T): - return (T+transpose(T))*0.5 - - -def transpose(T): - return T.T + return (T+T.T)*0.5 def polar_decomposition(T,requested): @@ -96,7 +85,7 @@ def polar_decomposition(T,requested): if 'R' in requested: output.append(R) if 'V' in requested: - output.append(np.dot(T,R.T)) + output.append(np.dot(T,R.T)) if 'U' in requested: output.append(np.dot(R.T,T)) @@ -123,20 +112,15 @@ class TestMechanics: for i,v in enumerate(np.reshape(vectorized(test_data),vectorized(test_data_flat).shape)): assert np.allclose(single(test_data_flat[i]),v) - @pytest.mark.parametrize('vectorized,single',[ - (mechanics.deviatoric_part, deviatoric_part), - (mechanics.eigenvalues , eigenvalues ), - (mechanics.eigenvectors , eigenvectors ), - (mechanics.left_stretch , left_stretch ), - (mechanics.maximum_shear , maximum_shear ), - (mechanics.Mises_strain , Mises_strain ), - (mechanics.Mises_stress , Mises_stress ), - (mechanics.right_stretch , right_stretch ), - (mechanics.rotational_part, rotational_part), - (mechanics.spherical_part , spherical_part ), - (mechanics.symmetric , symmetric ), - (mechanics.transpose , transpose ), - ]) + @pytest.mark.parametrize('vectorized,single',[(mechanics.deviatoric_part, deviatoric_part), + (mechanics.left_stretch , left_stretch ), + (mechanics.maximum_shear , maximum_shear ), + (mechanics.Mises_strain , Mises_strain ), + (mechanics.Mises_stress , Mises_stress ), + (mechanics.right_stretch , right_stretch ), + (mechanics.rotational_part, rotational_part), + (mechanics.spherical_part , spherical_part ), + ]) def test_vectorize_1_arg(self,vectorized,single): epsilon = np.random.rand(self.n,3,3) epsilon_vec = np.reshape(epsilon,(self.n//10,10,3,3)) @@ -171,7 +155,7 @@ class TestMechanics: def test_stress_measures(self,function): """Ensure that all stress measures are equivalent for no deformation.""" P = np.random.rand(self.n,3,3) - assert np.allclose(function(P,np.broadcast_to(np.eye(3),(self.n,3,3))),mechanics.symmetric(P)) + assert np.allclose(function(P,np.broadcast_to(np.eye(3),(self.n,3,3))),tensor.symmetric(P)) def test_deviatoric_part(self): I_n = np.broadcast_to(np.eye(3),(self.n,3,3)) @@ -237,9 +221,9 @@ class TestMechanics: def test_spherical_mapping(self): """Ensure that mapping to tensor is correct.""" x = np.random.rand(self.n,3,3) - tensor = mechanics.spherical_part(x,True) + tnsr = mechanics.spherical_part(x,True) scalar = mechanics.spherical_part(x) - assert np.allclose(np.linalg.det(tensor), + assert np.allclose(np.linalg.det(tnsr), scalar**3.0) def test_spherical_Mises(self): @@ -249,17 +233,6 @@ class TestMechanics: assert np.allclose(mechanics.Mises_strain(sph), 0.0) - def test_symmetric(self): - """Ensure that a symmetric tensor is half of the sum of a tensor and its transpose.""" - x = np.random.rand(self.n,3,3) - assert np.allclose(mechanics.symmetric(x)*2.0, - mechanics.transpose(x)+x) - - def test_transpose(self): - """Ensure that a symmetric tensor equals its transpose.""" - x = mechanics.symmetric(np.random.rand(self.n,3,3)) - assert np.allclose(mechanics.transpose(x), - x) def test_Mises(self): """Ensure that equivalent stress is 3/2 of equivalent strain.""" @@ -267,31 +240,7 @@ class TestMechanics: assert np.allclose(mechanics.Mises_stress(x)/mechanics.Mises_strain(x), 1.5) - def test_eigenvalues(self): - """Ensure that the characteristic polynomial can be solved.""" - A = mechanics.symmetric(np.random.rand(self.n,3,3)) - lambd = mechanics.eigenvalues(A) - s = np.random.randint(self.n) - for i in range(3): - assert np.allclose(np.linalg.det(A[s]-lambd[s,i]*np.eye(3)),.0) - - def test_eigenvalues_and_vectors(self): - """Ensure that eigenvalues and -vectors are the solution to the characteristic polynomial.""" - A = mechanics.symmetric(np.random.rand(self.n,3,3)) - lambd = mechanics.eigenvalues(A) - x = mechanics.eigenvectors(A) - s = np.random.randint(self.n) - for i in range(3): - assert np.allclose(np.dot(A[s]-lambd[s,i]*np.eye(3),x[s,:,i]),.0) - - def test_eigenvectors_RHS(self): - """Ensure that RHS coordinate system does only change sign of determinant.""" - A = mechanics.symmetric(np.random.rand(self.n,3,3)) - LRHS = np.linalg.det(mechanics.eigenvectors(A,RHS=False)) - RHS = np.linalg.det(mechanics.eigenvectors(A,RHS=True)) - assert np.allclose(np.abs(LRHS),RHS) - def test_spherical_no_shear(self): """Ensure that sherical stress has max shear of 0.0.""" - A = mechanics.spherical_part(mechanics.symmetric(np.random.rand(self.n,3,3)),True) + A = mechanics.spherical_part(tensor.symmetric(np.random.rand(self.n,3,3)),True) assert np.allclose(mechanics.maximum_shear(A),0.0) diff --git a/python/tests/test_tensor.py b/python/tests/test_tensor.py new file mode 100644 index 000000000..24fbe6126 --- /dev/null +++ b/python/tests/test_tensor.py @@ -0,0 +1,73 @@ +import pytest +import numpy as np + +from damask import tensor + + +def eigenvalues(T_sym): + return np.linalg.eigvalsh(symmetric(T_sym)) + +def eigenvectors(T_sym,RHS=False): + (u,v) = np.linalg.eigh(symmetric(T_sym)) + + if RHS: + if np.linalg.det(v) < 0.0: v[:,2] *= -1.0 + return v + +def symmetric(T): + return (T+transpose(T))*0.5 + +def transpose(T): + return T.T + + +class TestTensor: + + n = 1000 + c = np.random.randint(n) + + + @pytest.mark.parametrize('vectorized,single',[(tensor.eigenvalues , eigenvalues ), + (tensor.eigenvectors , eigenvectors ), + (tensor.symmetric , symmetric ), + (tensor.transpose , transpose ), + ]) + def test_vectorize_1_arg(self,vectorized,single): + epsilon = np.random.rand(self.n,3,3) + epsilon_vec = np.reshape(epsilon,(self.n//10,10,3,3)) + for i,v in enumerate(np.reshape(vectorized(epsilon_vec),vectorized(epsilon).shape)): + assert np.allclose(single(epsilon[i]),v) + + def test_symmetric(self): + """Ensure that a symmetric tensor is half of the sum of a tensor and its transpose.""" + x = np.random.rand(self.n,3,3) + assert np.allclose(tensor.symmetric(x)*2.0,tensor.transpose(x)+x) + + def test_transpose(self): + """Ensure that a symmetric tensor equals its transpose.""" + x = tensor.symmetric(np.random.rand(self.n,3,3)) + assert np.allclose(tensor.transpose(x),x) + + def test_eigenvalues(self): + """Ensure that the characteristic polynomial can be solved.""" + A = tensor.symmetric(np.random.rand(self.n,3,3)) + lambd = tensor.eigenvalues(A) + s = np.random.randint(self.n) + for i in range(3): + assert np.allclose(np.linalg.det(A[s]-lambd[s,i]*np.eye(3)),.0) + + def test_eigenvalues_and_vectors(self): + """Ensure that eigenvalues and -vectors are the solution to the characteristic polynomial.""" + A = tensor.symmetric(np.random.rand(self.n,3,3)) + lambd = tensor.eigenvalues(A) + x = tensor.eigenvectors(A) + s = np.random.randint(self.n) + for i in range(3): + assert np.allclose(np.dot(A[s]-lambd[s,i]*np.eye(3),x[s,:,i]),.0) + + def test_eigenvectors_RHS(self): + """Ensure that RHS coordinate system does only change sign of determinant.""" + A = tensor.symmetric(np.random.rand(self.n,3,3)) + LRHS = np.linalg.det(tensor.eigenvectors(A,RHS=False)) + RHS = np.linalg.det(tensor.eigenvectors(A,RHS=True)) + assert np.allclose(np.abs(LRHS),RHS) From b893967b68a828e5f8fee8475cd2ffab06fdb965 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 16 Nov 2020 01:01:32 +0100 Subject: [PATCH 17/35] more systematic names and extended docstrings --- python/damask/_result.py | 2 +- python/damask/mechanics.py | 186 +++++++++++++++++++++++++-------- python/damask/tensor.py | 37 +++++-- python/tests/test_Result.py | 4 +- python/tests/test_mechanics.py | 24 ++--- 5 files changed, 182 insertions(+), 71 deletions(-) diff --git a/python/damask/_result.py b/python/damask/_result.py index a5dcca8db..ad8c5ec0e 100644 --- a/python/damask/_result.py +++ b/python/damask/_result.py @@ -1024,7 +1024,7 @@ class Result: @staticmethod def _add_stretch_tensor(F,t): return { - 'data': (mechanics.left_stretch if t.upper() == 'V' else mechanics.right_stretch)(F['data']), + 'data': (mechanics.stretch_left if t.upper() == 'V' else mechanics.stretch_right)(F['data']), 'label': f"{t}({F['label']})", 'meta': { 'Unit': F['meta']['Unit'], diff --git a/python/damask/mechanics.py b/python/damask/mechanics.py index dca26b9de..9118ea875 100644 --- a/python/damask/mechanics.py +++ b/python/damask/mechanics.py @@ -5,18 +5,58 @@ from . import tensor import numpy as _np -def Cauchy(P,F): +def Cauchy_Green_deformation_left(F): """ - Return Cauchy stress calculated from first Piola-Kirchhoff stress and deformation gradient. - - Resulting tensor is symmetrized as the Cauchy stress needs to be symmetric. + Calculate left Cauchy-Green deformation tensor (Finger deformation tensor). Parameters ---------- F : numpy.ndarray of shape (...,3,3) Deformation gradient. + + Returns + ------- + B : numpy.ndarray of shape (...,3,3) + Left Cauchy-Green deformation tensor. + + """ + return _np.matmul(F,tensor.transpose(F)) + + +def Cauchy_Green_deformation_right(F): + """ + Calculate right Cauchy-Green deformation tensor. + + Parameters + ---------- + F : numpy.ndarray of shape (...,3,3) + Deformation gradient. + + Returns + ------- + C : numpy.ndarray of shape (...,3,3) + Right Cauchy-Green deformation tensor. + + """ + return _np.matmul(tensor.transpose(F),F) + +def Cauchy(P,F): + """ + Calculate the Cauchy (true) stress. + + Resulting tensor is symmetrized as the Cauchy stress needs to be symmetric. + + Parameters + ---------- P : numpy.ndarray of shape (...,3,3) First Piola-Kirchhoff stress. + F : numpy.ndarray of shape (...,3,3) + Deformation gradient. + + Returns + ------- + sigma : numpy.ndarray of shape (...,3,3) + Cauchy stress. """ sigma = _np.einsum('...,...ij,...kj->...ik',1.0/_np.linalg.det(F),P,F) @@ -25,39 +65,36 @@ def Cauchy(P,F): def deviatoric_part(T): """ - Return deviatoric part of a tensor. + Calculate deviatoric part of a tensor. Parameters ---------- T : numpy.ndarray of shape (...,3,3) Tensor of which the deviatoric part is computed. + Returns + ------- + T' : numpy.ndarray of shape (...,3,3) + Deviatoric part of T. + """ return T - _np.einsum('...ij,...->...ij',_np.eye(3),spherical_part(T)) -def left_stretch(T): - """ - Return the left stretch of a tensor. - - Parameters - ---------- - T : numpy.ndarray of shape (...,3,3) - Tensor of which the left stretch is computed. - - """ - return _polar_decomposition(T,'V')[0] - - def maximum_shear(T_sym): """ - Return the maximum shear component of a symmetric tensor. + Calculate the maximum shear component of a symmetric tensor. Parameters ---------- T_sym : numpy.ndarray of shape (...,3,3) Symmetric tensor of which the maximum shear is computed. + Returns + ------- + gamma_max : numpy.ndarray of shape (...) + Maximum shear of T_sym. + """ w = tensor.eigenvalues(T_sym) return (w[...,0] - w[...,2])*0.5 @@ -65,33 +102,46 @@ def maximum_shear(T_sym): def Mises_strain(epsilon): """ - Return the Mises equivalent of a strain tensor. + Calculate the Mises equivalent of a strain tensor. Parameters ---------- epsilon : numpy.ndarray of shape (...,3,3) Symmetric strain tensor of which the von Mises equivalent is computed. + Returns + ------- + epsilon_vM : numpy.ndarray of shape (...) + Von Mises equivalent strain of epsilon. + """ return _Mises(epsilon,2.0/3.0) def Mises_stress(sigma): """ - Return the Mises equivalent of a stress tensor. + Calculate the Mises equivalent of a stress tensor. Parameters ---------- sigma : numpy.ndarray of shape (...,3,3) Symmetric stress tensor of which the von Mises equivalent is computed. + Returns + ------- + sigma_vM : numpy.ndarray of shape (...) + Von Mises equivalent stress of sigma. + """ return _Mises(sigma,3.0/2.0) def PK2(P,F): """ - Calculate second Piola-Kirchhoff stress from first Piola-Kirchhoff stress and deformation gradient. + Calculate the second Piola-Kirchhoff stress. + + Resulting tensor is symmetrized as the second Piola-Kirchhoff stress + needs to be symmetric. Parameters ---------- @@ -100,47 +150,51 @@ def PK2(P,F): F : numpy.ndarray of shape (...,3,3) Deformation gradient. + Returns + ------- + S : numpy.ndarray of shape (...,3,3) + Second Piola-Kirchhoff stress. + """ S = _np.einsum('...jk,...kl->...jl',_np.linalg.inv(F),P) return tensor.symmetric(S) -def right_stretch(T): - """ - Return the right stretch of a tensor. - - Parameters - ---------- - T : numpy.ndarray of shape (...,3,3) - Tensor of which the right stretch is computed. - - """ - return _polar_decomposition(T,'U')[0] - - def rotational_part(T): """ - Return the rotational part of a tensor. + Calculate the rotational part of a tensor. Parameters ---------- T : numpy.ndarray of shape (...,3,3) Tensor of which the rotational part is computed. + Returns + ------- + R : numpy.ndarray of shape (...,3,3) + Rotational part. + """ return _polar_decomposition(T,'R')[0] def spherical_part(T,tensor=False): """ - Return spherical (hydrostatic) part of a tensor. + Calculate spherical (hydrostatic) part of a tensor. Parameters ---------- T : numpy.ndarray of shape (...,3,3) Tensor of which the hydrostatic part is computed. tensor : bool, optional - Map spherical part onto identity tensor. Default is false + Map spherical part onto identity tensor. Defaults to false + + Returns + ------- + p : numpy.ndarray of shape (...) + unless tensor == True: shape (...,3,3) + Spherical part of tensor T, e.g. the hydrostatic part/pressure + of a stress tensor. """ sph = _np.trace(T,axis2=-2,axis1=-1)/3.0 @@ -149,7 +203,7 @@ def spherical_part(T,tensor=False): def strain_tensor(F,t,m): """ - Return strain tensor calculated from deformation gradient. + Calculate strain tensor from deformation gradient. For details refer to https://en.wikipedia.org/wiki/Finite_strain_theory and https://de.wikipedia.org/wiki/Verzerrungstensor @@ -159,17 +213,21 @@ def strain_tensor(F,t,m): F : numpy.ndarray of shape (...,3,3) Deformation gradient. 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 Order of the strain. + Returns + ------- + epsilon : numpy.ndarray of shape (...,3,3) + Strain of F. + """ if t == 'V': - B = _np.matmul(F,tensor.transpose(F)) - w,n = _np.linalg.eigh(B) + w,n = _np.linalg.eigh(Cauchy_Green_deformation_left(F)) elif t == 'U': - C = _np.matmul(tensor.transpose(F),F) - w,n = _np.linalg.eigh(C) + w,n = _np.linalg.eigh(Cauchy_Green_deformation_right(F)) if m > 0.0: eps = 1.0/(2.0*abs(m)) * (+ _np.matmul(n,_np.einsum('...j,...kj->...jk',w**m,n)) @@ -184,9 +242,45 @@ def strain_tensor(F,t,m): return eps +def stretch_left(T): + """ + Calculate left stretch of a tensor. + + Parameters + ---------- + T : numpy.ndarray of shape (...,3,3) + Tensor of which the left stretch is computed. + + Returns + ------- + V : numpy.ndarray of shape (...,3,3) + Left stretch tensor from Polar decomposition of T. + + """ + return _polar_decomposition(T,'V')[0] + + +def stretch_right(T): + """ + Calculate right stretch of a tensor. + + Parameters + ---------- + T : numpy.ndarray of shape (...,3,3) + Tensor of which the right stretch is computed. + + Returns + ------- + U : numpy.ndarray of shape (...,3,3) + Left stretch tensor from Polar decomposition of T. + + """ + return _polar_decomposition(T,'U')[0] + + def _polar_decomposition(T,requested): """ - Singular value decomposition. + Perform singular value decomposition. Parameters ---------- @@ -197,7 +291,7 @@ def _polar_decomposition(T,requested): ‘V’ for left stretch tensor and ‘U’ for right stretch tensor. """ - u, s, vh = _np.linalg.svd(T) + u, _, vh = _np.linalg.svd(T) R = _np.einsum('...ij,...jk->...ik',u,vh) output = [] diff --git a/python/damask/tensor.py b/python/damask/tensor.py index 71617e32c..6e5bdaf85 100644 --- a/python/damask/tensor.py +++ b/python/damask/tensor.py @@ -13,58 +13,75 @@ import numpy as _np def symmetric(T): """ - Return the symmetrized tensor. + Symmetrize tensor. Parameters ---------- T : numpy.ndarray of shape (...,3,3) Tensor of which the symmetrized values are computed. + Returns + ------- + T_sym : numpy.ndarray of shape (...,3,3) + Symmetrized tensor T. + """ return (T+transpose(T))*0.5 def transpose(T): """ - Return the transpose of a tensor. + Transpose tensor. Parameters ---------- T : numpy.ndarray of shape (...,3,3) Tensor of which the transpose is computed. + Returns + ------- + T.T : numpy.ndarray of shape (...,3,3) + Transpose of tensor T. + """ return _np.swapaxes(T,axis2=-2,axis1=-1) def eigenvalues(T_sym): """ - Return the eigenvalues, i.e. principal components, of a symmetric tensor. - - The eigenvalues are sorted in ascending order, each repeated according to - its multiplicity. + Eigenvalues, i.e. principal components, of a symmetric tensor. Parameters ---------- T_sym : numpy.ndarray of shape (...,3,3) Symmetric tensor of which the eigenvalues are computed. + Returns + ------- + lambda : numpy.ndarray of shape (...,3) + Eigenvalues of T_sym sorted in ascending order, each repeated + according to its multiplicity. + """ return _np.linalg.eigvalsh(symmetric(T_sym)) def eigenvectors(T_sym,RHS=False): """ - Return eigenvectors of a symmetric tensor. - - The eigenvalues are sorted in ascending order of their associated eigenvalues. + Eigenvectors of a symmetric tensor. Parameters ---------- T_sym : numpy.ndarray of shape (...,3,3) Symmetric tensor of which the eigenvectors are computed. RHS: bool, optional - Enforce right-handed coordinate system. Default is False. + Enforce right-handed coordinate system. Defaults to False. + + Returns + ------- + x : numpy.ndarray of shape (...,3,3) + Eigenvectors of T_sym sorted in ascending order of their + associated eigenvalues. """ (u,v) = _np.linalg.eigh(symmetric(T_sym)) diff --git a/python/tests/test_Result.py b/python/tests/test_Result.py index b739bb2ca..89379ce53 100644 --- a/python/tests/test_Result.py +++ b/python/tests/test_Result.py @@ -292,7 +292,7 @@ class TestResult: default.add_stretch_tensor('F','U') loc = {'F': default.get_dataset_location('F'), 'U(F)': default.get_dataset_location('U(F)')} - in_memory = mechanics.right_stretch(default.read_dataset(loc['F'],0)) + in_memory = mechanics.stretch_right(default.read_dataset(loc['F'],0)) in_file = default.read_dataset(loc['U(F)'],0) assert np.allclose(in_memory,in_file) @@ -300,7 +300,7 @@ class TestResult: default.add_stretch_tensor('F','V') loc = {'F': default.get_dataset_location('F'), 'V(F)': default.get_dataset_location('V(F)')} - in_memory = mechanics.left_stretch(default.read_dataset(loc['F'],0)) + in_memory = mechanics.stretch_left(default.read_dataset(loc['F'],0)) in_file = default.read_dataset(loc['V(F)'],0) assert np.allclose(in_memory,in_file) diff --git a/python/tests/test_mechanics.py b/python/tests/test_mechanics.py index 4c0d68572..8af8d1a03 100644 --- a/python/tests/test_mechanics.py +++ b/python/tests/test_mechanics.py @@ -17,10 +17,6 @@ def eigenvalues(T_sym): return np.linalg.eigvalsh(symmetric(T_sym)) -def left_stretch(T): - return polar_decomposition(T,'V')[0] - - def maximum_shear(T_sym): w = eigenvalues(T_sym) return (w[0] - w[2])*0.5 @@ -39,10 +35,6 @@ def PK2(P,F): return symmetric(S) -def right_stretch(T): - return polar_decomposition(T,'U')[0] - - def rotational_part(T): return polar_decomposition(T,'R')[0] @@ -73,6 +65,14 @@ def strain_tensor(F,t,m): return eps +def stretch_left(T): + return polar_decomposition(T,'V')[0] + + +def stretch_right(T): + return polar_decomposition(T,'U')[0] + + def symmetric(T): return (T+T.T)*0.5 @@ -113,13 +113,13 @@ class TestMechanics: assert np.allclose(single(test_data_flat[i]),v) @pytest.mark.parametrize('vectorized,single',[(mechanics.deviatoric_part, deviatoric_part), - (mechanics.left_stretch , left_stretch ), (mechanics.maximum_shear , maximum_shear ), (mechanics.Mises_strain , Mises_strain ), (mechanics.Mises_stress , Mises_stress ), - (mechanics.right_stretch , right_stretch ), (mechanics.rotational_part, rotational_part), (mechanics.spherical_part , spherical_part ), + (mechanics.stretch_left , stretch_left ), + (mechanics.stretch_right , stretch_right ), ]) def test_vectorize_1_arg(self,vectorized,single): epsilon = np.random.rand(self.n,3,3) @@ -166,8 +166,8 @@ class TestMechanics: """F = RU = VR.""" F = np.broadcast_to(np.eye(3),[self.n,3,3])*np.random.rand(self.n,3,3) R = mechanics.rotational_part(F) - V = mechanics.left_stretch(F) - U = mechanics.right_stretch(F) + V = mechanics.stretch_left(F) + U = mechanics.stretch_right(F) assert np.allclose(np.matmul(R,U), np.matmul(V,R)) From 9b9d83d93c5b336aa95331b1c56774552bf00e5c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 16 Nov 2020 01:12:23 +0100 Subject: [PATCH 18/35] adding '_tensor' not needed --- python/damask/_result.py | 12 ++++++------ python/damask/mechanics.py | 11 ++++++++--- python/damask/seeds.py | 5 +---- python/tests/test_Result.py | 6 +++--- python/tests/test_mechanics.py | 22 +++++++++++----------- 5 files changed, 29 insertions(+), 27 deletions(-) diff --git a/python/damask/_result.py b/python/damask/_result.py index ad8c5ec0e..ca7d8ec71 100644 --- a/python/damask/_result.py +++ b/python/damask/_result.py @@ -991,21 +991,21 @@ class Result: @staticmethod - def _add_strain_tensor(F,t,m): + def _add_strain(F,t,m): return { - 'data': mechanics.strain_tensor(F['data'],t,m), + 'data': mechanics.strain(F['data'],t,m), 'label': f"epsilon_{t}^{m}({F['label']})", 'meta': { 'Unit': F['meta']['Unit'], 'Description': f"Strain tensor of {F['label']} ({F['meta']['Description']})", - 'Creator': 'add_strain_tensor' + 'Creator': 'add_strain' } } - def add_strain_tensor(self,F='F',t='V',m=0.0): + def add_strain(self,F='F',t='V',m=0.0): """ Add strain tensor of a deformation gradient. - For details refer to damask.mechanics.strain_tensor + For details refer to damask.mechanics.strain Parameters ---------- @@ -1018,7 +1018,7 @@ class Result: Order of the strain calculation. Defaults to ‘0.0’. """ - self._add_generic_pointwise(self._add_strain_tensor,{'F':F},{'t':t,'m':m}) + self._add_generic_pointwise(self._add_strain,{'F':F},{'t':t,'m':m}) @staticmethod diff --git a/python/damask/mechanics.py b/python/damask/mechanics.py index 9118ea875..a83f27370 100644 --- a/python/damask/mechanics.py +++ b/python/damask/mechanics.py @@ -40,9 +40,10 @@ def Cauchy_Green_deformation_right(F): """ return _np.matmul(tensor.transpose(F),F) + def Cauchy(P,F): """ - Calculate the Cauchy (true) stress. + Calculate the Cauchy stress (true stress). Resulting tensor is symmetrized as the Cauchy stress needs to be symmetric. @@ -201,9 +202,9 @@ def spherical_part(T,tensor=False): return _np.einsum('...jk,...->...jk',_np.eye(3),sph) if tensor else sph -def strain_tensor(F,t,m): +def strain(F,t,m): """ - Calculate strain tensor from deformation gradient. + Calculate strain tensor (Seth–Hill family). For details refer to https://en.wikipedia.org/wiki/Finite_strain_theory and https://de.wikipedia.org/wiki/Verzerrungstensor @@ -319,3 +320,7 @@ def _Mises(T_sym,s): """ d = deviatoric_part(T_sym) return _np.sqrt(s*_np.einsum('...jk->...',d**2.0)) + + +# for compatibility +strain_tensor = strain diff --git a/python/damask/seeds.py b/python/damask/seeds.py index c75dfb5e5..3c72952f7 100644 --- a/python/damask/seeds.py +++ b/python/damask/seeds.py @@ -1,7 +1,4 @@ -""" -Functionality for generation of seed points for Voronoi -or Laguerre tessellation. -""" +"""Functionality for generation of seed points for Voronoi or Laguerre tessellation.""" from scipy import spatial as _spatial import numpy as _np diff --git a/python/tests/test_Result.py b/python/tests/test_Result.py index 89379ce53..0613716f6 100644 --- a/python/tests/test_Result.py +++ b/python/tests/test_Result.py @@ -194,7 +194,7 @@ class TestResult: def test_add_Mises_strain(self,default): t = ['V','U'][np.random.randint(0,2)] m = np.random.random()*2.0 - 1.0 - default.add_strain_tensor('F',t,m) + default.add_strain('F',t,m) label = f'epsilon_{t}^{m}(F)' default.add_Mises(label) loc = {label :default.get_dataset_location(label), @@ -280,11 +280,11 @@ class TestResult: def test_add_strain(self,default): t = ['V','U'][np.random.randint(0,2)] m = np.random.random()*2.0 - 1.0 - default.add_strain_tensor('F',t,m) + default.add_strain('F',t,m) label = f'epsilon_{t}^{m}(F)' loc = {'F': default.get_dataset_location('F'), label: default.get_dataset_location(label)} - in_memory = mechanics.strain_tensor(default.read_dataset(loc['F'],0),t,m) + in_memory = mechanics.strain(default.read_dataset(loc['F'],0),t,m) in_file = default.read_dataset(loc[label],0) assert np.allclose(in_memory,in_file) diff --git a/python/tests/test_mechanics.py b/python/tests/test_mechanics.py index 8af8d1a03..78b441e18 100644 --- a/python/tests/test_mechanics.py +++ b/python/tests/test_mechanics.py @@ -44,7 +44,7 @@ def spherical_part(T,tensor=False): return sph if not tensor else np.eye(3)*sph -def strain_tensor(F,t,m): +def strain(F,t,m): if t == 'V': B = np.matmul(F,F.T) @@ -140,8 +140,8 @@ class TestMechanics: assert np.allclose(single(P[i],F[i]),v) - @pytest.mark.parametrize('vectorized,single',[(mechanics.strain_tensor,strain_tensor)]) - def test_vectorize_strain_tensor(self,vectorized,single): + @pytest.mark.parametrize('vectorized,single',[(mechanics.strain,strain)]) + def test_vectorize_strain(self,vectorized,single): F = np.random.rand(self.n,3,3) F_vec = np.reshape(F,(self.n//10,10,3,3)) t = ['V','U'][np.random.randint(0,2)] @@ -172,25 +172,25 @@ class TestMechanics: np.matmul(V,R)) @pytest.mark.parametrize('m',[0.0,np.random.random()*10.,np.random.random()*-10.]) - def test_strain_tensor_no_rotation(self,m): + def test_strain_no_rotation(self,m): """Ensure that left and right stretch give same results for no rotation.""" F = np.broadcast_to(np.eye(3),[self.n,3,3])*np.random.rand(self.n,3,3) - assert np.allclose(mechanics.strain_tensor(F,'U',m), - mechanics.strain_tensor(F,'V',m)) + assert np.allclose(mechanics.strain(F,'U',m), + mechanics.strain(F,'V',m)) @pytest.mark.parametrize('m',[0.0,np.random.random()*2.5,np.random.random()*-2.5]) - def test_strain_tensor_rotation_equivalence(self,m): + def test_strain_rotation_equivalence(self,m): """Ensure that left and right strain differ only by a rotation.""" F = np.broadcast_to(np.eye(3),[self.n,3,3]) + (np.random.rand(self.n,3,3)*0.5 - 0.25) - assert np.allclose(np.linalg.det(mechanics.strain_tensor(F,'U',m)), - np.linalg.det(mechanics.strain_tensor(F,'V',m))) + assert np.allclose(np.linalg.det(mechanics.strain(F,'U',m)), + np.linalg.det(mechanics.strain(F,'V',m))) @pytest.mark.parametrize('m',[0.0,np.random.random(),np.random.random()*-1.]) @pytest.mark.parametrize('t',['V','U']) - def test_strain_tensor_rotation(self,m,t): + def test_strain_rotation(self,m,t): """Ensure that pure rotation results in no strain.""" F = mechanics.rotational_part(np.random.rand(self.n,3,3)) - assert np.allclose(mechanics.strain_tensor(F,t,m), + assert np.allclose(mechanics.strain(F,t,m), 0.0) def test_rotation_determinant(self): From 5ebde607a29d669ef787b7f0ebbb08fa4991c178 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 16 Nov 2020 07:12:37 +0100 Subject: [PATCH 19/35] einsum simplifications --- python/damask/_orientation.py | 14 +++++++------- python/damask/_rotation.py | 2 +- python/damask/mechanics.py | 26 ++++++++++++-------------- python/tests/test_mechanics.py | 7 +++---- 4 files changed, 23 insertions(+), 26 deletions(-) diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index 9ae1d2ed3..ee2462287 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -1103,7 +1103,7 @@ class Orientation(Rotation): (np.array(hkil),np.array([[1,0,0,0], [0,1,0,0], [0,0,0,1]])) - return np.einsum('il,...l->...i',basis,axis) + return np.einsum('il,...l',basis,axis) @classmethod @@ -1133,7 +1133,7 @@ class Orientation(Rotation): [ 0, 1, 0], [-1,-1, 0], [ 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): @@ -1157,7 +1157,7 @@ class Orientation(Rotation): axis,basis = (np.array(direction),self.basis_reciprocal.T) \ if plane is None else \ (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): @@ -1183,9 +1183,9 @@ class Orientation(Rotation): if hkl is None else \ (np.array(hkl),self.basis_reciprocal) 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 - np.einsum('il,...l->...i',basis,axis)) + np.einsum('il,...l',basis,axis)) def to_pole(self,uvw=None,hkl=None,with_symmetry=False): @@ -1227,8 +1227,8 @@ class Orientation(Rotation): """ 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 = np.einsum('...i,...j->...ij',d/np.linalg.norm(d,axis=-1,keepdims=True), - p/np.linalg.norm(p,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)) return ~self.broadcast_to( self.shape+P.shape[:-2],mode='right') \ @ np.broadcast_to(P,self.shape+P.shape) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index e60b14ecb..4cf8cfcca 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -553,7 +553,7 @@ class Rotation: orthonormal = False # contains stretch if not orthonormal: (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)): raise ValueError('Orientation matrix has determinant ≠ 1.') if not np.all(np.isclose(np.einsum('...i,...i',om[...,0],om[...,1]), 0.0)) \ diff --git a/python/damask/mechanics.py b/python/damask/mechanics.py index a83f27370..ef2434e00 100644 --- a/python/damask/mechanics.py +++ b/python/damask/mechanics.py @@ -60,7 +60,7 @@ def Cauchy(P,F): 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) @@ -79,7 +79,7 @@ def deviatoric_part(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): @@ -157,7 +157,7 @@ def PK2(P,F): 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) @@ -199,7 +199,7 @@ def spherical_part(T,tensor=False): """ 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): @@ -231,14 +231,12 @@ def strain(F,t,m): w,n = _np.linalg.eigh(Cauchy_Green_deformation_right(F)) if m > 0.0: - eps = 1.0/(2.0*abs(m)) * (+ _np.matmul(n,_np.einsum('...j,...kj->...jk',w**m,n)) - - _np.einsum('...jk->...jk',_np.eye(3))) + eps = 1.0/(2.0*abs(m)) * (+ _np.einsum('...j,...kj,...lj',w**m,n,n) - _np.eye(3)) elif m < 0.0: - eps = 1.0/(2.0*abs(m)) * (- _np.matmul(n,_np.einsum('...j,...kj->...jk',w**m,n)) - + _np.einsum('...jk->...jk',_np.eye(3))) + eps = 1.0/(2.0*abs(m)) * (- _np.einsum('...j,...kj,...lj',w**m,n,n) + _np.eye(3)) 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 @@ -293,22 +291,22 @@ def _polar_decomposition(T,requested): """ u, _, vh = _np.linalg.svd(T) - R = _np.einsum('...ij,...jk->...ik',u,vh) + R = _np.einsum('...ij,...jk',u,vh) output = [] if 'R' in requested: output.append(R) 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: - output.append(_np.einsum('...ji,...jk->...ik',R,T)) + output.append(_np.einsum('...ji,...jk',R,T)) return tuple(output) 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 ---------- @@ -319,7 +317,7 @@ def _Mises(T_sym,s): """ 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 diff --git a/python/tests/test_mechanics.py b/python/tests/test_mechanics.py index 78b441e18..0b4393ed3 100644 --- a/python/tests/test_mechanics.py +++ b/python/tests/test_mechanics.py @@ -127,9 +127,8 @@ class TestMechanics: for i,v in enumerate(np.reshape(vectorized(epsilon_vec),vectorized(epsilon).shape)): assert np.allclose(single(epsilon[i]),v) - @pytest.mark.parametrize('vectorized,single',[ - (mechanics.Cauchy,Cauchy), - (mechanics.PK2 ,PK2 ) + @pytest.mark.parametrize('vectorized,single',[(mechanics.Cauchy,Cauchy), + (mechanics.PK2 ,PK2 ) ]) def test_vectorize_2_arg(self,vectorized,single): P = np.random.rand(self.n,3,3) @@ -227,7 +226,7 @@ class TestMechanics: scalar**3.0) 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) sph = mechanics.spherical_part(x,True) assert np.allclose(mechanics.Mises_strain(sph), From 6bedd84759b51e79525e6673abcf1a84e6ac284b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 17 Nov 2020 22:56:22 +0100 Subject: [PATCH 20/35] self-explanatory names --- python/damask/_geom.py | 2 +- python/damask/_orientation.py | 22 +++++++------- python/damask/_result.py | 32 ++++++++++++--------- python/damask/_rotation.py | 14 ++++----- python/damask/mechanics.py | 22 +++++++------- python/tests/test_Geom.py | 2 +- python/tests/test_Orientation.py | 22 +++++++------- python/tests/test_Result.py | 44 ++++++++++++++-------------- python/tests/test_Rotation.py | 42 +++++++++++++-------------- python/tests/test_mechanics.py | 49 ++++++++++++++++---------------- 10 files changed, 127 insertions(+), 124 deletions(-) diff --git a/python/damask/_geom.py b/python/damask/_geom.py index 3d73861ad..5eca451c6 100644 --- a/python/damask/_geom.py +++ b/python/damask/_geom.py @@ -753,7 +753,7 @@ class Geom: if fill is None: fill = np.nanmax(self.material) + 1 dtype = float if np.isnan(fill) or int(fill) != fill or self.material.dtype==np.float else int - Eulers = R.as_Eulers(degrees=True) + Eulers = R.as_Euler_angles(degrees=True) material_in = self.material.copy() # These rotations are always applied in the reference coordinate system, i.e. (z,x,z) not (z,x',z'') diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index ee2462287..a840a46db 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -86,8 +86,8 @@ class Orientation(Rotation): >>> damask.Orientation.from_random(shape=(3,5),lattice='tetragonal').reduced Disorientation between two specific orientations of hexagonal symmetry: - >>> a = damask.Orientation.from_Eulers(phi=[123,32,21],degrees=True,lattice='hexagonal') - >>> b = damask.Orientation.from_Eulers(phi=[104,11,87],degrees=True,lattice='hexagonal') + >>> a = damask.Orientation.from_Euler_angles(phi=[123,32,21],degrees=True,lattice='hexagonal') + >>> b = damask.Orientation.from_Euler_angles(phi=[104,11,87],degrees=True,lattice='hexagonal') >>> a.disorientation(b) Inverse pole figure color of the e_3 direction for a crystal in "Cube" orientation with cubic symmetry: @@ -95,7 +95,7 @@ class Orientation(Rotation): >>> o.IPF_color(o.to_SST(np.array([0,0,1]))) Schmid matrix (in lab frame) of slip systems of a face-centered cubic crystal in "Goss" orientation: - >>> damask.Orientation.from_Eulers(phi=[0,45,0],degrees=True,lattice='cF').Schmid('slip') + >>> damask.Orientation.from_Euler_angles(phi=[0,45,0],degrees=True,lattice='cF').Schmid('slip') """ @@ -291,9 +291,9 @@ class Orientation(Rotation): @classmethod - @extended_docstring(Rotation.from_Eulers) - def from_Eulers(cls,**kwargs): - return cls(rotation=Rotation.from_Eulers(**kwargs),**kwargs) + @extended_docstring(Rotation.from_Euler_angles) + def from_Euler_angles(cls,**kwargs): + return cls(rotation=Rotation.from_Euler_angles(**kwargs),**kwargs) @classmethod @@ -315,9 +315,9 @@ class Orientation(Rotation): @classmethod - @extended_docstring(Rotation.from_Rodrigues) - def from_Rodrigues(cls,**kwargs): - return cls(rotation=Rotation.from_Rodrigues(**kwargs),**kwargs) + @extended_docstring(Rotation.from_Rodrigues_vector) + def from_Rodrigues_vector(cls,**kwargs): + return cls(rotation=Rotation.from_Rodrigues_vector(**kwargs),**kwargs) @classmethod @@ -496,7 +496,7 @@ class Orientation(Rotation): if self.family is None: raise ValueError('Missing crystal symmetry') - rho_abs = np.abs(self.as_Rodrigues(vector=True)) + rho_abs = np.abs(self.as_Rodrigues_vector(vector=True)) with np.errstate(invalid='ignore'): # using '*'/prod for 'and' @@ -539,7 +539,7 @@ class Orientation(Rotation): if self.family is None: raise ValueError('Missing crystal symmetry') - rho = self.as_Rodrigues(vector=True) + rho = self.as_Rodrigues_vector(vector=True) with np.errstate(invalid='ignore'): if self.family == 'cubic': diff --git a/python/damask/_result.py b/python/damask/_result.py index ca7d8ec71..9859f8c4d 100644 --- a/python/damask/_result.py +++ b/python/damask/_result.py @@ -597,19 +597,19 @@ class Result: @staticmethod - def _add_Cauchy(P,F): + def _add_stress_Cauchy(P,F): return { - 'data': mechanics.Cauchy(P['data'],F['data']), + 'data': mechanics.stress_Cauchy(P['data'],F['data']), 'label': 'sigma', 'meta': { 'Unit': P['meta']['Unit'], 'Description': "Cauchy stress calculated " f"from {P['label']} ({P['meta']['Description']})" f" and {F['label']} ({F['meta']['Description']})", - 'Creator': 'add_Cauchy' + 'Creator': 'add_stress_Cauchy' } } - def add_Cauchy(self,P='P',F='F'): + def add_stress_Cauchy(self,P='P',F='F'): """ Add Cauchy stress calculated from first Piola-Kirchhoff stress and deformation gradient. @@ -621,7 +621,7 @@ class Result: Label of the dataset containing the deformation gradient. Defaults to ‘F’. """ - self._add_generic_pointwise(self._add_Cauchy,{'P':P,'F':F}) + self._add_generic_pointwise(self._add_stress_Cauchy,{'P':P,'F':F}) @staticmethod @@ -798,7 +798,7 @@ class Result: @staticmethod - def _add_Mises(T_sym,kind): + def _add_equivalent_Mises(T_sym,kind): k = kind if k is None: if T_sym['meta']['Unit'] == '1': @@ -809,7 +809,8 @@ class Result: raise ValueError('invalid von Mises kind {kind}') return { - 'data': (mechanics.Mises_strain if k=='strain' else mechanics.Mises_stress)(T_sym['data']), + 'data': (mechanics.equivalent_strain_Mises if k=='strain' else \ + mechanics.equivalent_stress_Mises)(T_sym['data']), 'label': f"{T_sym['label']}_vM", 'meta': { 'Unit': T_sym['meta']['Unit'], @@ -817,7 +818,7 @@ class Result: 'Creator': 'add_Mises' } } - def add_Mises(self,T_sym,kind=None): + def add_equivalent_Mises(self,T_sym,kind=None): """ Add the equivalent Mises stress or strain of a symmetric tensor. @@ -830,7 +831,7 @@ class Result: it is selected based on the unit of the dataset ('1' -> strain, 'Pa' -> stress'). """ - self._add_generic_pointwise(self._add_Mises,{'T_sym':T_sym},{'kind':kind}) + self._add_generic_pointwise(self._add_equivalent_Mises,{'T_sym':T_sym},{'kind':kind}) @staticmethod @@ -872,19 +873,19 @@ class Result: @staticmethod - def _add_PK2(P,F): + def _add_stress_second_Piola_Kirchhoff(P,F): return { - 'data': mechanics.PK2(P['data'],F['data']), + 'data': mechanics.stress_second_Piola_Kirchhoff(P['data'],F['data']), 'label': 'S', 'meta': { 'Unit': P['meta']['Unit'], 'Description': "2. Piola-Kirchhoff stress calculated " f"from {P['label']} ({P['meta']['Description']})" f" and {F['label']} ({F['meta']['Description']})", - 'Creator': 'add_PK2' + 'Creator': 'add_stress_second_Piola_Kirchhoff' } } - def add_PK2(self,P='P',F='F'): + def add_stress_second_Piola_Kirchhoff(self,P='P',F='F'): """ Add second Piola-Kirchhoff stress calculated from first Piola-Kirchhoff stress and deformation gradient. @@ -896,7 +897,7 @@ class Result: Label of deformation gradient dataset. Defaults to ‘F’. """ - self._add_generic_pointwise(self._add_PK2,{'P':P,'F':F}) + self._add_generic_pointwise(self._add_stress_second_Piola_Kirchhoff,{'P':P,'F':F}) # The add_pole functionality needs discussion. @@ -1297,3 +1298,6 @@ class Result: v.add(u,'u') v.save(f'{self.fname.stem}_inc{inc[3:].zfill(N_digits)}') + + +Result.add_PK2 = Result.add_stress_second_Piola_Kirchhoff diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 4cf8cfcca..e695ce407 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -66,12 +66,12 @@ class Rotation: def __repr__(self): """Represent rotation as unit quaternion, rotation matrix, and Bunge-Euler angles.""" - return 'Quaternions:\n'+str(self.quaternion) \ + return 'As quaternions:\n'+str(self.quaternion) \ if self.quaternion.shape != (4,) else \ '\n'.join([ 'Quaternion: (real={:.3f}, imag=<{:+.3f}, {:+.3f}, {:+.3f}>)'.format(*(self.quaternion)), 'Matrix:\n{}'.format(np.round(self.as_matrix(),8)), - 'Bunge Eulers / deg: ({:3.2f}, {:3.2f}, {:3.2f})'.format(*self.as_Eulers(degrees=True)), + 'Bunge Eulers / deg: ({:3.2f}, {:3.2f}, {:3.2f})'.format(*self.as_Euler_angles(degrees=True)), ]) @@ -314,7 +314,7 @@ class Rotation: """ return self.quaternion.copy() - def as_Eulers(self, + def as_Euler_angles(self, degrees = False): """ Represent as Bunge-Euler angles. @@ -372,7 +372,7 @@ class Rotation: """ return Rotation._qu2om(self.quaternion) - def as_Rodrigues(self, + def as_Rodrigues_vector(self, vector = False): """ Represent as Rodrigues-Frank vector with separated axis and angle argument. @@ -462,7 +462,7 @@ class Rotation: return Rotation(qu) @staticmethod - def from_Eulers(phi, + def from_Euler_angles(phi, degrees = False, **kwargs): """ @@ -606,7 +606,7 @@ class Rotation: @staticmethod - def from_Rodrigues(rho, + def from_Rodrigues_vector(rho, normalize = False, P = -1, **kwargs): @@ -784,7 +784,7 @@ class Rotation: dg = 1.0 if fractions else _dg(Eulers,degrees) dV_V = dg * np.maximum(0.0,weights.squeeze()) - return Rotation.from_Eulers(Eulers[util.hybrid_IA(dV_V,N,seed)],degrees) + return Rotation.from_Euler_angles(Eulers[util.hybrid_IA(dV_V,N,seed)],degrees) @staticmethod diff --git a/python/damask/mechanics.py b/python/damask/mechanics.py index ef2434e00..bab41a18c 100644 --- a/python/damask/mechanics.py +++ b/python/damask/mechanics.py @@ -5,7 +5,7 @@ from . import tensor import numpy as _np -def Cauchy_Green_deformation_left(F): +def deformation_Cauchy_Green_left(F): """ Calculate left Cauchy-Green deformation tensor (Finger deformation tensor). @@ -23,7 +23,7 @@ def Cauchy_Green_deformation_left(F): return _np.matmul(F,tensor.transpose(F)) -def Cauchy_Green_deformation_right(F): +def deformation_Cauchy_Green_right(F): """ Calculate right Cauchy-Green deformation tensor. @@ -41,7 +41,7 @@ def Cauchy_Green_deformation_right(F): return _np.matmul(tensor.transpose(F),F) -def Cauchy(P,F): +def stress_Cauchy(P,F): """ Calculate the Cauchy stress (true stress). @@ -101,7 +101,7 @@ def maximum_shear(T_sym): return (w[...,0] - w[...,2])*0.5 -def Mises_strain(epsilon): +def equivalent_strain_Mises(epsilon): """ Calculate the Mises equivalent of a strain tensor. @@ -116,10 +116,10 @@ def Mises_strain(epsilon): Von Mises equivalent strain of epsilon. """ - return _Mises(epsilon,2.0/3.0) + return _equivalent_Mises(epsilon,2.0/3.0) -def Mises_stress(sigma): +def equivalent_stress_Mises(sigma): """ Calculate the Mises equivalent of a stress tensor. @@ -134,10 +134,10 @@ def Mises_stress(sigma): Von Mises equivalent stress of sigma. """ - return _Mises(sigma,3.0/2.0) + return _equivalent_Mises(sigma,3.0/2.0) -def PK2(P,F): +def stress_second_Piola_Kirchhoff(P,F): """ Calculate the second Piola-Kirchhoff stress. @@ -226,9 +226,9 @@ def strain(F,t,m): """ if t == 'V': - w,n = _np.linalg.eigh(Cauchy_Green_deformation_left(F)) + w,n = _np.linalg.eigh(deformation_Cauchy_Green_left(F)) elif t == 'U': - w,n = _np.linalg.eigh(Cauchy_Green_deformation_right(F)) + w,n = _np.linalg.eigh(deformation_Cauchy_Green_right(F)) if m > 0.0: eps = 1.0/(2.0*abs(m)) * (+ _np.einsum('...j,...kj,...lj',w**m,n,n) - _np.eye(3)) @@ -304,7 +304,7 @@ def _polar_decomposition(T,requested): return tuple(output) -def _Mises(T_sym,s): +def _equivalent_Mises(T_sym,s): """ Base equation for Mises equivalent of a stress or strain tensor. diff --git a/python/tests/test_Geom.py b/python/tests/test_Geom.py index a9b9c91b5..b9b8d33b2 100644 --- a/python/tests/test_Geom.py +++ b/python/tests/test_Geom.py @@ -223,7 +223,7 @@ class TestGeom: @pytest.mark.parametrize('Eulers',[[32.0,68.0,21.0], [0.0,32.0,240.0]]) def test_rotate(self,default,update,reference_dir,Eulers): - modified = default.rotate(Rotation.from_Eulers(Eulers,degrees=True)) + modified = default.rotate(Rotation.from_Euler_angles(Eulers,degrees=True)) tag = f'Eulers_{util.srepr(Eulers,"-")}' reference = reference_dir/f'rotate_{tag}.vtr' if update: modified.save(reference) diff --git a/python/tests/test_Orientation.py b/python/tests/test_Orientation.py index fe35bda14..8a3f85e78 100644 --- a/python/tests/test_Orientation.py +++ b/python/tests/test_Orientation.py @@ -16,7 +16,7 @@ def reference_dir(reference_dir_base): @pytest.fixture def set_of_rodrigues(set_of_quaternions): - return Rotation(set_of_quaternions).as_Rodrigues()[:200] + return Rotation(set_of_quaternions).as_Rodrigues_vector()[:200] class TestOrientation: @@ -90,8 +90,8 @@ class TestOrientation: assert np.all(Orientation.from_quaternion(q=np.array([1,0,0,0]),lattice='triclinic').as_matrix() == np.eye(3)) - def test_from_Eulers(self): - assert np.all(Orientation.from_Eulers(phi=np.zeros(3),lattice='triclinic').as_matrix() + def test_from_Euler_angles(self): + assert np.all(Orientation.from_Euler_angles(phi=np.zeros(3),lattice='triclinic').as_matrix() == np.eye(3)) def test_from_axis_angle(self): @@ -106,8 +106,8 @@ class TestOrientation: assert np.all(Orientation.from_matrix(R=np.eye(3),lattice='triclinic').as_matrix() == np.eye(3)) - def test_from_Rodrigues(self): - assert np.all(Orientation.from_Rodrigues(rho=np.array([0,0,1,0]),lattice='triclinic').as_matrix() + def test_from_Rodrigues_vector(self): + assert np.all(Orientation.from_Rodrigues_vector(rho=np.array([0,0,1,0]),lattice='triclinic').as_matrix() == np.eye(3)) def test_from_homochoric(self): @@ -208,7 +208,7 @@ class TestOrientation: @pytest.mark.parametrize('lattice',Orientation.crystal_families) def test_disorientation360(self,lattice): o_1 = Orientation(Rotation(),lattice) - o_2 = Orientation.from_Eulers(lattice=lattice,phi=[360,0,0],degrees=True) + o_2 = Orientation.from_Euler_angles(lattice=lattice,phi=[360,0,0],degrees=True) assert np.allclose((o_1.disorientation(o_2)).as_matrix(),np.eye(3)) @pytest.mark.parametrize('lattice',Orientation.crystal_families) @@ -275,16 +275,16 @@ class TestOrientation: @pytest.mark.parametrize('lattice',Orientation.crystal_families) def test_in_FZ_vectorization(self,set_of_rodrigues,lattice): - result = Orientation.from_Rodrigues(rho=set_of_rodrigues.reshape((50,4,-1)),lattice=lattice).in_FZ.reshape(-1) + result = Orientation.from_Rodrigues_vector(rho=set_of_rodrigues.reshape((50,4,-1)),lattice=lattice).in_FZ.reshape(-1) for r,rho in zip(result,set_of_rodrigues[:len(result)]): - assert r == Orientation.from_Rodrigues(rho=rho,lattice=lattice).in_FZ + assert r == Orientation.from_Rodrigues_vector(rho=rho,lattice=lattice).in_FZ @pytest.mark.parametrize('lattice',Orientation.crystal_families) def test_in_disorientation_FZ_vectorization(self,set_of_rodrigues,lattice): - result = Orientation.from_Rodrigues(rho=set_of_rodrigues.reshape((50,4,-1)), + result = Orientation.from_Rodrigues_vector(rho=set_of_rodrigues.reshape((50,4,-1)), lattice=lattice).in_disorientation_FZ.reshape(-1) for r,rho in zip(result,set_of_rodrigues[:len(result)]): - assert r == Orientation.from_Rodrigues(rho=rho,lattice=lattice).in_disorientation_FZ + assert r == Orientation.from_Rodrigues_vector(rho=rho,lattice=lattice).in_disorientation_FZ @pytest.mark.parametrize('proper',[True,False]) @pytest.mark.parametrize('lattice',Orientation.crystal_families) @@ -417,7 +417,7 @@ class TestOrientation: def test_relationship_reference(self,update,reference_dir,model,lattice): reference = reference_dir/f'{lattice}_{model}.txt' o = Orientation(lattice=lattice) - eu = o.related(model).as_Eulers(degrees=True) + eu = o.related(model).as_Euler_angles(degrees=True) if update: coords = np.array([(1,i+1) for i,x in enumerate(eu)]) Table(eu,{'Eulers':(3,)})\ diff --git a/python/tests/test_Result.py b/python/tests/test_Result.py index 0613716f6..46b947cdc 100644 --- a/python/tests/test_Result.py +++ b/python/tests/test_Result.py @@ -121,13 +121,13 @@ class TestResult: in_file = default.read_dataset(loc['x'],0) assert np.allclose(in_memory,in_file) - def test_add_Cauchy(self,default): - default.add_Cauchy('P','F') + def test_add_stress_Cauchy(self,default): + default.add_stress_Cauchy('P','F') loc = {'F': default.get_dataset_location('F'), 'P': default.get_dataset_location('P'), 'sigma':default.get_dataset_location('sigma')} - in_memory = mechanics.Cauchy(default.read_dataset(loc['P'],0), - default.read_dataset(loc['F'],0)) + in_memory = mechanics.stress_Cauchy(default.read_dataset(loc['P'],0), + default.read_dataset(loc['F'],0)) in_file = default.read_dataset(loc['sigma'],0) assert np.allclose(in_memory,in_file) @@ -149,7 +149,7 @@ class TestResult: @pytest.mark.parametrize('eigenvalue,function',[('max',np.amax),('min',np.amin)]) def test_add_eigenvalue(self,default,eigenvalue,function): - default.add_Cauchy('P','F') + default.add_stress_Cauchy('P','F') default.add_eigenvalue('sigma',eigenvalue) loc = {'sigma' :default.get_dataset_location('sigma'), 'lambda':default.get_dataset_location(f'lambda_{eigenvalue}(sigma)')} @@ -159,7 +159,7 @@ class TestResult: @pytest.mark.parametrize('eigenvalue,idx',[('max',2),('mid',1),('min',0)]) def test_add_eigenvector(self,default,eigenvalue,idx): - default.add_Cauchy('P','F') + default.add_stress_Cauchy('P','F') default.add_eigenvector('sigma',eigenvalue) loc = {'sigma' :default.get_dataset_location('sigma'), 'v(sigma)':default.get_dataset_location(f'v_{eigenvalue}(sigma)')} @@ -183,7 +183,7 @@ class TestResult: assert np.allclose(in_memory,in_file) def test_add_maximum_shear(self,default): - default.add_Cauchy('P','F') + default.add_stress_Cauchy('P','F') default.add_maximum_shear('sigma') loc = {'sigma' :default.get_dataset_location('sigma'), 'max_shear(sigma)':default.get_dataset_location('max_shear(sigma)')} @@ -196,34 +196,34 @@ class TestResult: m = np.random.random()*2.0 - 1.0 default.add_strain('F',t,m) label = f'epsilon_{t}^{m}(F)' - default.add_Mises(label) + default.add_equivalent_Mises(label) loc = {label :default.get_dataset_location(label), label+'_vM':default.get_dataset_location(label+'_vM')} - in_memory = mechanics.Mises_strain(default.read_dataset(loc[label],0)).reshape(-1,1) + in_memory = mechanics.equivalent_strain_Mises(default.read_dataset(loc[label],0)).reshape(-1,1) in_file = default.read_dataset(loc[label+'_vM'],0) assert np.allclose(in_memory,in_file) def test_add_Mises_stress(self,default): - default.add_Cauchy('P','F') - default.add_Mises('sigma') + default.add_stress_Cauchy('P','F') + default.add_equivalent_Mises('sigma') loc = {'sigma' :default.get_dataset_location('sigma'), 'sigma_vM':default.get_dataset_location('sigma_vM')} - in_memory = mechanics.Mises_stress(default.read_dataset(loc['sigma'],0)).reshape(-1,1) + in_memory = mechanics.equivalent_stress_Mises(default.read_dataset(loc['sigma'],0)).reshape(-1,1) in_file = default.read_dataset(loc['sigma_vM'],0) assert np.allclose(in_memory,in_file) def test_add_Mises_invalid(self,default): - default.add_Cauchy('P','F') + default.add_stress_Cauchy('P','F') default.add_calculation('sigma_y','#sigma#',unit='y') - default.add_Mises('sigma_y') + default.add_equivalent_Mises('sigma_y') assert default.get_dataset_location('sigma_y_vM') == [] def test_add_Mises_stress_strain(self,default): - default.add_Cauchy('P','F') + default.add_stress_Cauchy('P','F') default.add_calculation('sigma_y','#sigma#',unit='y') default.add_calculation('sigma_x','#sigma#',unit='x') - default.add_Mises('sigma_y',kind='strain') - default.add_Mises('sigma_x',kind='stress') + default.add_equivalent_Mises('sigma_y',kind='strain') + default.add_equivalent_Mises('sigma_x',kind='stress') loc = {'y' :default.get_dataset_location('sigma_y_vM'), 'x' :default.get_dataset_location('sigma_x_vM')} assert not np.allclose(default.read_dataset(loc['y'],0),default.read_dataset(loc['x'],0)) @@ -236,13 +236,13 @@ class TestResult: in_file = default.read_dataset(loc['|F|_1'],0) assert np.allclose(in_memory,in_file) - def test_add_PK2(self,default): - default.add_PK2('P','F') + def test_add_stress_second_Piola_Kirchhoff(self,default): + default.add_stress_second_Piola_Kirchhoff('P','F') loc = {'F':default.get_dataset_location('F'), 'P':default.get_dataset_location('P'), 'S':default.get_dataset_location('S')} - in_memory = mechanics.PK2(default.read_dataset(loc['P'],0), - default.read_dataset(loc['F'],0)) + in_memory = mechanics.stress_second_Piola_Kirchhoff(default.read_dataset(loc['P'],0), + default.read_dataset(loc['F'],0)) in_file = default.read_dataset(loc['S'],0) assert np.allclose(in_memory,in_file) @@ -312,7 +312,7 @@ class TestResult: def test_add_overwrite(self,default,overwrite): default.pick('times',default.times_in_range(0,np.inf)[-1]) - default.add_Cauchy() + default.add_stress_Cauchy() loc = default.get_dataset_location('sigma') with h5py.File(default.fname,'r') as f: # h5py3 compatibility diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index 9361008e5..ad4e4e702 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -522,7 +522,7 @@ class TestRotation: def test_Eulers_internal(self,set_of_rotations,forward,backward): """Ensure invariance of conversion from Euler angles and back.""" for rot in set_of_rotations: - m = rot.as_Eulers() + m = rot.as_Euler_angles() o = backward(forward(m)) u = np.array([np.pi*2,np.pi,np.pi*2]) ok = np.allclose(m,o,atol=atol) @@ -559,7 +559,7 @@ class TestRotation: """Ensure invariance of conversion from Rodrigues-Frank vector and back.""" cutoff = np.tan(np.pi*.5*(1.-1e-4)) for rot in set_of_rotations: - m = rot.as_Rodrigues() + m = rot.as_Rodrigues_vector() o = backward(forward(m)) ok = np.allclose(np.clip(m,None,cutoff),np.clip(o,None,cutoff),atol=atol) ok = ok or np.isclose(m[3],0.0,atol=atol) @@ -626,7 +626,7 @@ class TestRotation: (Rotation._eu2ro,eu2ro)]) def test_Eulers_vectorization(self,set_of_rotations,vectorized,single): """Check vectorized implementation for Euler angles against single point calculation.""" - eu = np.array([rot.as_Eulers() for rot in set_of_rotations]) + eu = np.array([rot.as_Euler_angles() for rot in set_of_rotations]) vectorized(eu.reshape(eu.shape[0]//2,-1,3)) co = vectorized(eu) for e,c in zip(eu,co): @@ -649,7 +649,7 @@ class TestRotation: (Rotation._ro2ho,ro2ho)]) def test_Rodrigues_vectorization(self,set_of_rotations,vectorized,single): """Check vectorized implementation for Rodrigues-Frank vector against single point calculation.""" - ro = np.array([rot.as_Rodrigues() for rot in set_of_rotations]) + ro = np.array([rot.as_Rodrigues_vector() for rot in set_of_rotations]) vectorized(ro.reshape(ro.shape[0]//2,-1,4)) co = vectorized(ro) for r,c in zip(ro,co): @@ -687,7 +687,7 @@ class TestRotation: def test_Eulers(self,set_of_rotations,degrees): for rot in set_of_rotations: m = rot.as_quaternion() - o = Rotation.from_Eulers(rot.as_Eulers(degrees),degrees).as_quaternion() + o = Rotation.from_Euler_angles(rot.as_Euler_angles(degrees),degrees).as_quaternion() ok = np.allclose(m,o,atol=atol) if np.isclose(rot.as_quaternion()[0],0.0,atol=atol): ok |= np.allclose(m*-1.,o,atol=atol) @@ -699,8 +699,8 @@ class TestRotation: def test_axis_angle(self,set_of_rotations,degrees,normalize,P): c = np.array([P*-1,P*-1,P*-1,1.]) for rot in set_of_rotations: - m = rot.as_Eulers() - o = Rotation.from_axis_angle(rot.as_axis_angle(degrees)*c,degrees,normalize,P).as_Eulers() + m = rot.as_Euler_angles() + o = Rotation.from_axis_angle(rot.as_axis_angle(degrees)*c,degrees,normalize,P).as_Euler_angles() u = np.array([np.pi*2,np.pi,np.pi*2]) ok = np.allclose(m,o,atol=atol) ok |= np.allclose(np.where(np.isclose(m,u),m-u,m),np.where(np.isclose(o,u),o-u,o),atol=atol) @@ -726,7 +726,7 @@ class TestRotation: c = np.array([P*-1,P*-1,P*-1,1.]) for rot in set_of_rotations: m = rot.as_matrix() - o = Rotation.from_Rodrigues(rot.as_Rodrigues()*c,normalize,P).as_matrix() + o = Rotation.from_Rodrigues_vector(rot.as_Rodrigues_vector()*c,normalize,P).as_matrix() ok = np.allclose(m,o,atol=atol) assert ok and np.isclose(np.linalg.det(o),1.0), f'{m},{o}' @@ -734,8 +734,8 @@ class TestRotation: def test_homochoric(self,set_of_rotations,P): cutoff = np.tan(np.pi*.5*(1.-1e-4)) for rot in set_of_rotations: - m = rot.as_Rodrigues() - o = Rotation.from_homochoric(rot.as_homochoric()*P*-1,P).as_Rodrigues() + m = rot.as_Rodrigues_vector() + o = Rotation.from_homochoric(rot.as_homochoric()*P*-1,P).as_Rodrigues_vector() ok = np.allclose(np.clip(m,None,cutoff),np.clip(o,None,cutoff),atol=atol) ok = ok or np.isclose(m[3],0.0,atol=atol) assert ok and np.isclose(np.linalg.norm(o[:3]),1.0), f'{m},{o},{rot.as_quaternion()}' @@ -819,10 +819,10 @@ class TestRotation: assert r == r.flatten(order).reshape(shape,order) @pytest.mark.parametrize('function',[Rotation.from_quaternion, - Rotation.from_Eulers, + Rotation.from_Euler_angles, Rotation.from_axis_angle, Rotation.from_matrix, - Rotation.from_Rodrigues, + Rotation.from_Rodrigues_vector, Rotation.from_homochoric, Rotation.from_cubochoric]) def test_invalid_shape(self,function): @@ -832,7 +832,7 @@ class TestRotation: @pytest.mark.parametrize('fr,to',[(Rotation.from_quaternion,'as_quaternion'), (Rotation.from_axis_angle,'as_axis_angle'), - (Rotation.from_Rodrigues, 'as_Rodrigues'), + (Rotation.from_Rodrigues_vector, 'as_Rodrigues_vector'), (Rotation.from_homochoric,'as_homochoric'), (Rotation.from_cubochoric,'as_cubochoric')]) def test_invalid_P(self,fr,to): @@ -852,13 +852,13 @@ class TestRotation: @pytest.mark.parametrize('function,invalid',[(Rotation.from_quaternion, np.array([-1,0,0,0])), (Rotation.from_quaternion, np.array([1,1,1,0])), - (Rotation.from_Eulers, np.array([1,4,0])), + (Rotation.from_Euler_angles, np.array([1,4,0])), (Rotation.from_axis_angle, np.array([1,0,0,4])), (Rotation.from_axis_angle, np.array([1,1,0,1])), (Rotation.from_matrix, np.random.rand(3,3)), (Rotation.from_matrix, np.array([[1,1,0],[1,2,0],[0,0,1]])), - (Rotation.from_Rodrigues, np.array([1,0,0,-1])), - (Rotation.from_Rodrigues, np.array([1,1,0,1])), + (Rotation.from_Rodrigues_vector, np.array([1,0,0,-1])), + (Rotation.from_Rodrigues_vector, np.array([1,1,0,1])), (Rotation.from_homochoric, np.array([2,2,2])), (Rotation.from_cubochoric, np.array([1.1,0,0])) ]) def test_invalid_value(self,function,invalid): @@ -904,8 +904,8 @@ class TestRotation: def test_rotate_360deg(self,data): phi_1 = np.random.random() * np.pi phi_2 = 2*np.pi - phi_1 - R_1 = Rotation.from_Eulers(np.array([phi_1,0.,0.])) - R_2 = Rotation.from_Eulers(np.array([0.,0.,phi_2])) + R_1 = Rotation.from_Euler_angles(np.array([phi_1,0.,0.])) + R_2 = Rotation.from_Euler_angles(np.array([0.,0.,phi_2])) assert np.allclose(data,R_2@(R_1@data)) @pytest.mark.parametrize('pwr',[-10,0,1,2.5,np.pi,np.random.random()]) @@ -951,7 +951,7 @@ class TestRotation: def test_misorientation360(self): R_1 = Rotation() - R_2 = Rotation.from_Eulers([360,0,0],degrees=True) + R_2 = Rotation.from_Euler_angles([360,0,0],degrees=True) assert np.allclose(R_1.misorientation(R_2).as_matrix(),np.eye(3)) @pytest.mark.parametrize('angle',[10,20,30,40,50,60,70,80,90,100,120]) @@ -1014,7 +1014,7 @@ class TestRotation: Eulers = grid_filters.cell_coord0(steps,limits) Eulers = np.radians(Eulers) if not degrees else Eulers - Eulers_r = Rotation.from_ODF(weights,Eulers.reshape(-1,3,order='F'),N,degrees,fractions).as_Eulers(True) + Eulers_r = Rotation.from_ODF(weights,Eulers.reshape(-1,3,order='F'),N,degrees,fractions).as_Euler_angles(True) weights_r = np.histogramdd(Eulers_r,steps,rng)[0].flatten(order='F')/N * np.sum(weights) if fractions: assert np.sqrt(((weights_r - weights) ** 2).mean()) < 4 @@ -1032,7 +1032,7 @@ class TestRotation: Eulers = grid_filters.node_coord0(steps,limits)[:-1,:-1,:-1] Eulers = np.radians(Eulers) if not degrees else Eulers - Eulers_r = Rotation.from_ODF(weights,Eulers.reshape(-1,3,order='F'),N,degrees).as_Eulers(True) + Eulers_r = Rotation.from_ODF(weights,Eulers.reshape(-1,3,order='F'),N,degrees).as_Euler_angles(True) weights_r = np.histogramdd(Eulers_r,steps,rng)[0].flatten(order='F')/N * np.sum(weights) assert np.sqrt(((weights_r - weights) ** 2).mean()) < 5 diff --git a/python/tests/test_mechanics.py b/python/tests/test_mechanics.py index 0b4393ed3..4a61f5b1e 100644 --- a/python/tests/test_mechanics.py +++ b/python/tests/test_mechanics.py @@ -4,7 +4,7 @@ import numpy as np from damask import tensor from damask import mechanics -def Cauchy(P,F): +def stress_Cauchy(P,F): sigma = 1.0/np.linalg.det(F) * np.dot(P,F.T) return symmetric(sigma) @@ -22,15 +22,15 @@ def maximum_shear(T_sym): return (w[0] - w[2])*0.5 -def Mises_strain(epsilon): - return Mises(epsilon,2.0/3.0) +def equivalent_strain_Mises(epsilon): + return equivalent_Mises(epsilon,2.0/3.0) -def Mises_stress(sigma): - return Mises(sigma,3.0/2.0) +def equivalent_stress_Mises(sigma): + return equivalent_Mises(sigma,3.0/2.0) -def PK2(P,F): +def stress_second_Piola_Kirchhoff(P,F): S = np.dot(np.linalg.inv(F),P) return symmetric(S) @@ -91,9 +91,8 @@ def polar_decomposition(T,requested): return tuple(output) -def Mises(T_sym,s): - d = deviatoric_part(T_sym) - return np.sqrt(s*(np.sum(d**2.0))) +def equivalent_Mises(T_sym,s): + return np.sqrt(s*(np.sum(deviatoric_part(T_sym)**2.0))) class TestMechanics: @@ -112,14 +111,14 @@ class TestMechanics: for i,v in enumerate(np.reshape(vectorized(test_data),vectorized(test_data_flat).shape)): assert np.allclose(single(test_data_flat[i]),v) - @pytest.mark.parametrize('vectorized,single',[(mechanics.deviatoric_part, deviatoric_part), - (mechanics.maximum_shear , maximum_shear ), - (mechanics.Mises_strain , Mises_strain ), - (mechanics.Mises_stress , Mises_stress ), - (mechanics.rotational_part, rotational_part), - (mechanics.spherical_part , spherical_part ), - (mechanics.stretch_left , stretch_left ), - (mechanics.stretch_right , stretch_right ), + @pytest.mark.parametrize('vectorized,single',[(mechanics.deviatoric_part, deviatoric_part), + (mechanics.maximum_shear, maximum_shear), + (mechanics.equivalent_stress_Mises, equivalent_stress_Mises), + (mechanics.equivalent_strain_Mises, equivalent_strain_Mises), + (mechanics.rotational_part, rotational_part), + (mechanics.spherical_part, spherical_part), + (mechanics.stretch_left, stretch_left), + (mechanics.stretch_right, stretch_right), ]) def test_vectorize_1_arg(self,vectorized,single): epsilon = np.random.rand(self.n,3,3) @@ -127,8 +126,8 @@ class TestMechanics: for i,v in enumerate(np.reshape(vectorized(epsilon_vec),vectorized(epsilon).shape)): assert np.allclose(single(epsilon[i]),v) - @pytest.mark.parametrize('vectorized,single',[(mechanics.Cauchy,Cauchy), - (mechanics.PK2 ,PK2 ) + @pytest.mark.parametrize('vectorized,single',[(mechanics.stress_Cauchy, stress_Cauchy), + (mechanics.stress_second_Piola_Kirchhoff, stress_second_Piola_Kirchhoff) ]) def test_vectorize_2_arg(self,vectorized,single): P = np.random.rand(self.n,3,3) @@ -148,8 +147,8 @@ class TestMechanics: for i,v in enumerate(np.reshape(vectorized(F_vec,t,m),vectorized(F,t,m).shape)): assert np.allclose(single(F[i],t,m),v) - @pytest.mark.parametrize('function',[mechanics.Cauchy, - mechanics.PK2, + @pytest.mark.parametrize('function',[mechanics.stress_Cauchy, + mechanics.stress_second_Piola_Kirchhoff, ]) def test_stress_measures(self,function): """Ensure that all stress measures are equivalent for no deformation.""" @@ -212,8 +211,8 @@ class TestMechanics: def test_deviatoric_Mises(self): """Ensure that Mises equivalent stress depends only on deviatoric part.""" x = np.random.rand(self.n,3,3) - full = mechanics.Mises_stress(x) - dev = mechanics.Mises_stress(mechanics.deviatoric_part(x)) + full = mechanics.equivalent_stress_Mises(x) + dev = mechanics.equivalent_stress_Mises(mechanics.deviatoric_part(x)) assert np.allclose(full, dev) @@ -229,14 +228,14 @@ class TestMechanics: """Ensure that Mises equivalent strain of spherical strain is 0.""" x = np.random.rand(self.n,3,3) sph = mechanics.spherical_part(x,True) - assert np.allclose(mechanics.Mises_strain(sph), + assert np.allclose(mechanics.equivalent_strain_Mises(sph), 0.0) def test_Mises(self): """Ensure that equivalent stress is 3/2 of equivalent strain.""" x = np.random.rand(self.n,3,3) - assert np.allclose(mechanics.Mises_stress(x)/mechanics.Mises_strain(x), + assert np.allclose(mechanics.equivalent_stress_Mises(x)/mechanics.equivalent_strain_Mises(x), 1.5) def test_spherical_no_shear(self): From 9a1e8e3c381d3b3a648688833eb2a6ecdda0d4ba Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 18 Nov 2020 13:28:53 +0100 Subject: [PATCH 21/35] for the transition period --- python/damask/_rotation.py | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index e695ce407..b9d4ae644 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -463,8 +463,8 @@ class Rotation: @staticmethod def from_Euler_angles(phi, - degrees = False, - **kwargs): + degrees = False, + **kwargs): """ Initialize from Bunge-Euler angles. @@ -607,9 +607,9 @@ class Rotation: @staticmethod def from_Rodrigues_vector(rho, - normalize = False, - P = -1, - **kwargs): + normalize = False, + P = -1, + **kwargs): """ Initialize from Rodrigues-Frank vector. @@ -1440,3 +1440,9 @@ class Rotation: np.where(np.maximum(np.abs(xyz[...,1]),np.abs(xyz[...,2])) <= np.abs(xyz[...,0]),1,2)) return order[direction][p] + + +# for compatibility with deprecated tests +Rotation.from_Eulers = Rotation.from_Euler_angles +Rotation.as_Eulers = Rotation.as_Euler_angles +Rotation.from_Rodrigues = Rotation.from_Rodrigues_vector From 1c9028d0b180aad2aa2e3261d15f0b9bfae469de Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 18 Nov 2020 14:01:52 +0100 Subject: [PATCH 22/35] simplified --- python/damask/mechanics.py | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/python/damask/mechanics.py b/python/damask/mechanics.py index bab41a18c..d1a208e3c 100644 --- a/python/damask/mechanics.py +++ b/python/damask/mechanics.py @@ -60,8 +60,7 @@ def stress_Cauchy(P,F): Cauchy stress. """ - sigma = _np.einsum('...,...ij,...kj',1.0/_np.linalg.det(F),P,F) - return tensor.symmetric(sigma) + return tensor.symmetric(_np.einsum('...,...ij,...kj',1.0/_np.linalg.det(F),P,F)) def deviatoric_part(T): @@ -79,7 +78,7 @@ def deviatoric_part(T): Deviatoric part of T. """ - return T - _np.einsum('...ij,...',_np.eye(3),spherical_part(T)) + return T - spherical_part(T,tensor=True) def maximum_shear(T_sym): @@ -157,8 +156,7 @@ def stress_second_Piola_Kirchhoff(P,F): Second Piola-Kirchhoff stress. """ - S = _np.einsum('...jk,...kl',_np.linalg.inv(F),P) - return tensor.symmetric(S) + return tensor.symmetric(_np.einsum('...ij,...jk',_np.linalg.inv(F),P)) def rotational_part(T): From 86dfb58732f489d4068b4782f5ae8ee8efb9a38d Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Wed, 18 Nov 2020 14:18:05 -0500 Subject: [PATCH 23/35] adjust to renaming of from_Eulers/Rodrigues --- PRIVATE | 2 +- processing/post/addOrientations.py | 14 +++++++------- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/PRIVATE b/PRIVATE index 2d00aa541..dd29f2cc9 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 2d00aa541f071dbfc200f32e358d324995a061f5 +Subproject commit dd29f2cc9ad2c5310840fabbfb7b5ab3ffadf07a diff --git a/processing/post/addOrientations.py b/processing/post/addOrientations.py index 6a02cca08..cb5361599 100755 --- a/processing/post/addOrientations.py +++ b/processing/post/addOrientations.py @@ -115,10 +115,10 @@ for name in filenames: if options.eulers is not None: label = options.eulers print(np.max(table.get(options.eulers),axis=0)) - o = damask.Rotation.from_Eulers(table.get(options.eulers), options.degrees) + o = damask.Rotation.from_Euler_angles(table.get(options.eulers), options.degrees) elif options.rodrigues is not None: label = options.rodrigues - o = damask.Rotation.from_Rodrigues(table.get(options.rodrigues)) + o = damask.Rotation.from_Rodrigues_vector(table.get(options.rodrigues)) elif options.matrix is not None: label = options.matrix o = damask.Rotation.from_matrix(table.get(options.matrix).reshape(-1,3,3)) @@ -137,14 +137,14 @@ for name in filenames: if 'rodrigues' in options.output: - table = table.add('ro({})'.format(label),o.as_Rodrigues(), scriptID+' '+' '.join(sys.argv[1:])) + table = table.add('ro({})'.format(label),o.as_Rodrigues_vector(), scriptID+' '+' '.join(sys.argv[1:])) if 'eulers' in options.output: - table = table.add('eu({})'.format(label),o.as_Eulers(options.degrees), scriptID+' '+' '.join(sys.argv[1:])) + table = table.add('eu({})'.format(label),o.as_Euler_angles(options.degrees),scriptID+' '+' '.join(sys.argv[1:])) if 'quaternion' in options.output: - table = table.add('qu({})'.format(label),o.as_quaternion(), scriptID+' '+' '.join(sys.argv[1:])) + table = table.add('qu({})'.format(label),o.as_quaternion(), scriptID+' '+' '.join(sys.argv[1:])) if 'matrix' in options.output: - table = table.add('om({})'.format(label),o.as_matrix(), scriptID+' '+' '.join(sys.argv[1:])) + table = table.add('om({})'.format(label),o.as_matrix(), scriptID+' '+' '.join(sys.argv[1:])) if 'axisangle' in options.output: - table = table.add('om({})'.format(label),o.as_axisangle(options.degrees), scriptID+' '+' '.join(sys.argv[1:])) + table = table.add('om({})'.format(label),o.as_axis_angle(options.degrees), scriptID+' '+' '.join(sys.argv[1:])) table.save((sys.stdout if name is None else name), legacy=True) From c74e57f225d085f20db17631d16725a3fce82ec3 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 18 Nov 2020 16:18:57 +0100 Subject: [PATCH 24/35] vtk error handling is not very helpful --- python/damask/_vtk.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/python/damask/_vtk.py b/python/damask/_vtk.py index 942be4f1f..8104c3cca 100644 --- a/python/damask/_vtk.py +++ b/python/damask/_vtk.py @@ -1,3 +1,4 @@ +import os import multiprocessing as mp from pathlib import Path @@ -165,6 +166,9 @@ class VTK: else: raise TypeError(f'Unknown file extension {ext}') + if not os.path.isfile(fname): + raise FileNotFoundError(f'No such file: {fname}') + reader.SetFileName(str(fname)) reader.Update() vtk_data = reader.GetOutput() From 903c185ee642d4212943293c3863a7b1b69ce65d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 19 Nov 2020 10:39:41 +0100 Subject: [PATCH 25/35] distinguish 'file not found' and 'invalid file' --- python/damask/_vtk.py | 5 ++--- python/damask/mechanics.py | 3 +++ python/tests/test_Geom.py | 2 +- python/tests/test_VTK.py | 15 ++++++++++----- 4 files changed, 16 insertions(+), 9 deletions(-) diff --git a/python/damask/_vtk.py b/python/damask/_vtk.py index 8104c3cca..8e6cf490c 100644 --- a/python/damask/_vtk.py +++ b/python/damask/_vtk.py @@ -139,6 +139,8 @@ class VTK: vtkUnstructuredGrid, and vtkPolyData. """ + if not os.path.isfile(fname): # vtk has a strange error handling + raise FileNotFoundError(f'No such file: {fname}') ext = Path(fname).suffix if ext == '.vtk' or dataset_type is not None: reader = vtk.vtkGenericDataObjectReader() @@ -166,9 +168,6 @@ class VTK: else: raise TypeError(f'Unknown file extension {ext}') - if not os.path.isfile(fname): - raise FileNotFoundError(f'No such file: {fname}') - reader.SetFileName(str(fname)) reader.Update() vtk_data = reader.GetOutput() diff --git a/python/damask/mechanics.py b/python/damask/mechanics.py index d1a208e3c..c8d865885 100644 --- a/python/damask/mechanics.py +++ b/python/damask/mechanics.py @@ -299,6 +299,9 @@ def _polar_decomposition(T,requested): if 'U' in requested: output.append(_np.einsum('...ji,...jk',R,T)) + if len(output) == 0: + raise ValueError('Output needs to be out of V,R,U') + return tuple(output) diff --git a/python/tests/test_Geom.py b/python/tests/test_Geom.py index 3b7ef61ce..10220559b 100644 --- a/python/tests/test_Geom.py +++ b/python/tests/test_Geom.py @@ -55,7 +55,7 @@ class TestGeom: def test_invalid_vtr(self,tmp_path): v = VTK.from_rectilinear_grid(np.random.randint(5,10,3)*2,np.random.random(3) + 1.0) - v.save(tmp_path/'no_materialpoint.vtr') + v.save(tmp_path/'no_materialpoint.vtr',parallel=False) with pytest.raises(ValueError): Geom.load(tmp_path/'no_materialpoint.vtr') diff --git a/python/tests/test_VTK.py b/python/tests/test_VTK.py index 48dce4707..af4ecf918 100644 --- a/python/tests/test_VTK.py +++ b/python/tests/test_VTK.py @@ -100,12 +100,17 @@ class TestVTK: v = VTK.from_poly_data(points) v.save(tmp_path/fname) - @pytest.mark.parametrize('name,dataset_type',[('this_file_does_not_exist.vtk', None), - ('this_file_does_not_exist.vtk','vtk'), - ('this_file_does_not_exist.vtx', None)]) - def test_invalid_dataset_type(self,name,dataset_type): + @pytest.mark.parametrize('fname,dataset_type',[('a_file.vtk', None), + ('a_file.vtk','vtk'), + ('a_file.vtx', None)]) + def test_invalid_dataset_type(self,tmp_path,fname,dataset_type): + open(tmp_path/fname,'a').close() with pytest.raises(TypeError): - VTK.load(name,dataset_type) + VTK.load(tmp_path/fname,dataset_type) + + def test_file_not_found(self): + with pytest.raises(FileNotFoundError): + VTK.load('/dev/null') def test_add_extension(self,tmp_path,default): default.save(tmp_path/'default.txt',parallel=False) From 5a5dd24687020964f1fe9f6fbeee42ea3329c1f2 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 19 Nov 2020 10:50:14 +0100 Subject: [PATCH 26/35] removed aliases --- PRIVATE | 2 +- python/damask/_result.py | 3 --- python/damask/mechanics.py | 4 ---- 3 files changed, 1 insertion(+), 8 deletions(-) diff --git a/PRIVATE b/PRIVATE index dd29f2cc9..3a4f76fee 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit dd29f2cc9ad2c5310840fabbfb7b5ab3ffadf07a +Subproject commit 3a4f76fee6cbaa8520190fcab1845713adcb1f85 diff --git a/python/damask/_result.py b/python/damask/_result.py index c1753113e..cae3298f7 100644 --- a/python/damask/_result.py +++ b/python/damask/_result.py @@ -1296,6 +1296,3 @@ class Result: v.add(u,'u') v.save(f'{self.fname.stem}_inc{inc[3:].zfill(N_digits)}') - - -Result.add_PK2 = Result.add_stress_second_Piola_Kirchhoff diff --git a/python/damask/mechanics.py b/python/damask/mechanics.py index c8d865885..37e21945f 100644 --- a/python/damask/mechanics.py +++ b/python/damask/mechanics.py @@ -319,7 +319,3 @@ def _equivalent_Mises(T_sym,s): """ d = deviatoric_part(T_sym) return _np.sqrt(s*_np.sum(d**2.0,axis=(-1,-2))) - - -# for compatibility -strain_tensor = strain From 1c07152b96cb700e31c27e3b2f84923f44e3c576 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 19 Nov 2020 14:05:59 +0100 Subject: [PATCH 27/35] sorted alphabetically --- python/damask/mechanics.py | 127 +++++++++++++++++++------------------ python/damask/tensor.py | 72 ++++++++++----------- 2 files changed, 100 insertions(+), 99 deletions(-) diff --git a/python/damask/mechanics.py b/python/damask/mechanics.py index 37e21945f..1531c66d3 100644 --- a/python/damask/mechanics.py +++ b/python/damask/mechanics.py @@ -1,6 +1,6 @@ """Finite-strain continuum mechanics.""" -from . import tensor +from . import tensor as _tensor import numpy as _np @@ -17,10 +17,10 @@ def deformation_Cauchy_Green_left(F): Returns ------- B : numpy.ndarray of shape (...,3,3) - Left Cauchy-Green deformation tensor. + Left Cauchy-Green deformation _tensor. """ - return _np.matmul(F,tensor.transpose(F)) + return _np.matmul(F,_tensor.transpose(F)) def deformation_Cauchy_Green_right(F): @@ -35,32 +35,10 @@ def deformation_Cauchy_Green_right(F): Returns ------- C : numpy.ndarray of shape (...,3,3) - Right Cauchy-Green deformation tensor. + Right Cauchy-Green deformation _tensor. """ - return _np.matmul(tensor.transpose(F),F) - - -def stress_Cauchy(P,F): - """ - Calculate the Cauchy stress (true stress). - - Resulting tensor is symmetrized as the Cauchy stress needs to be symmetric. - - Parameters - ---------- - P : numpy.ndarray of shape (...,3,3) - First Piola-Kirchhoff stress. - F : numpy.ndarray of shape (...,3,3) - Deformation gradient. - - Returns - ------- - sigma : numpy.ndarray of shape (...,3,3) - Cauchy stress. - - """ - return tensor.symmetric(_np.einsum('...,...ij,...kj',1.0/_np.linalg.det(F),P,F)) + return _np.matmul(_tensor.transpose(F),F) def deviatoric_part(T): @@ -81,25 +59,6 @@ def deviatoric_part(T): return T - spherical_part(T,tensor=True) -def maximum_shear(T_sym): - """ - Calculate the maximum shear component of a symmetric tensor. - - Parameters - ---------- - T_sym : numpy.ndarray of shape (...,3,3) - Symmetric tensor of which the maximum shear is computed. - - Returns - ------- - gamma_max : numpy.ndarray of shape (...) - Maximum shear of T_sym. - - """ - w = tensor.eigenvalues(T_sym) - return (w[...,0] - w[...,2])*0.5 - - def equivalent_strain_Mises(epsilon): """ Calculate the Mises equivalent of a strain tensor. @@ -136,27 +95,23 @@ def equivalent_stress_Mises(sigma): return _equivalent_Mises(sigma,3.0/2.0) -def stress_second_Piola_Kirchhoff(P,F): +def maximum_shear(T_sym): """ - Calculate the second Piola-Kirchhoff stress. - - Resulting tensor is symmetrized as the second Piola-Kirchhoff stress - needs to be symmetric. + Calculate the maximum shear component of a symmetric tensor. Parameters ---------- - P : numpy.ndarray of shape (...,3,3) - First Piola-Kirchhoff stress. - F : numpy.ndarray of shape (...,3,3) - Deformation gradient. + T_sym : numpy.ndarray of shape (...,3,3) + Symmetric tensor of which the maximum shear is computed. Returns ------- - S : numpy.ndarray of shape (...,3,3) - Second Piola-Kirchhoff stress. + gamma_max : numpy.ndarray of shape (...) + Maximum shear of T_sym. """ - return tensor.symmetric(_np.einsum('...ij,...jk',_np.linalg.inv(F),P)) + w = _tensor.eigenvalues(T_sym) + return (w[...,0] - w[...,2])*0.5 def rotational_part(T): @@ -186,14 +141,14 @@ def spherical_part(T,tensor=False): T : numpy.ndarray of shape (...,3,3) Tensor of which the hydrostatic part is computed. tensor : bool, optional - Map spherical part onto identity tensor. Defaults to false + Map spherical part onto identity _tensor. Defaults to false Returns ------- p : numpy.ndarray of shape (...) unless tensor == True: shape (...,3,3) Spherical part of tensor T, e.g. the hydrostatic part/pressure - of a stress tensor. + of a stress _tensor. """ sph = _np.trace(T,axis2=-2,axis1=-1)/3.0 @@ -213,7 +168,7 @@ def strain(F,t,m): Deformation gradient. t : {‘V’, ‘U’} Type of the polar decomposition, ‘V’ for left stretch tensor - and ‘U’ for right stretch tensor. + and ‘U’ for right stretch _tensor. m : float Order of the strain. @@ -239,9 +194,55 @@ def strain(F,t,m): return eps + +def stress_Cauchy(P,F): + """ + Calculate the Cauchy stress (true stress). + + Resulting tensor is symmetrized as the Cauchy stress needs to be symmetric. + + Parameters + ---------- + P : numpy.ndarray of shape (...,3,3) + First Piola-Kirchhoff stress. + F : numpy.ndarray of shape (...,3,3) + Deformation gradient. + + Returns + ------- + sigma : numpy.ndarray of shape (...,3,3) + Cauchy stress. + + """ + return _tensor.symmetric(_np.einsum('...,...ij,...kj',1.0/_np.linalg.det(F),P,F)) + + +def stress_second_Piola_Kirchhoff(P,F): + """ + Calculate the second Piola-Kirchhoff stress. + + Resulting tensor is symmetrized as the second Piola-Kirchhoff stress + needs to be symmetric. + + Parameters + ---------- + P : numpy.ndarray of shape (...,3,3) + First Piola-Kirchhoff stress. + F : numpy.ndarray of shape (...,3,3) + Deformation gradient. + + Returns + ------- + S : numpy.ndarray of shape (...,3,3) + Second Piola-Kirchhoff stress. + + """ + return _tensor.symmetric(_np.einsum('...ij,...jk',_np.linalg.inv(F),P)) + + def stretch_left(T): """ - Calculate left stretch of a tensor. + Calculate left stretch of a _tensor. Parameters ---------- @@ -259,7 +260,7 @@ def stretch_left(T): def stretch_right(T): """ - Calculate right stretch of a tensor. + Calculate right stretch of a _tensor. Parameters ---------- diff --git a/python/damask/tensor.py b/python/damask/tensor.py index 6e5bdaf85..4b22560b6 100644 --- a/python/damask/tensor.py +++ b/python/damask/tensor.py @@ -11,42 +11,6 @@ to operate on numpy.ndarrays of shape (...,3,3). import numpy as _np -def symmetric(T): - """ - Symmetrize tensor. - - Parameters - ---------- - T : numpy.ndarray of shape (...,3,3) - Tensor of which the symmetrized values are computed. - - Returns - ------- - T_sym : numpy.ndarray of shape (...,3,3) - Symmetrized tensor T. - - """ - return (T+transpose(T))*0.5 - - -def transpose(T): - """ - Transpose tensor. - - Parameters - ---------- - T : numpy.ndarray of shape (...,3,3) - Tensor of which the transpose is computed. - - Returns - ------- - T.T : numpy.ndarray of shape (...,3,3) - Transpose of tensor T. - - """ - return _np.swapaxes(T,axis2=-2,axis1=-1) - - def eigenvalues(T_sym): """ Eigenvalues, i.e. principal components, of a symmetric tensor. @@ -89,3 +53,39 @@ def eigenvectors(T_sym,RHS=False): if RHS: v[_np.linalg.det(v) < 0.0,:,2] *= -1.0 return v + + +def symmetric(T): + """ + Symmetrize tensor. + + Parameters + ---------- + T : numpy.ndarray of shape (...,3,3) + Tensor of which the symmetrized values are computed. + + Returns + ------- + T_sym : numpy.ndarray of shape (...,3,3) + Symmetrized tensor T. + + """ + return (T+transpose(T))*0.5 + + +def transpose(T): + """ + Transpose tensor. + + Parameters + ---------- + T : numpy.ndarray of shape (...,3,3) + Tensor of which the transpose is computed. + + Returns + ------- + T.T : numpy.ndarray of shape (...,3,3) + Transpose of tensor T. + + """ + return _np.swapaxes(T,axis2=-2,axis1=-1) From 894a8de9f9a163082a81547a3cc8229fbc7e8500 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 19 Nov 2020 14:31:14 +0100 Subject: [PATCH 28/35] avoid name clash: as_Rodrigues_vector(vector = ...) --- python/damask/_orientation.py | 4 ++-- python/damask/_rotation.py | 9 ++++----- 2 files changed, 6 insertions(+), 7 deletions(-) diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index 5f0889135..5595a9693 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -469,7 +469,7 @@ class Orientation(Rotation): if self.family is None: raise ValueError('Missing crystal symmetry') - rho_abs = np.abs(self.as_Rodrigues_vector(vector=True)) + rho_abs = np.abs(self.as_Rodrigues_vector(compact=True)) with np.errstate(invalid='ignore'): # using '*'/prod for 'and' @@ -512,7 +512,7 @@ class Orientation(Rotation): if self.family is None: raise ValueError('Missing crystal symmetry') - rho = self.as_Rodrigues_vector(vector=True) + rho = self.as_Rodrigues_vector(compact=True) with np.errstate(invalid='ignore'): if self.family == 'cubic': diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index a5e534244..fc713401b 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -357,7 +357,7 @@ class Rotation: return Rotation._qu2om(self.quaternion) def as_Rodrigues_vector(self, - vector = False): + compact = False): """ Represent as Rodrigues-Frank vector with separated axis and angle argument. @@ -375,7 +375,7 @@ class Rotation: """ ro = Rotation._qu2ro(self.quaternion) - if vector: + if compact: with np.errstate(invalid='ignore'): return ro[...,:3]*ro[...,3:4] else: @@ -595,13 +595,12 @@ class Rotation: P = -1, **kwargs): """ - Initialize from Rodrigues-Frank vector. + Initialize from Rodrigues-Frank vector (angle separated from axis). Parameters ---------- rho : numpy.ndarray of shape (...,4) - Rodrigues-Frank vector (angle separated from axis). - (n_1, n_2, n_3, tan(ω/2)), ǀnǀ = 1 and ω ∈ [0,π]. + Rodrigues-Frank vector. (n_1, n_2, n_3, tan(ω/2)), ǀnǀ = 1 and ω ∈ [0,π]. normalize : boolean, optional Allow ǀnǀ ≠ 1. Defaults to False. P : int ∈ {-1,1}, optional From 6e5cb6013252ddebf037f6e7eb42cc4125e874b4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 19 Nov 2020 14:38:54 +0100 Subject: [PATCH 29/35] general tensor functionality --- python/damask/__init__.py | 10 +++++- python/damask/_result.py | 4 +-- python/damask/mechanics.py | 53 ++++------------------------- python/damask/tensor.py | 41 +++++++++++++++++++++++ python/tests/test_Result.py | 4 +-- python/tests/test_mechanics.py | 61 ++++++---------------------------- python/tests/test_tensor.py | 34 ++++++++++++++++--- 7 files changed, 101 insertions(+), 106 deletions(-) diff --git a/python/damask/__init__.py b/python/damask/__init__.py index f3d5ad2fb..91ca09b29 100644 --- a/python/damask/__init__.py +++ b/python/damask/__init__.py @@ -1,4 +1,12 @@ -"""Tools for pre and post processing of DAMASK simulations.""" +""" +Tools for pre and post processing of DAMASK simulations. + +Modules that contain only one class (of the same name), +are prefixed by a '_'. For example, '_colormap' contains +a class called 'Colormap' which is imported as 'damask.Colormap'. + +""" + from pathlib import Path as _Path import re as _re diff --git a/python/damask/_result.py b/python/damask/_result.py index cae3298f7..8bacebb90 100644 --- a/python/damask/_result.py +++ b/python/damask/_result.py @@ -649,7 +649,7 @@ class Result: @staticmethod def _add_deviator(T): return { - 'data': mechanics.deviatoric_part(T['data']), + 'data': tensor.deviatoric(T['data']), 'label': f"s_{T['label']}", 'meta': { 'Unit': T['meta']['Unit'], @@ -968,7 +968,7 @@ class Result: @staticmethod def _add_spherical(T): return { - 'data': mechanics.spherical_part(T['data']), + 'data': tensor.spherical(T['data'],False), 'label': f"p_{T['label']}", 'meta': { 'Unit': T['meta']['Unit'], diff --git a/python/damask/mechanics.py b/python/damask/mechanics.py index 1531c66d3..2be93b91b 100644 --- a/python/damask/mechanics.py +++ b/python/damask/mechanics.py @@ -17,7 +17,7 @@ def deformation_Cauchy_Green_left(F): Returns ------- B : numpy.ndarray of shape (...,3,3) - Left Cauchy-Green deformation _tensor. + Left Cauchy-Green deformation tensor. """ return _np.matmul(F,_tensor.transpose(F)) @@ -35,30 +35,12 @@ def deformation_Cauchy_Green_right(F): Returns ------- C : numpy.ndarray of shape (...,3,3) - Right Cauchy-Green deformation _tensor. + Right Cauchy-Green deformation tensor. """ return _np.matmul(_tensor.transpose(F),F) -def deviatoric_part(T): - """ - Calculate deviatoric part of a tensor. - - Parameters - ---------- - T : numpy.ndarray of shape (...,3,3) - Tensor of which the deviatoric part is computed. - - Returns - ------- - T' : numpy.ndarray of shape (...,3,3) - Deviatoric part of T. - - """ - return T - spherical_part(T,tensor=True) - - def equivalent_strain_Mises(epsilon): """ Calculate the Mises equivalent of a strain tensor. @@ -132,29 +114,6 @@ def rotational_part(T): return _polar_decomposition(T,'R')[0] -def spherical_part(T,tensor=False): - """ - Calculate spherical (hydrostatic) part of a tensor. - - Parameters - ---------- - T : numpy.ndarray of shape (...,3,3) - Tensor of which the hydrostatic part is computed. - tensor : bool, optional - Map spherical part onto identity _tensor. Defaults to false - - Returns - ------- - p : numpy.ndarray of shape (...) - unless tensor == True: shape (...,3,3) - Spherical part of tensor T, e.g. the hydrostatic part/pressure - of a stress _tensor. - - """ - sph = _np.trace(T,axis2=-2,axis1=-1)/3.0 - return _np.einsum('...jk,...',_np.eye(3),sph) if tensor else sph - - def strain(F,t,m): """ Calculate strain tensor (Seth–Hill family). @@ -168,7 +127,7 @@ def strain(F,t,m): Deformation gradient. t : {‘V’, ‘U’} Type of the polar decomposition, ‘V’ for left stretch tensor - and ‘U’ for right stretch _tensor. + and ‘U’ for right stretch tensor. m : float Order of the strain. @@ -242,7 +201,7 @@ def stress_second_Piola_Kirchhoff(P,F): def stretch_left(T): """ - Calculate left stretch of a _tensor. + Calculate left stretch of a tensor. Parameters ---------- @@ -260,7 +219,7 @@ def stretch_left(T): def stretch_right(T): """ - Calculate right stretch of a _tensor. + Calculate right stretch of a tensor. Parameters ---------- @@ -318,5 +277,5 @@ def _equivalent_Mises(T_sym,s): Scaling factor (2/3 for strain, 3/2 for stress). """ - d = deviatoric_part(T_sym) + d = _tensor.deviatoric(T_sym) return _np.sqrt(s*_np.sum(d**2.0,axis=(-1,-2))) diff --git a/python/damask/tensor.py b/python/damask/tensor.py index 4b22560b6..3ee551d02 100644 --- a/python/damask/tensor.py +++ b/python/damask/tensor.py @@ -11,6 +11,24 @@ to operate on numpy.ndarrays of shape (...,3,3). import numpy as _np +def deviatoric(T): + """ + Calculate deviatoric part of a tensor. + + Parameters + ---------- + T : numpy.ndarray of shape (...,3,3) + Tensor of which the deviatoric part is computed. + + Returns + ------- + T' : numpy.ndarray of shape (...,3,3) + Deviatoric part of T. + + """ + return T - spherical(T,tensor=True) + + def eigenvalues(T_sym): """ Eigenvalues, i.e. principal components, of a symmetric tensor. @@ -55,6 +73,29 @@ def eigenvectors(T_sym,RHS=False): return v +def spherical(T,tensor=True): + """ + Calculate spherical part of a tensor. + + + Parameters + ---------- + T : numpy.ndarray of shape (...,3,3) + Tensor of which the spherical part is computed. + tensor : bool, optional + Map spherical part onto identity tensor. Defaults to True. + + Returns + ------- + p : numpy.ndarray of shape (...,3,3) + unless tensor == False: shape (...,) + Spherical part of tensor T. p is an isotropic tensor. + + """ + sph = _np.trace(T,axis2=-2,axis1=-1)/3.0 + return _np.einsum('...jk,...',_np.eye(3),sph) if tensor else sph + + def symmetric(T): """ Symmetrize tensor. diff --git a/python/tests/test_Result.py b/python/tests/test_Result.py index 46b947cdc..fb79f4db7 100644 --- a/python/tests/test_Result.py +++ b/python/tests/test_Result.py @@ -143,7 +143,7 @@ class TestResult: default.add_deviator('P') loc = {'P' :default.get_dataset_location('P'), 's_P':default.get_dataset_location('s_P')} - in_memory = mechanics.deviatoric_part(default.read_dataset(loc['P'],0)) + in_memory = tensor.deviatoric(default.read_dataset(loc['P'],0)) in_file = default.read_dataset(loc['s_P'],0) assert np.allclose(in_memory,in_file) @@ -273,7 +273,7 @@ class TestResult: default.add_spherical('P') loc = {'P': default.get_dataset_location('P'), 'p_P': default.get_dataset_location('p_P')} - in_memory = mechanics.spherical_part(default.read_dataset(loc['P'],0)).reshape(-1,1) + in_memory = tensor.spherical(default.read_dataset(loc['P'],0),False).reshape(-1,1) in_file = default.read_dataset(loc['p_P'],0) assert np.allclose(in_memory,in_file) diff --git a/python/tests/test_mechanics.py b/python/tests/test_mechanics.py index 4a61f5b1e..4912aa962 100644 --- a/python/tests/test_mechanics.py +++ b/python/tests/test_mechanics.py @@ -9,10 +9,6 @@ def stress_Cauchy(P,F): return symmetric(sigma) -def deviatoric_part(T): - return T - np.eye(3)*spherical_part(T) - - def eigenvalues(T_sym): return np.linalg.eigvalsh(symmetric(T_sym)) @@ -39,11 +35,6 @@ def rotational_part(T): return polar_decomposition(T,'R')[0] -def spherical_part(T,tensor=False): - sph = np.trace(T)/3.0 - return sph if not tensor else np.eye(3)*sph - - def strain(F,t,m): if t == 'V': @@ -92,31 +83,20 @@ def polar_decomposition(T,requested): return tuple(output) def equivalent_Mises(T_sym,s): - return np.sqrt(s*(np.sum(deviatoric_part(T_sym)**2.0))) + return np.sqrt(s*(np.sum(deviatoric(T_sym)**2.0))) +def deviatoric(T): + return T - np.eye(3)*np.trace(T)/3.0 class TestMechanics: n = 1000 c = np.random.randint(n) - - @pytest.mark.parametrize('vectorized,single',[(mechanics.deviatoric_part, deviatoric_part), - (mechanics.spherical_part, spherical_part) - ]) - def test_vectorize_1_arg_(self,vectorized,single): - print("done") - test_data_flat = np.random.rand(self.n,3,3) - test_data = np.reshape(test_data_flat,(self.n//10,10,3,3)) - for i,v in enumerate(np.reshape(vectorized(test_data),vectorized(test_data_flat).shape)): - assert np.allclose(single(test_data_flat[i]),v) - - @pytest.mark.parametrize('vectorized,single',[(mechanics.deviatoric_part, deviatoric_part), - (mechanics.maximum_shear, maximum_shear), + @pytest.mark.parametrize('vectorized,single',[(mechanics.maximum_shear, maximum_shear), (mechanics.equivalent_stress_Mises, equivalent_stress_Mises), (mechanics.equivalent_strain_Mises, equivalent_strain_Mises), (mechanics.rotational_part, rotational_part), - (mechanics.spherical_part, spherical_part), (mechanics.stretch_left, stretch_left), (mechanics.stretch_right, stretch_right), ]) @@ -155,11 +135,6 @@ class TestMechanics: P = np.random.rand(self.n,3,3) assert np.allclose(function(P,np.broadcast_to(np.eye(3),(self.n,3,3))),tensor.symmetric(P)) - def test_deviatoric_part(self): - I_n = np.broadcast_to(np.eye(3),(self.n,3,3)) - r = np.logical_not(I_n)*np.random.rand(self.n,3,3) - assert np.allclose(mechanics.deviatoric_part(I_n+r),r) - def test_polar_decomposition(self): """F = RU = VR.""" F = np.broadcast_to(np.eye(3),[self.n,3,3])*np.random.rand(self.n,3,3) @@ -201,34 +176,20 @@ class TestMechanics: assert np.allclose(np.abs(np.linalg.det(mechanics.rotational_part(x))), 1.0) - def test_spherical_deviatoric_part(self): - """Ensure that full tensor is sum of spherical and deviatoric part.""" - x = np.random.rand(self.n,3,3) - sph = mechanics.spherical_part(x,True) - assert np.allclose(sph + mechanics.deviatoric_part(x), - x) - def test_deviatoric_Mises(self): """Ensure that Mises equivalent stress depends only on deviatoric part.""" x = np.random.rand(self.n,3,3) full = mechanics.equivalent_stress_Mises(x) - dev = mechanics.equivalent_stress_Mises(mechanics.deviatoric_part(x)) + dev = mechanics.equivalent_stress_Mises(tensor.deviatoric(x)) assert np.allclose(full, dev) - def test_spherical_mapping(self): - """Ensure that mapping to tensor is correct.""" + @pytest.mark.parametrize('Mises_equivalent',[mechanics.equivalent_strain_Mises, + mechanics.equivalent_stress_Mises]) + def test_spherical_Mises(self,Mises_equivalent): + """Ensure that Mises equivalent strain/stress of spherical strain is 0.""" x = np.random.rand(self.n,3,3) - tnsr = mechanics.spherical_part(x,True) - scalar = mechanics.spherical_part(x) - assert np.allclose(np.linalg.det(tnsr), - scalar**3.0) - - def test_spherical_Mises(self): - """Ensure that Mises equivalent strain of spherical strain is 0.""" - x = np.random.rand(self.n,3,3) - sph = mechanics.spherical_part(x,True) - assert np.allclose(mechanics.equivalent_strain_Mises(sph), + assert np.allclose(Mises_equivalent(tensor.spherical(x,True)), 0.0) @@ -240,5 +201,5 @@ class TestMechanics: def test_spherical_no_shear(self): """Ensure that sherical stress has max shear of 0.0.""" - A = mechanics.spherical_part(tensor.symmetric(np.random.rand(self.n,3,3)),True) + A = tensor.spherical(tensor.symmetric(np.random.rand(self.n,3,3)),True) assert np.allclose(mechanics.maximum_shear(A),0.0) diff --git a/python/tests/test_tensor.py b/python/tests/test_tensor.py index 24fbe6126..58194130e 100644 --- a/python/tests/test_tensor.py +++ b/python/tests/test_tensor.py @@ -3,6 +3,8 @@ import numpy as np from damask import tensor +def deviatoric(T): + return T - spherical(T) def eigenvalues(T_sym): return np.linalg.eigvalsh(symmetric(T_sym)) @@ -20,6 +22,10 @@ def symmetric(T): def transpose(T): return T.T +def spherical(T,tensor=True): + sph = np.trace(T)/3.0 + return sph if not tensor else np.eye(3)*sph + class TestTensor: @@ -27,10 +33,12 @@ class TestTensor: c = np.random.randint(n) - @pytest.mark.parametrize('vectorized,single',[(tensor.eigenvalues , eigenvalues ), - (tensor.eigenvectors , eigenvectors ), - (tensor.symmetric , symmetric ), - (tensor.transpose , transpose ), + @pytest.mark.parametrize('vectorized,single',[(tensor.deviatoric, deviatoric), + (tensor.eigenvalues, eigenvalues), + (tensor.eigenvectors, eigenvectors), + (tensor.symmetric, symmetric), + (tensor.transpose, transpose), + (tensor.spherical, spherical), ]) def test_vectorize_1_arg(self,vectorized,single): epsilon = np.random.rand(self.n,3,3) @@ -71,3 +79,21 @@ class TestTensor: LRHS = np.linalg.det(tensor.eigenvectors(A,RHS=False)) RHS = np.linalg.det(tensor.eigenvectors(A,RHS=True)) assert np.allclose(np.abs(LRHS),RHS) + + def test_spherical_deviatoric_part(self): + """Ensure that full tensor is sum of spherical and deviatoric part.""" + x = np.random.rand(self.n,3,3) + assert np.allclose(tensor.spherical(x,True) + tensor.deviatoric(x), + x) + def test_spherical_mapping(self): + """Ensure that mapping to tensor is correct.""" + x = np.random.rand(self.n,3,3) + tnsr = tensor.spherical(x,True) + scalar = tensor.spherical(x,False) + assert np.allclose(np.linalg.det(tnsr), + scalar**3.0) + + def test_deviatoric(self): + I_n = np.broadcast_to(np.eye(3),(self.n,3,3)) + r = np.logical_not(I_n)*np.random.rand(self.n,3,3) + assert np.allclose(tensor.deviatoric(I_n+r),r) From a4b5c2a537461b7c2842a988ddbe1fcf58e8c048 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 19 Nov 2020 15:19:15 +0100 Subject: [PATCH 30/35] shorter (but still self-explanatory) name --- python/damask/_result.py | 10 +++++----- python/damask/mechanics.py | 2 +- python/damask/tensor.py | 1 - python/tests/test_Result.py | 6 +++--- python/tests/test_mechanics.py | 10 +++++----- 5 files changed, 14 insertions(+), 15 deletions(-) diff --git a/python/damask/_result.py b/python/damask/_result.py index 8bacebb90..99eed6d00 100644 --- a/python/damask/_result.py +++ b/python/damask/_result.py @@ -942,17 +942,17 @@ class Result: @staticmethod - def _add_rotational_part(F): + def _add_rotational(F): return { - 'data': mechanics.rotational_part(F['data']), + 'data': mechanics.rotational(F['data']), 'label': f"R({F['label']})", 'meta': { 'Unit': F['meta']['Unit'], 'Description': f"Rotational part of {F['label']} ({F['meta']['Description']})", - 'Creator': 'add_rotational_part' + 'Creator': 'add_rotational' } } - def add_rotational_part(self,F): + def add_rotational(self,F): """ Add rotational part of a deformation gradient. @@ -962,7 +962,7 @@ class Result: Label of deformation gradient dataset. """ - self._add_generic_pointwise(self._add_rotational_part,{'F':F}) + self._add_generic_pointwise(self._add_rotational,{'F':F}) @staticmethod diff --git a/python/damask/mechanics.py b/python/damask/mechanics.py index 2be93b91b..817d5aa9c 100644 --- a/python/damask/mechanics.py +++ b/python/damask/mechanics.py @@ -96,7 +96,7 @@ def maximum_shear(T_sym): return (w[...,0] - w[...,2])*0.5 -def rotational_part(T): +def rotational(T): """ Calculate the rotational part of a tensor. diff --git a/python/damask/tensor.py b/python/damask/tensor.py index 3ee551d02..00fc7e72a 100644 --- a/python/damask/tensor.py +++ b/python/damask/tensor.py @@ -77,7 +77,6 @@ def spherical(T,tensor=True): """ Calculate spherical part of a tensor. - Parameters ---------- T : numpy.ndarray of shape (...,3,3) diff --git a/python/tests/test_Result.py b/python/tests/test_Result.py index fb79f4db7..dea66e1ef 100644 --- a/python/tests/test_Result.py +++ b/python/tests/test_Result.py @@ -261,11 +261,11 @@ class TestResult: in_file = default.read_dataset(loc['pole']) assert np.allclose(in_memory,in_file) - def test_add_rotational_part(self,default): - default.add_rotational_part('F') + def test_add_rotational(self,default): + default.add_rotational('F') loc = {'F': default.get_dataset_location('F'), 'R(F)': default.get_dataset_location('R(F)')} - in_memory = mechanics.rotational_part(default.read_dataset(loc['F'],0)) + in_memory = mechanics.rotational(default.read_dataset(loc['F'],0)) in_file = default.read_dataset(loc['R(F)'],0) assert np.allclose(in_memory,in_file) diff --git a/python/tests/test_mechanics.py b/python/tests/test_mechanics.py index 4912aa962..94ab042ed 100644 --- a/python/tests/test_mechanics.py +++ b/python/tests/test_mechanics.py @@ -31,7 +31,7 @@ def stress_second_Piola_Kirchhoff(P,F): return symmetric(S) -def rotational_part(T): +def rotational(T): return polar_decomposition(T,'R')[0] @@ -96,7 +96,7 @@ class TestMechanics: @pytest.mark.parametrize('vectorized,single',[(mechanics.maximum_shear, maximum_shear), (mechanics.equivalent_stress_Mises, equivalent_stress_Mises), (mechanics.equivalent_strain_Mises, equivalent_strain_Mises), - (mechanics.rotational_part, rotational_part), + (mechanics.rotational, rotational), (mechanics.stretch_left, stretch_left), (mechanics.stretch_right, stretch_right), ]) @@ -138,7 +138,7 @@ class TestMechanics: def test_polar_decomposition(self): """F = RU = VR.""" F = np.broadcast_to(np.eye(3),[self.n,3,3])*np.random.rand(self.n,3,3) - R = mechanics.rotational_part(F) + R = mechanics.rotational(F) V = mechanics.stretch_left(F) U = mechanics.stretch_right(F) assert np.allclose(np.matmul(R,U), @@ -162,7 +162,7 @@ class TestMechanics: @pytest.mark.parametrize('t',['V','U']) def test_strain_rotation(self,m,t): """Ensure that pure rotation results in no strain.""" - F = mechanics.rotational_part(np.random.rand(self.n,3,3)) + F = mechanics.rotational(np.random.rand(self.n,3,3)) assert np.allclose(mechanics.strain(F,t,m), 0.0) @@ -173,7 +173,7 @@ class TestMechanics: Should be +1, but random F might contain a reflection. """ x = np.random.rand(self.n,3,3) - assert np.allclose(np.abs(np.linalg.det(mechanics.rotational_part(x))), + assert np.allclose(np.abs(np.linalg.det(mechanics.rotational(x))), 1.0) def test_deviatoric_Mises(self): From a87596cefcdf47a3a66e3af6a7b5f6b4bf2b28be Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 19 Nov 2020 22:36:19 +0100 Subject: [PATCH 31/35] return rotation type (ensures proper rotation) --- python/damask/_orientation.py | 4 ++-- python/damask/_result.py | 2 +- python/damask/_rotation.py | 19 +++++++++++-------- python/damask/mechanics.py | 9 +++++---- python/tests/test_Result.py | 2 +- python/tests/test_Rotation.py | 10 ++++++++++ python/tests/test_mechanics.py | 20 ++++++++++++++++---- 7 files changed, 46 insertions(+), 20 deletions(-) diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index 5595a9693..d09a1fb9d 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -956,9 +956,9 @@ class Orientation(Rotation): """ if self.family is None or other.family is None: - raise ValueError('Missing crystal symmetry') + raise ValueError('missing crystal symmetry') if self.family != other.family: - raise NotImplementedError('Disorientation between different crystal families not supported yet.') + raise NotImplementedError('disorientation between different crystal families') blend = util.shapeblender(self.shape,other.shape) s = self.equivalent diff --git a/python/damask/_result.py b/python/damask/_result.py index 99eed6d00..06598e13b 100644 --- a/python/damask/_result.py +++ b/python/damask/_result.py @@ -944,7 +944,7 @@ class Result: @staticmethod def _add_rotational(F): return { - 'data': mechanics.rotational(F['data']), + 'data': mechanics.rotational(F['data']).as_matrix(), 'label': f"R({F['label']})", 'meta': { 'Unit': F['meta']['Unit'], diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index fc713401b..5fb8109c3 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -61,18 +61,21 @@ class Rotation: elif np.array(rotation).shape[-1] == 4: self.quaternion = np.array(rotation) else: - raise ValueError('"rotation" is neither a Rotation nor a quaternion') + raise TypeError('"rotation" is neither a Rotation nor a quaternion') def __repr__(self): """Represent rotation as unit quaternion, rotation matrix, and Bunge-Euler angles.""" - return 'As quaternions:\n'+str(self.quaternion) \ - if self.quaternion.shape != (4,) else \ - '\n'.join([ - 'Quaternion: (real={:.3f}, imag=<{:+.3f}, {:+.3f}, {:+.3f}>)'.format(*(self.quaternion)), - 'Matrix:\n{}'.format(np.round(self.as_matrix(),8)), - 'Bunge Eulers / deg: ({:3.2f}, {:3.2f}, {:3.2f})'.format(*self.as_Euler_angles(degrees=True)), - ]) + if self == Rotation(): + return 'Rotation()' + else: + return f'Quaternions {self.shape}:\n'+str(self.quaternion) \ + if self.quaternion.shape != (4,) else \ + '\n'.join([ + 'Quaternion: (real={:.3f}, imag=<{:+.3f}, {:+.3f}, {:+.3f}>)'.format(*(self.quaternion)), + 'Matrix:\n{}'.format(np.round(self.as_matrix(),8)), + 'Bunge Eulers / deg: ({:3.2f}, {:3.2f}, {:3.2f})'.format(*self.as_Euler_angles(degrees=True)), + ]) # ToDo: Check difference __copy__ vs __deepcopy__ diff --git a/python/damask/mechanics.py b/python/damask/mechanics.py index 817d5aa9c..c401a3032 100644 --- a/python/damask/mechanics.py +++ b/python/damask/mechanics.py @@ -1,6 +1,7 @@ """Finite-strain continuum mechanics.""" from . import tensor as _tensor +from . import _rotation import numpy as _np @@ -107,11 +108,11 @@ def rotational(T): Returns ------- - R : numpy.ndarray of shape (...,3,3) - Rotational part. + R : damask.Rotation of shape (...) + Rotational part of the vector. """ - return _polar_decomposition(T,'R')[0] + return _rotation.Rotation.from_matrix(_polar_decomposition(T,'R')[0]) def strain(F,t,m): @@ -260,7 +261,7 @@ def _polar_decomposition(T,requested): output.append(_np.einsum('...ji,...jk',R,T)) if len(output) == 0: - raise ValueError('Output needs to be out of V,R,U') + raise ValueError('output needs to be out of V, R, U') return tuple(output) diff --git a/python/tests/test_Result.py b/python/tests/test_Result.py index dea66e1ef..ce34c48de 100644 --- a/python/tests/test_Result.py +++ b/python/tests/test_Result.py @@ -265,7 +265,7 @@ class TestResult: default.add_rotational('F') loc = {'F': default.get_dataset_location('F'), 'R(F)': default.get_dataset_location('R(F)')} - in_memory = mechanics.rotational(default.read_dataset(loc['F'],0)) + in_memory = mechanics.rotational(default.read_dataset(loc['F'],0)).as_matrix() in_file = default.read_dataset(loc['R(F)'],0) assert np.allclose(in_memory,in_file) diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index 352d7cf37..8b99f34ec 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -682,6 +682,9 @@ class TestRotation: for v,o in zip(vec,ori): assert np.allclose(func(v,normalize=True).as_quaternion(),o.as_quaternion()) + def test_invalid_init(self): + with pytest.raises(TypeError): + Rotation(np.ones(3)) @pytest.mark.parametrize('degrees',[True,False]) def test_Eulers(self,set_of_rotations,degrees): @@ -831,6 +834,13 @@ class TestRotation: with pytest.raises(ValueError): function(invalid_shape) + def test_invalid_shape_parallel(self): + invalid_a = np.random.random(np.random.randint(8,32,(3))) + invalid_b = np.random.random(np.random.randint(8,32,(3))) + with pytest.raises(ValueError): + Rotation.from_parallel(invalid_a,invalid_b) + + @pytest.mark.parametrize('fr,to',[(Rotation.from_quaternion,'as_quaternion'), (Rotation.from_axis_angle,'as_axis_angle'), (Rotation.from_Rodrigues_vector, 'as_Rodrigues_vector'), diff --git a/python/tests/test_mechanics.py b/python/tests/test_mechanics.py index 94ab042ed..b4672dbbf 100644 --- a/python/tests/test_mechanics.py +++ b/python/tests/test_mechanics.py @@ -3,6 +3,7 @@ import numpy as np from damask import tensor from damask import mechanics +from damask import Rotation def stress_Cauchy(P,F): sigma = 1.0/np.linalg.det(F) * np.dot(P,F.T) @@ -96,7 +97,6 @@ class TestMechanics: @pytest.mark.parametrize('vectorized,single',[(mechanics.maximum_shear, maximum_shear), (mechanics.equivalent_stress_Mises, equivalent_stress_Mises), (mechanics.equivalent_strain_Mises, equivalent_strain_Mises), - (mechanics.rotational, rotational), (mechanics.stretch_left, stretch_left), (mechanics.stretch_right, stretch_right), ]) @@ -106,6 +106,14 @@ class TestMechanics: for i,v in enumerate(np.reshape(vectorized(epsilon_vec),vectorized(epsilon).shape)): assert np.allclose(single(epsilon[i]),v) + def test_vectorize_rotation(self): + epsilon = Rotation.from_random(self.n).as_matrix() + epsilon_vec = np.reshape(epsilon,(self.n//10,10,3,3)) + for i,v in enumerate(np.reshape(mechanics.rotational(epsilon_vec).as_matrix(), + mechanics.rotational(epsilon).as_matrix().shape)): + assert np.allclose(rotational(epsilon[i]),v) + + @pytest.mark.parametrize('vectorized,single',[(mechanics.stress_Cauchy, stress_Cauchy), (mechanics.stress_second_Piola_Kirchhoff, stress_second_Piola_Kirchhoff) ]) @@ -138,7 +146,7 @@ class TestMechanics: def test_polar_decomposition(self): """F = RU = VR.""" F = np.broadcast_to(np.eye(3),[self.n,3,3])*np.random.rand(self.n,3,3) - R = mechanics.rotational(F) + R = mechanics.rotational(F).as_matrix() V = mechanics.stretch_left(F) U = mechanics.stretch_right(F) assert np.allclose(np.matmul(R,U), @@ -162,7 +170,7 @@ class TestMechanics: @pytest.mark.parametrize('t',['V','U']) def test_strain_rotation(self,m,t): """Ensure that pure rotation results in no strain.""" - F = mechanics.rotational(np.random.rand(self.n,3,3)) + F = Rotation.from_random(self.n).as_matrix() assert np.allclose(mechanics.strain(F,t,m), 0.0) @@ -173,7 +181,7 @@ class TestMechanics: Should be +1, but random F might contain a reflection. """ x = np.random.rand(self.n,3,3) - assert np.allclose(np.abs(np.linalg.det(mechanics.rotational(x))), + assert np.allclose(np.abs(np.linalg.det(mechanics._polar_decomposition(x,'R')[0])), 1.0) def test_deviatoric_Mises(self): @@ -203,3 +211,7 @@ class TestMechanics: """Ensure that sherical stress has max shear of 0.0.""" A = tensor.spherical(tensor.symmetric(np.random.rand(self.n,3,3)),True) assert np.allclose(mechanics.maximum_shear(A),0.0) + + def test_invalid_decomposition(self): + with pytest.raises(ValueError): + mechanics._polar_decomposition(np.random.rand(10,3,3),'A') From 2de3a6a2058b45d69233cc8f414cb3fea36d649b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 19 Nov 2020 22:46:52 +0100 Subject: [PATCH 32/35] fits better to stretch_left/right --- python/damask/_result.py | 10 +++++----- python/damask/mechanics.py | 2 +- python/tests/test_Result.py | 6 +++--- python/tests/test_mechanics.py | 10 +++++----- 4 files changed, 14 insertions(+), 14 deletions(-) diff --git a/python/damask/_result.py b/python/damask/_result.py index 06598e13b..c75968f72 100644 --- a/python/damask/_result.py +++ b/python/damask/_result.py @@ -942,17 +942,17 @@ class Result: @staticmethod - def _add_rotational(F): + def _add_rotation(F): return { - 'data': mechanics.rotational(F['data']).as_matrix(), + 'data': mechanics.rotation(F['data']).as_matrix(), 'label': f"R({F['label']})", 'meta': { 'Unit': F['meta']['Unit'], 'Description': f"Rotational part of {F['label']} ({F['meta']['Description']})", - 'Creator': 'add_rotational' + 'Creator': 'add_rotation' } } - def add_rotational(self,F): + def add_rotation(self,F): """ Add rotational part of a deformation gradient. @@ -962,7 +962,7 @@ class Result: Label of deformation gradient dataset. """ - self._add_generic_pointwise(self._add_rotational,{'F':F}) + self._add_generic_pointwise(self._add_rotation,{'F':F}) @staticmethod diff --git a/python/damask/mechanics.py b/python/damask/mechanics.py index c401a3032..66f5c9916 100644 --- a/python/damask/mechanics.py +++ b/python/damask/mechanics.py @@ -97,7 +97,7 @@ def maximum_shear(T_sym): return (w[...,0] - w[...,2])*0.5 -def rotational(T): +def rotation(T): """ Calculate the rotational part of a tensor. diff --git a/python/tests/test_Result.py b/python/tests/test_Result.py index ce34c48de..07140c20e 100644 --- a/python/tests/test_Result.py +++ b/python/tests/test_Result.py @@ -261,11 +261,11 @@ class TestResult: in_file = default.read_dataset(loc['pole']) assert np.allclose(in_memory,in_file) - def test_add_rotational(self,default): - default.add_rotational('F') + def test_add_rotation(self,default): + default.add_rotation('F') loc = {'F': default.get_dataset_location('F'), 'R(F)': default.get_dataset_location('R(F)')} - in_memory = mechanics.rotational(default.read_dataset(loc['F'],0)).as_matrix() + in_memory = mechanics.rotation(default.read_dataset(loc['F'],0)).as_matrix() in_file = default.read_dataset(loc['R(F)'],0) assert np.allclose(in_memory,in_file) diff --git a/python/tests/test_mechanics.py b/python/tests/test_mechanics.py index b4672dbbf..667747f8e 100644 --- a/python/tests/test_mechanics.py +++ b/python/tests/test_mechanics.py @@ -32,7 +32,7 @@ def stress_second_Piola_Kirchhoff(P,F): return symmetric(S) -def rotational(T): +def rotation(T): return polar_decomposition(T,'R')[0] @@ -109,9 +109,9 @@ class TestMechanics: def test_vectorize_rotation(self): epsilon = Rotation.from_random(self.n).as_matrix() epsilon_vec = np.reshape(epsilon,(self.n//10,10,3,3)) - for i,v in enumerate(np.reshape(mechanics.rotational(epsilon_vec).as_matrix(), - mechanics.rotational(epsilon).as_matrix().shape)): - assert np.allclose(rotational(epsilon[i]),v) + for i,v in enumerate(np.reshape(mechanics.rotation(epsilon_vec).as_matrix(), + mechanics.rotation(epsilon).as_matrix().shape)): + assert np.allclose(rotation(epsilon[i]),v) @pytest.mark.parametrize('vectorized,single',[(mechanics.stress_Cauchy, stress_Cauchy), @@ -146,7 +146,7 @@ class TestMechanics: def test_polar_decomposition(self): """F = RU = VR.""" F = np.broadcast_to(np.eye(3),[self.n,3,3])*np.random.rand(self.n,3,3) - R = mechanics.rotational(F).as_matrix() + R = mechanics.rotation(F).as_matrix() V = mechanics.stretch_left(F) U = mechanics.stretch_right(F) assert np.allclose(np.matmul(R,U), From cca8649d10ae8f94f9d0f1dbf433696667d8c4f9 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Sun, 22 Nov 2020 11:38:39 +0100 Subject: [PATCH 33/35] updated PRIVATE repo --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index 3a4f76fee..5b2a0eb3e 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 3a4f76fee6cbaa8520190fcab1845713adcb1f85 +Subproject commit 5b2a0eb3e5ea0991d54ed9753804f45947cf61b4 From c58f9f23a18a5794d758a0b8f50cd6cd07123438 Mon Sep 17 00:00:00 2001 From: Sharan Roongta Date: Sun, 22 Nov 2020 14:32:32 +0100 Subject: [PATCH 34/35] read old version while resolving merge conflicts --- python/damask/_orientation.py | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index 2fc3c311e..c5e853212 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -854,6 +854,9 @@ class Orientation(Rotation): if vector.shape[-1] != 3: raise ValueError('Input is not a field of three-dimensional vectors.') + vector_ = self.to_SST(vector,proper) if in_SST else \ + self @ np.broadcast_to(vector,self.shape+(3,)) + if self.family == 'cubic': basis = {'improper':np.array([ [-1. , 0. , 1. ], [ np.sqrt(2.) , -np.sqrt(2.) , 0. ], @@ -887,23 +890,23 @@ class Orientation(Rotation): [ 0., 1., 0.] ]), } else: # direct exit for unspecified symmetry - return np.zeros_like(vector) + return np.zeros_like(vector_) if proper: components_proper = np.around(np.einsum('...ji,...i', - np.broadcast_to(basis['proper'], vector.shape+(3,)), - vector), 12) + np.broadcast_to(basis['proper'], vector_.shape+(3,)), + vector_), 12) components_improper = np.around(np.einsum('...ji,...i', - np.broadcast_to(basis['improper'], vector.shape+(3,)), - vector), 12) + np.broadcast_to(basis['improper'], vector_.shape+(3,)), + vector_), 12) in_SST = np.all(components_proper >= 0.0,axis=-1) \ | np.all(components_improper >= 0.0,axis=-1) components = np.where((in_SST & np.all(components_proper >= 0.0,axis=-1))[...,np.newaxis], components_proper,components_improper) else: components = np.around(np.einsum('...ji,...i', - np.broadcast_to(basis['improper'], vector.shape+(3,)), - np.block([vector[...,:2],np.abs(vector[...,2:3])])), 12) + np.broadcast_to(basis['improper'], vector_.shape+(3,)), + np.block([vector_[...,:2],np.abs(vector_[...,2:3])])), 12) in_SST = np.all(components >= 0.0,axis=-1) From 5ab2847b3631fbf9e376068a16c735e8ecb9328a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 23 Nov 2020 20:06:34 +0100 Subject: [PATCH 35/35] polishing --- python/damask/__init__.py | 2 +- python/damask/_orientation.py | 2 +- python/damask/_vtk.py | 2 +- python/tests/conftest.py | 2 +- python/tests/test_Orientation.py | 6 +++--- python/tests/test_Rotation.py | 16 ++++++++-------- 6 files changed, 15 insertions(+), 15 deletions(-) diff --git a/python/damask/__init__.py b/python/damask/__init__.py index 91ca09b29..a1ccfa3f6 100644 --- a/python/damask/__init__.py +++ b/python/damask/__init__.py @@ -15,6 +15,7 @@ with open(_Path(__file__).parent/_Path('VERSION')) as _f: version = _re.sub(r'^v','',_f.readline().strip()) __version__ = version +# make classes directly accessible as damask.Class from ._environment import Environment as _ # noqa environment = _() from . import util # noqa @@ -24,7 +25,6 @@ from . import mechanics # noqa from . import solver # noqa from . import grid_filters # noqa from . import lattice # noqa -# make classes directly accessible as damask.Class from ._rotation import Rotation # noqa from ._orientation import Orientation # noqa from ._table import Table # noqa diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index c5e853212..05301561f 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -851,7 +851,7 @@ class Orientation(Rotation): ... } """ - if vector.shape[-1] != 3: + if np.array(vector).shape[-1] != 3: raise ValueError('Input is not a field of three-dimensional vectors.') vector_ = self.to_SST(vector,proper) if in_SST else \ diff --git a/python/damask/_vtk.py b/python/damask/_vtk.py index 8e6cf490c..2f4d63791 100644 --- a/python/damask/_vtk.py +++ b/python/damask/_vtk.py @@ -140,7 +140,7 @@ class VTK: """ if not os.path.isfile(fname): # vtk has a strange error handling - raise FileNotFoundError(f'No such file: {fname}') + raise FileNotFoundError(f'no such file: {fname}') ext = Path(fname).suffix if ext == '.vtk' or dataset_type is not None: reader = vtk.vtkGenericDataObjectReader() diff --git a/python/tests/conftest.py b/python/tests/conftest.py index 9410ba15a..719fa8e61 100644 --- a/python/tests/conftest.py +++ b/python/tests/conftest.py @@ -82,7 +82,7 @@ def set_of_quaternions(): return qu - n = 1100 + n = 600 scatter=1.e-2 specials = np.array([ [1.0, 0.0, 0.0, 0.0], diff --git a/python/tests/test_Orientation.py b/python/tests/test_Orientation.py index e0f8a1570..edc5a79d8 100644 --- a/python/tests/test_Orientation.py +++ b/python/tests/test_Orientation.py @@ -16,7 +16,7 @@ def reference_dir(reference_dir_base): @pytest.fixture def set_of_rodrigues(set_of_quaternions): - return Rotation(set_of_quaternions).as_Rodrigues_vector()[:200] + return Rotation(set_of_quaternions).as_Rodrigues_vector() class TestOrientation: @@ -275,13 +275,13 @@ class TestOrientation: @pytest.mark.parametrize('lattice',Orientation.crystal_families) def test_in_FZ_vectorization(self,set_of_rodrigues,lattice): - result = Orientation.from_Rodrigues_vector(rho=set_of_rodrigues.reshape((50,4,-1)),lattice=lattice).in_FZ.reshape(-1) + result = Orientation.from_Rodrigues_vector(rho=set_of_rodrigues.reshape((-1,4,4)),lattice=lattice).in_FZ.reshape(-1) for r,rho in zip(result,set_of_rodrigues[:len(result)]): assert r == Orientation.from_Rodrigues_vector(rho=rho,lattice=lattice).in_FZ @pytest.mark.parametrize('lattice',Orientation.crystal_families) def test_in_disorientation_FZ_vectorization(self,set_of_rodrigues,lattice): - result = Orientation.from_Rodrigues_vector(rho=set_of_rodrigues.reshape((50,4,-1)), + result = Orientation.from_Rodrigues_vector(rho=set_of_rodrigues.reshape((-1,4,4)), lattice=lattice).in_disorientation_FZ.reshape(-1) for r,rho in zip(result,set_of_rodrigues[:len(result)]): assert r == Orientation.from_Rodrigues_vector(rho=rho,lattice=lattice).in_disorientation_FZ diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index 8b99f34ec..cc22aba4d 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -862,16 +862,16 @@ class TestRotation: @pytest.mark.parametrize('function,invalid',[(Rotation.from_quaternion, np.array([-1,0,0,0])), - (Rotation.from_quaternion, np.array([1,1,1,0])), - (Rotation.from_Euler_angles, np.array([1,4,0])), - (Rotation.from_axis_angle, np.array([1,0,0,4])), - (Rotation.from_axis_angle, np.array([1,1,0,1])), - (Rotation.from_matrix, np.random.rand(3,3)), - (Rotation.from_matrix, np.array([[1,1,0],[1,2,0],[0,0,1]])), + (Rotation.from_quaternion, np.array([1,1,1,0])), + (Rotation.from_Euler_angles, np.array([1,4,0])), + (Rotation.from_axis_angle, np.array([1,0,0,4])), + (Rotation.from_axis_angle, np.array([1,1,0,1])), + (Rotation.from_matrix, np.random.rand(3,3)), + (Rotation.from_matrix, np.array([[1,1,0],[1,2,0],[0,0,1]])), (Rotation.from_Rodrigues_vector, np.array([1,0,0,-1])), (Rotation.from_Rodrigues_vector, np.array([1,1,0,1])), - (Rotation.from_homochoric, np.array([2,2,2])), - (Rotation.from_cubochoric, np.array([1.1,0,0])) ]) + (Rotation.from_homochoric, np.array([2,2,2])), + (Rotation.from_cubochoric, np.array([1.1,0,0])) ]) def test_invalid_value(self,function,invalid): with pytest.raises(ValueError): function(invalid)