From 6e5cb6013252ddebf037f6e7eb42cc4125e874b4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 19 Nov 2020 14:38:54 +0100 Subject: [PATCH] general tensor functionality --- python/damask/__init__.py | 10 +++++- python/damask/_result.py | 4 +-- python/damask/mechanics.py | 53 ++++------------------------- python/damask/tensor.py | 41 +++++++++++++++++++++++ python/tests/test_Result.py | 4 +-- python/tests/test_mechanics.py | 61 ++++++---------------------------- python/tests/test_tensor.py | 34 ++++++++++++++++--- 7 files changed, 101 insertions(+), 106 deletions(-) diff --git a/python/damask/__init__.py b/python/damask/__init__.py index f3d5ad2fb..91ca09b29 100644 --- a/python/damask/__init__.py +++ b/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 diff --git a/python/damask/_result.py b/python/damask/_result.py index cae3298f7..8bacebb90 100644 --- a/python/damask/_result.py +++ b/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'], diff --git a/python/damask/mechanics.py b/python/damask/mechanics.py index 1531c66d3..2be93b91b 100644 --- a/python/damask/mechanics.py +++ b/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))) diff --git a/python/damask/tensor.py b/python/damask/tensor.py index 4b22560b6..3ee551d02 100644 --- a/python/damask/tensor.py +++ b/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. diff --git a/python/tests/test_Result.py b/python/tests/test_Result.py index 46b947cdc..fb79f4db7 100644 --- a/python/tests/test_Result.py +++ b/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) diff --git a/python/tests/test_mechanics.py b/python/tests/test_mechanics.py index 4a61f5b1e..4912aa962 100644 --- a/python/tests/test_mechanics.py +++ b/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.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), + @pytest.mark.parametrize('vectorized,single',[(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): - """Ensure that mapping to tensor is correct.""" + @pytest.mark.parametrize('Mises_equivalent',[mechanics.equivalent_strain_Mises, + 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) - 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), + assert np.allclose(Mises_equivalent(tensor.spherical(x,True)), 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) diff --git a/python/tests/test_tensor.py b/python/tests/test_tensor.py index 24fbe6126..58194130e 100644 --- a/python/tests/test_tensor.py +++ b/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 ), - (tensor.eigenvectors , eigenvectors ), - (tensor.symmetric , symmetric ), - (tensor.transpose , transpose ), + @pytest.mark.parametrize('vectorized,single',[(tensor.deviatoric, deviatoric), + (tensor.eigenvalues, eigenvalues), + (tensor.eigenvectors, eigenvectors), + (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)