diff --git a/python/damask/mechanics.py b/python/damask/mechanics.py index 08f5ca125..307f1d83d 100644 --- a/python/damask/mechanics.py +++ b/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): diff --git a/python/tests/test_mechanics.py b/python/tests/test_mechanics.py index aab92bef3..1ea1f2bba 100644 --- a/python/tests/test_mechanics.py +++ b/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."""