From b893967b68a828e5f8fee8475cd2ffab06fdb965 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 16 Nov 2020 01:01:32 +0100 Subject: [PATCH] more systematic names and extended docstrings --- python/damask/_result.py | 2 +- python/damask/mechanics.py | 186 +++++++++++++++++++++++++-------- python/damask/tensor.py | 37 +++++-- python/tests/test_Result.py | 4 +- python/tests/test_mechanics.py | 24 ++--- 5 files changed, 182 insertions(+), 71 deletions(-) diff --git a/python/damask/_result.py b/python/damask/_result.py index a5dcca8db..ad8c5ec0e 100644 --- a/python/damask/_result.py +++ b/python/damask/_result.py @@ -1024,7 +1024,7 @@ class Result: @staticmethod def _add_stretch_tensor(F,t): return { - 'data': (mechanics.left_stretch if t.upper() == 'V' else mechanics.right_stretch)(F['data']), + 'data': (mechanics.stretch_left if t.upper() == 'V' else mechanics.stretch_right)(F['data']), 'label': f"{t}({F['label']})", 'meta': { 'Unit': F['meta']['Unit'], diff --git a/python/damask/mechanics.py b/python/damask/mechanics.py index dca26b9de..9118ea875 100644 --- a/python/damask/mechanics.py +++ b/python/damask/mechanics.py @@ -5,18 +5,58 @@ from . import tensor import numpy as _np -def Cauchy(P,F): +def Cauchy_Green_deformation_left(F): """ - Return Cauchy stress calculated from first Piola-Kirchhoff stress and deformation gradient. - - Resulting tensor is symmetrized as the Cauchy stress needs to be symmetric. + Calculate left Cauchy-Green deformation tensor (Finger deformation tensor). Parameters ---------- F : numpy.ndarray of shape (...,3,3) Deformation gradient. + + Returns + ------- + B : numpy.ndarray of shape (...,3,3) + Left Cauchy-Green deformation tensor. + + """ + return _np.matmul(F,tensor.transpose(F)) + + +def Cauchy_Green_deformation_right(F): + """ + Calculate right Cauchy-Green deformation tensor. + + Parameters + ---------- + F : numpy.ndarray of shape (...,3,3) + Deformation gradient. + + Returns + ------- + C : numpy.ndarray of shape (...,3,3) + Right Cauchy-Green deformation tensor. + + """ + return _np.matmul(tensor.transpose(F),F) + +def Cauchy(P,F): + """ + Calculate the Cauchy (true) stress. + + Resulting tensor is symmetrized as the Cauchy stress needs to be symmetric. + + Parameters + ---------- P : numpy.ndarray of shape (...,3,3) First Piola-Kirchhoff stress. + F : numpy.ndarray of shape (...,3,3) + Deformation gradient. + + Returns + ------- + sigma : numpy.ndarray of shape (...,3,3) + Cauchy stress. """ sigma = _np.einsum('...,...ij,...kj->...ik',1.0/_np.linalg.det(F),P,F) @@ -25,39 +65,36 @@ def Cauchy(P,F): def deviatoric_part(T): """ - Return deviatoric part of a tensor. + 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 - _np.einsum('...ij,...->...ij',_np.eye(3),spherical_part(T)) -def left_stretch(T): - """ - Return the left stretch of a tensor. - - Parameters - ---------- - T : numpy.ndarray of shape (...,3,3) - Tensor of which the left stretch is computed. - - """ - return _polar_decomposition(T,'V')[0] - - def maximum_shear(T_sym): """ - Return the maximum shear component of a symmetric tensor. + Calculate the maximum shear component of a symmetric tensor. Parameters ---------- T_sym : numpy.ndarray of shape (...,3,3) Symmetric tensor of which the maximum shear is computed. + Returns + ------- + gamma_max : numpy.ndarray of shape (...) + Maximum shear of T_sym. + """ w = tensor.eigenvalues(T_sym) return (w[...,0] - w[...,2])*0.5 @@ -65,33 +102,46 @@ def maximum_shear(T_sym): def Mises_strain(epsilon): """ - Return the Mises equivalent of a strain tensor. + Calculate the Mises equivalent of a strain tensor. Parameters ---------- epsilon : numpy.ndarray of shape (...,3,3) Symmetric strain tensor of which the von Mises equivalent is computed. + Returns + ------- + epsilon_vM : numpy.ndarray of shape (...) + Von Mises equivalent strain of epsilon. + """ return _Mises(epsilon,2.0/3.0) def Mises_stress(sigma): """ - Return the Mises equivalent of a stress tensor. + Calculate the Mises equivalent of a stress tensor. Parameters ---------- sigma : numpy.ndarray of shape (...,3,3) Symmetric stress tensor of which the von Mises equivalent is computed. + Returns + ------- + sigma_vM : numpy.ndarray of shape (...) + Von Mises equivalent stress of sigma. + """ return _Mises(sigma,3.0/2.0) def PK2(P,F): """ - Calculate second Piola-Kirchhoff stress from first Piola-Kirchhoff stress and deformation gradient. + Calculate the second Piola-Kirchhoff stress. + + Resulting tensor is symmetrized as the second Piola-Kirchhoff stress + needs to be symmetric. Parameters ---------- @@ -100,47 +150,51 @@ def PK2(P,F): F : numpy.ndarray of shape (...,3,3) Deformation gradient. + Returns + ------- + S : numpy.ndarray of shape (...,3,3) + Second Piola-Kirchhoff stress. + """ S = _np.einsum('...jk,...kl->...jl',_np.linalg.inv(F),P) return tensor.symmetric(S) -def right_stretch(T): - """ - Return the right stretch of a tensor. - - Parameters - ---------- - T : numpy.ndarray of shape (...,3,3) - Tensor of which the right stretch is computed. - - """ - return _polar_decomposition(T,'U')[0] - - def rotational_part(T): """ - Return the rotational part of a tensor. + Calculate the rotational part of a tensor. Parameters ---------- T : numpy.ndarray of shape (...,3,3) Tensor of which the rotational part is computed. + Returns + ------- + R : numpy.ndarray of shape (...,3,3) + Rotational part. + """ return _polar_decomposition(T,'R')[0] def spherical_part(T,tensor=False): """ - Return spherical (hydrostatic) part of a tensor. + 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. Default is false + 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 @@ -149,7 +203,7 @@ def spherical_part(T,tensor=False): def strain_tensor(F,t,m): """ - Return strain tensor calculated from deformation gradient. + Calculate strain tensor from deformation gradient. For details refer to https://en.wikipedia.org/wiki/Finite_strain_theory and https://de.wikipedia.org/wiki/Verzerrungstensor @@ -159,17 +213,21 @@ def strain_tensor(F,t,m): F : numpy.ndarray of shape (...,3,3) Deformation gradient. t : {‘V’, ‘U’} - Type of the polar decomposition, ‘V’ for left stretch tensor and ‘U’ for right stretch tensor. + Type of the polar decomposition, ‘V’ for left stretch tensor + and ‘U’ for right stretch tensor. m : float Order of the strain. + Returns + ------- + epsilon : numpy.ndarray of shape (...,3,3) + Strain of F. + """ if t == 'V': - B = _np.matmul(F,tensor.transpose(F)) - w,n = _np.linalg.eigh(B) + w,n = _np.linalg.eigh(Cauchy_Green_deformation_left(F)) elif t == 'U': - C = _np.matmul(tensor.transpose(F),F) - w,n = _np.linalg.eigh(C) + w,n = _np.linalg.eigh(Cauchy_Green_deformation_right(F)) if m > 0.0: eps = 1.0/(2.0*abs(m)) * (+ _np.matmul(n,_np.einsum('...j,...kj->...jk',w**m,n)) @@ -184,9 +242,45 @@ def strain_tensor(F,t,m): return eps +def stretch_left(T): + """ + Calculate left stretch of a tensor. + + Parameters + ---------- + T : numpy.ndarray of shape (...,3,3) + Tensor of which the left stretch is computed. + + Returns + ------- + V : numpy.ndarray of shape (...,3,3) + Left stretch tensor from Polar decomposition of T. + + """ + return _polar_decomposition(T,'V')[0] + + +def stretch_right(T): + """ + Calculate right stretch of a tensor. + + Parameters + ---------- + T : numpy.ndarray of shape (...,3,3) + Tensor of which the right stretch is computed. + + Returns + ------- + U : numpy.ndarray of shape (...,3,3) + Left stretch tensor from Polar decomposition of T. + + """ + return _polar_decomposition(T,'U')[0] + + def _polar_decomposition(T,requested): """ - Singular value decomposition. + Perform singular value decomposition. Parameters ---------- @@ -197,7 +291,7 @@ def _polar_decomposition(T,requested): ‘V’ for left stretch tensor and ‘U’ for right stretch tensor. """ - u, s, vh = _np.linalg.svd(T) + u, _, vh = _np.linalg.svd(T) R = _np.einsum('...ij,...jk->...ik',u,vh) output = [] diff --git a/python/damask/tensor.py b/python/damask/tensor.py index 71617e32c..6e5bdaf85 100644 --- a/python/damask/tensor.py +++ b/python/damask/tensor.py @@ -13,58 +13,75 @@ import numpy as _np def symmetric(T): """ - Return the symmetrized tensor. + Symmetrize tensor. Parameters ---------- T : numpy.ndarray of shape (...,3,3) Tensor of which the symmetrized values are computed. + Returns + ------- + T_sym : numpy.ndarray of shape (...,3,3) + Symmetrized tensor T. + """ return (T+transpose(T))*0.5 def transpose(T): """ - Return the transpose of a tensor. + Transpose tensor. Parameters ---------- T : numpy.ndarray of shape (...,3,3) Tensor of which the transpose is computed. + Returns + ------- + T.T : numpy.ndarray of shape (...,3,3) + Transpose of tensor T. + """ 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. + Eigenvalues, i.e. principal components, of a symmetric tensor. Parameters ---------- T_sym : numpy.ndarray of shape (...,3,3) Symmetric tensor of which the eigenvalues are computed. + Returns + ------- + lambda : numpy.ndarray of shape (...,3) + Eigenvalues of T_sym sorted in ascending order, each repeated + according to its multiplicity. + """ 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. + Eigenvectors of a symmetric tensor. 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. + Enforce right-handed coordinate system. Defaults to False. + + Returns + ------- + x : numpy.ndarray of shape (...,3,3) + Eigenvectors of T_sym sorted in ascending order of their + associated eigenvalues. """ (u,v) = _np.linalg.eigh(symmetric(T_sym)) diff --git a/python/tests/test_Result.py b/python/tests/test_Result.py index b739bb2ca..89379ce53 100644 --- a/python/tests/test_Result.py +++ b/python/tests/test_Result.py @@ -292,7 +292,7 @@ class TestResult: default.add_stretch_tensor('F','U') loc = {'F': default.get_dataset_location('F'), 'U(F)': default.get_dataset_location('U(F)')} - in_memory = mechanics.right_stretch(default.read_dataset(loc['F'],0)) + in_memory = mechanics.stretch_right(default.read_dataset(loc['F'],0)) in_file = default.read_dataset(loc['U(F)'],0) assert np.allclose(in_memory,in_file) @@ -300,7 +300,7 @@ class TestResult: default.add_stretch_tensor('F','V') loc = {'F': default.get_dataset_location('F'), 'V(F)': default.get_dataset_location('V(F)')} - in_memory = mechanics.left_stretch(default.read_dataset(loc['F'],0)) + in_memory = mechanics.stretch_left(default.read_dataset(loc['F'],0)) in_file = default.read_dataset(loc['V(F)'],0) assert np.allclose(in_memory,in_file) diff --git a/python/tests/test_mechanics.py b/python/tests/test_mechanics.py index 4c0d68572..8af8d1a03 100644 --- a/python/tests/test_mechanics.py +++ b/python/tests/test_mechanics.py @@ -17,10 +17,6 @@ def eigenvalues(T_sym): return np.linalg.eigvalsh(symmetric(T_sym)) -def left_stretch(T): - return polar_decomposition(T,'V')[0] - - def maximum_shear(T_sym): w = eigenvalues(T_sym) return (w[0] - w[2])*0.5 @@ -39,10 +35,6 @@ def PK2(P,F): return symmetric(S) -def right_stretch(T): - return polar_decomposition(T,'U')[0] - - def rotational_part(T): return polar_decomposition(T,'R')[0] @@ -73,6 +65,14 @@ def strain_tensor(F,t,m): return eps +def stretch_left(T): + return polar_decomposition(T,'V')[0] + + +def stretch_right(T): + return polar_decomposition(T,'U')[0] + + def symmetric(T): return (T+T.T)*0.5 @@ -113,13 +113,13 @@ class TestMechanics: assert np.allclose(single(test_data_flat[i]),v) @pytest.mark.parametrize('vectorized,single',[(mechanics.deviatoric_part, deviatoric_part), - (mechanics.left_stretch , left_stretch ), (mechanics.maximum_shear , maximum_shear ), (mechanics.Mises_strain , Mises_strain ), (mechanics.Mises_stress , Mises_stress ), - (mechanics.right_stretch , right_stretch ), (mechanics.rotational_part, rotational_part), (mechanics.spherical_part , spherical_part ), + (mechanics.stretch_left , stretch_left ), + (mechanics.stretch_right , stretch_right ), ]) def test_vectorize_1_arg(self,vectorized,single): epsilon = np.random.rand(self.n,3,3) @@ -166,8 +166,8 @@ class TestMechanics: """F = RU = VR.""" F = np.broadcast_to(np.eye(3),[self.n,3,3])*np.random.rand(self.n,3,3) R = mechanics.rotational_part(F) - V = mechanics.left_stretch(F) - U = mechanics.right_stretch(F) + V = mechanics.stretch_left(F) + U = mechanics.stretch_right(F) assert np.allclose(np.matmul(R,U), np.matmul(V,R))