some insights from continuum mechanics formulated as test

python/damask/mechanics.py

@@ -92,21 +92,27 @@ def deviatoric_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):
+def spherical_part(x,tensor=False):
     """
     Return spherical (hydrostatic) part of a tensor.
     
-    A single scalar is returned, i.e. the hydrostatic part is not mapped on the 3rd order identity
-    matrix. 
-      
     Parameters
     ----------
     x : numpy.array of shape (:,3,3) or (3,3)
       Tensor of which the hydrostatic part is computed.
+    tensor : bool, optional
+      Map spherical part onto identity tensor. Default is false
 
     """
-    return np.trace(x)/3.0 if np.shape(x) == (3,3) else \
-           np.trace(x,axis1=1,axis2=2)/3.0
+    if x.shape == (3,3):
+        sph = np.trace(x)/3.0
+        return sph if not tensor else np.eye(3)*sph
+    else:
+        sph = np.trace(x,axis1=1,axis2=2)/3.0
+        if not tensor:
+            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):

python/tests/test_mechanics.py

@@ -30,8 +30,8 @@ class TestMechanics:
 
     def test_vectorize_spherical_part(self):
          x = np.random.random((self.n,3,3))
-         assert np.allclose(mechanics.spherical_part(x)[self.c],
-                            mechanics.spherical_part(x[self.c]))
+         assert np.allclose(mechanics.spherical_part(x,True)[self.c],
+                            mechanics.spherical_part(x[self.c],True))
 
 
     def test_vectorize_Mises_stress(self):
@@ -94,6 +94,15 @@ class TestMechanics:
          assert np.allclose(mechanics.Cauchy(np.broadcast_to(np.eye(3),(self.n,3,3)),P),
                             mechanics.symmetric(P))
 
+    def test_polar_decomposition(self):
+         """F = RU = VR."""
+         F = np.broadcast_to(np.eye(3),[self.n,3,3])*np.random.random((self.n,3,3))
+         R = mechanics.rotational_part(F)
+         V = mechanics.left_stretch(F)
+         U = mechanics.right_stretch(F)
+         assert np.allclose(np.matmul(R,U),
+                            np.matmul(V,R))
+   
 
     def test_strain_tensor_no_rotation(self):
          """Ensure that left and right stretch give same results for no rotation."""
@@ -102,6 +111,12 @@ class TestMechanics:
          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.random.random((self.n,3,3))
+         m = np.random.random()*5.0-2.5
+         assert np.allclose(np.linalg.det(mechanics.strain_tensor(F,'U',m)),
+                            np.linalg.det(mechanics.strain_tensor(F,'V',m)))
 
     def test_strain_tensor_rotation(self):
          """Ensure that pure rotation results in no strain."""
@@ -111,15 +126,46 @@ class TestMechanics:
          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.
+
+         Should be +1, but random F might contain a reflection.
+         """
+         x = np.random.random((self.n,3,3))
+         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.random((self.n,3,3))
-         sph = np.broadcast_to(np.eye(3),(self.n,3,3))\
-             * np.repeat(mechanics.spherical_part(x),9).reshape(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.random((self.n,3,3))
+         full = mechanics.Mises_stress(x)
+         dev  = mechanics.Mises_stress(mechanics.deviatoric_part(x))
+         assert np.allclose(full,
+                            dev)
+
+    def test_spherical_mapping(self):
+         """Ensure that mapping to tensor is correct."""
+         x = np.random.random((self.n,3,3))
+         tensor = mechanics.spherical_part(x,True)
+         scalar = mechanics.spherical_part(x)
+         assert np.allclose(np.linalg.det(tensor),
+                            scalar**3.0)
+
+    def test_spherical_Mises(self):
+         """Ensure that Mises equivalent strrain of spherical strain is 0."""
+         x = np.random.random((self.n,3,3))
+         sph = mechanics.spherical_part(x,True)
+         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."""