separating general tensor math from mechanics operations

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

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

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',

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

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,10 +165,10 @@ 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:
@@ -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.

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
 

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

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)
 

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))
+        eps = 1.0/(2.0*abs(m)) * (+ np.matmul(n,np.einsum('j,kj->jk',w**m,n))
-                                  - np.broadcast_to(np.eye(3),[F_.shape[0],3,3]))
+                                  - 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('j,kj->jk',w**m,n))
-                                  + np.broadcast_to(np.eye(3),[F_.shape[0],3,3]))
+                                  + 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
+    return (T+T.T)*0.5
-
-
-def transpose(T):
-    return T.T
 
 
 def polar_decomposition(T,requested):
@@ -123,10 +112,7 @@ 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',[
+    @pytest.mark.parametrize('vectorized,single',[(mechanics.deviatoric_part, deviatoric_part),
-        (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   ),
@@ -134,8 +120,6 @@ class TestMechanics:
                                                   (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,vectorized,single):
         epsilon     = np.random.rand(self.n,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)

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)