diff --git a/python/damask/mechanics.py b/python/damask/mechanics.py index 307f1d83d..bf23edb40 100644 --- a/python/damask/mechanics.py +++ b/python/damask/mechanics.py @@ -3,9 +3,9 @@ import numpy as np def Cauchy(F,P): """ Return Cauchy stress calculated from 1. Piola-Kirchhoff stress and deformation gradient. - + Resulting tensor is symmetrized as the Cauchy stress needs to be symmetric. - + Parameters ---------- F : numpy.array of shape (:,3,3) or (3,3) @@ -24,7 +24,7 @@ def Cauchy(F,P): def PK2(F,P): """ Return 2. Piola-Kirchhoff stress calculated from 1. Piola-Kirchhoff stress and deformation gradient. - + Parameters ---------- F : numpy.array of shape (:,3,3) or (3,3) @@ -37,16 +37,16 @@ def PK2(F,P): S = np.dot(np.linalg.inv(F),P) else: S = np.einsum('ijk,ikl->ijl',np.linalg.inv(F),P) - return S + return symmetric(S) def strain_tensor(F,t,m): """ Return strain tensor calculated from deformation gradient. - + For details refer to https://en.wikipedia.org/wiki/Finite_strain_theory and https://de.wikipedia.org/wiki/Verzerrungstensor - + Parameters ---------- F : numpy.array of shape (:,3,3) or (3,3) @@ -64,16 +64,16 @@ def strain_tensor(F,t,m): 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)) + 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])) 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])) else: eps = np.matmul(n,np.einsum('ij,ikj->ijk',0.5*np.log(w),n)) - + return eps.reshape((3,3)) if np.shape(F) == (3,3) else \ eps @@ -81,7 +81,7 @@ def strain_tensor(F,t,m): def deviatoric_part(x): """ Return deviatoric part of a tensor. - + Parameters ---------- x : numpy.array of shape (:,3,3) or (3,3) @@ -89,13 +89,13 @@ def deviatoric_part(x): """ return x - np.eye(3)*spherical_part(x) if np.shape(x) == (3,3) else \ - x - np.einsum('ijk,i->ijk',np.broadcast_to(np.eye(3),[x.shape[0],3,3]),spherical_part(x)) + x - np.einsum('ijk,i->ijk',np.broadcast_to(np.eye(3),[x.shape[0],3,3]),spherical_part(x)) def spherical_part(x,tensor=False): """ Return spherical (hydrostatic) part of a tensor. - + Parameters ---------- x : numpy.array of shape (:,3,3) or (3,3) @@ -113,12 +113,12 @@ def spherical_part(x,tensor=False): return sph else: return np.einsum('ijk,i->ijk',np.broadcast_to(np.eye(3),(x.shape[0],3,3)),sph) - - + + def Mises_stress(sigma): """ Return the Mises equivalent of a stress tensor. - + Parameters ---------- sigma : numpy.array of shape (:,3,3) or (3,3) @@ -128,12 +128,12 @@ def Mises_stress(sigma): s = deviatoric_part(sigma) return np.sqrt(3.0/2.0*(np.sum(s**2.0))) if np.shape(sigma) == (3,3) else \ np.sqrt(3.0/2.0*np.einsum('ijk->i',s**2.0)) - - + + def Mises_strain(epsilon): """ Return the Mises equivalent of a strain tensor. - + Parameters ---------- epsilon : numpy.array of shape (:,3,3) or (3,3) @@ -148,7 +148,7 @@ def Mises_strain(epsilon): def symmetric(x): """ Return the symmetrized tensor. - + Parameters ---------- x : numpy.array of shape (:,3,3) or (3,3) @@ -161,40 +161,54 @@ def symmetric(x): def maximum_shear(x): """ Return the maximum shear component of a symmetric tensor. - + Parameters ---------- x : numpy.array of shape (:,3,3) or (3,3) Symmetric tensor of which the maximum shear is computed. """ - w = np.linalg.eigvalsh(symmetric(x)) # eigenvalues in ascending order - return (w[2] - w[0])*0.5 if np.shape(x) == (3,3) else \ - (w[:,2] - w[:,0])*0.5 - - -def principal_components(x): + w = eigenvalues(x) + return (w[0] - w[2])*0.5 if np.shape(x) == (3,3) else \ + (w[:,0] - w[:,2])*0.5 + + +def eigenvalues(x): """ - Return the principal components of a symmetric tensor. - - The principal components (eigenvalues) are sorted in descending order, each repeated according to + 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 ---------- x : numpy.array of shape (:,3,3) or (3,3) - Symmetric tensor of which the principal compontents are computed. + Symmetric tensor of which the eigenvalues are computed. """ - w = np.linalg.eigvalsh(symmetric(x)) # eigenvalues in ascending order - return w[::-1] if np.shape(x) == (3,3) else \ - w[:,::-1] - - + return np.linalg.eigvalsh(symmetric(x)) + + +def eigenvectors(x): + """ + Return eigenvectors of a symmetric tensor. + + The eigenvalues are sorted in ascending order of their associated eigenvalues. + + Parameters + ---------- + x : numpy.array of shape (:,3,3) or (3,3) + Symmetric tensor of which the eigenvectors are computed. + + """ + (u,v) = np.linalg.eigh(symmetric(x)) + return v + + def transpose(x): """ Return the transpose of a tensor. - + Parameters ---------- x : numpy.array of shape (:,3,3) or (3,3) @@ -208,7 +222,7 @@ def transpose(x): def rotational_part(x): """ Return the rotational part of a tensor. - + Parameters ---------- x : numpy.array of shape (:,3,3) or (3,3) @@ -221,7 +235,7 @@ def rotational_part(x): def left_stretch(x): """ Return the left stretch of a tensor. - + Parameters ---------- x : numpy.array of shape (:,3,3) or (3,3) @@ -229,12 +243,12 @@ def left_stretch(x): """ return __polar_decomposition(x,'V')[0] - - + + def right_stretch(x): """ Return the right stretch of a tensor. - + Parameters ---------- x : numpy.array of shape (:,3,3) or (3,3) @@ -247,20 +261,20 @@ def right_stretch(x): def __polar_decomposition(x,requested): """ Singular value decomposition. - + Parameters ---------- x : numpy.array of shape (:,3,3) or (3,3) Tensor of which the singular values are computed. requested : iterable of str - Requested outputs: ‘R’ for the rotation tensor, + Requested outputs: ‘R’ for the rotation tensor, ‘V’ for left stretch tensor and ‘U’ for right stretch tensor. """ u, s, vh = np.linalg.svd(x) R = np.dot(u,vh) if np.shape(x) == (3,3) else \ np.einsum('ijk,ikl->ijl',u,vh) - + output = [] if 'R' in requested: output.append(R) @@ -268,5 +282,5 @@ def __polar_decomposition(x,requested): output.append(np.dot(x,R.T) if np.shape(x) == (3,3) else np.einsum('ijk,ilk->ijl',x,R)) if 'U' in requested: output.append(np.dot(R.T,x) if np.shape(x) == (3,3) else np.einsum('ikj,ikl->ijl',R,x)) - + return tuple(output) diff --git a/python/tests/test_mechanics.py b/python/tests/test_mechanics.py index 9e1d9bc0c..483452a12 100644 --- a/python/tests/test_mechanics.py +++ b/python/tests/test_mechanics.py @@ -2,10 +2,10 @@ import numpy as np from damask import mechanics class TestMechanics: - + n = 1000 c = np.random.randint(n) - + def test_vectorize_Cauchy(self): P = np.random.random((self.n,3,3)) @@ -58,10 +58,23 @@ class TestMechanics: mechanics.maximum_shear(x[self.c])) - def test_vectorize_principal_components(self): + def test_vectorize_eigenvalues(self): x = np.random.random((self.n,3,3)) - assert np.allclose(mechanics.principal_components(x)[self.c], - mechanics.principal_components(x[self.c])) + assert np.allclose(mechanics.eigenvalues(x)[self.c], + mechanics.eigenvalues(x[self.c])) + + + def test_vectorize_eigenvectors(self): + x = np.random.random((self.n,3,3)) + assert np.allclose(mechanics.eigenvectors(x)[self.c], + mechanics.eigenvectors(x[self.c])) + + + def test_vectorize_PK2(self): + F = np.random.random((self.n,3,3)) + P = np.random.random((self.n,3,3)) + assert np.allclose(mechanics.PK2(F,P)[self.c], + mechanics.PK2(F[self.c],P[self.c])) def test_vectorize_transpose(self): @@ -102,7 +115,14 @@ class TestMechanics: U = mechanics.right_stretch(F) assert np.allclose(np.matmul(R,U), np.matmul(V,R)) - + + + def test_PK2(self): + """Ensure 2. Piola-Kirchhoff stress is symmetrized 1. Piola-Kirchhoff stress for no deformation.""" + P = np.random.random((self.n,3,3)) + assert np.allclose(mechanics.PK2(np.broadcast_to(np.eye(3),(self.n,3,3)),P), + mechanics.symmetric(P)) + def test_strain_tensor_no_rotation(self): """Ensure that left and right stretch give same results for no rotation.""" @@ -110,7 +130,7 @@ class TestMechanics: m = np.random.random()*20.0-10.0 assert np.allclose(mechanics.strain_tensor(F,'U',m), mechanics.strain_tensor(F,'V',m)) - + def test_strain_tensor_rotation_equivalence(self): """Ensure that left and right strain differ only by a rotation.""" F = np.broadcast_to(np.eye(3),[self.n,3,3]) + (np.random.random((self.n,3,3))*0.5 - 0.25) @@ -125,7 +145,7 @@ class TestMechanics: m = np.random.random()*2.0 - 1.0 assert np.allclose(mechanics.strain_tensor(F,t,m), 0.0) - + def test_rotation_determinant(self): """ Ensure that the determinant of the rotational part is +- 1. @@ -186,3 +206,22 @@ class TestMechanics: x = np.random.random((self.n,3,3)) 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.random((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.random((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)