general tensor functionality

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
 

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

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)))

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.

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)
 

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.maximum_shear,           maximum_shear),
-    @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),
                                                   (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):
+    @pytest.mark.parametrize('Mises_equivalent',[mechanics.equivalent_strain_Mises,
-        """Ensure that mapping to tensor is correct."""
+                                                 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)
+        assert np.allclose(Mises_equivalent(tensor.spherical(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),
                            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)

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    ),
+    @pytest.mark.parametrize('vectorized,single',[(tensor.deviatoric,   deviatoric),
-                                                  (tensor.eigenvectors   , eigenvectors   ),
+                                                  (tensor.eigenvalues,  eigenvalues),
-                                                  (tensor.symmetric      , symmetric      ),
+                                                  (tensor.eigenvectors, eigenvectors),
-                                                  (tensor.transpose      , transpose      ),
+                                                  (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)