eigenvalues is more specific name than principal components

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
+    w = eigenvalues(x)
-    return (w[2] - w[0])*0.5 if np.shape(x) == (3,3) else \
+    return (w[0] - w[2])*0.5 if np.shape(x) == (3,3) else \
-           (w[:,2] - w[:,0])*0.5
+           (w[:,0] - w[:,2])*0.5
-    
+
-    
+
-def principal_components(x):
+def eigenvalues(x):
     """
-    Return the principal components of a symmetric tensor.
+    Return the eigenvalues, i.e. principal components, of a symmetric tensor.
-    
+
-    The principal components (eigenvalues) are sorted in descending order, each repeated according to
+    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 np.linalg.eigvalsh(symmetric(x))
-    return w[::-1] if np.shape(x) == (3,3) else \
+
-           w[:,::-1]
+
-    
+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)

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],
+         assert np.allclose(mechanics.eigenvalues(x)[self.c],
-                            mechanics.principal_components(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)