diff --git a/python/damask/__init__.py b/python/damask/__init__.py index 9d0350337..f3d5ad2fb 100644 --- a/python/damask/__init__.py +++ b/python/damask/__init__.py @@ -7,15 +7,16 @@ with open(_Path(__file__).parent/_Path('VERSION')) as _f: version = _re.sub(r'^v','',_f.readline().strip()) __version__ = version -# make classes directly accessible as damask.Class from ._environment import Environment as _ # noqa environment = _() from . import util # noqa from . import seeds # noqa +from . import tensor # noqa from . import mechanics # noqa from . import solver # noqa from . import grid_filters # noqa from . import lattice # noqa +# make classes directly accessible as damask.Class from ._rotation import Rotation # noqa from ._orientation import Orientation # noqa from ._table import Table # noqa diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index e0a00e0a0..9ae1d2ed3 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -2,7 +2,7 @@ import numpy as np from . import Rotation from . import util -from . import mechanics +from . import tensor __parameter_doc__ = \ """lattice : str @@ -362,7 +362,7 @@ class Orientation(Rotation): x = o.to_frame(uvw=uvw) z = o.to_frame(hkl=hkl) om = np.stack([x,np.cross(z,x),z],axis=-2) - return o.copy(rotation=Rotation.from_matrix(mechanics.transpose(om/np.linalg.norm(om,axis=-1,keepdims=True)))) + return o.copy(rotation=Rotation.from_matrix(tensor.transpose(om/np.linalg.norm(om,axis=-1,keepdims=True)))) @property diff --git a/python/damask/_result.py b/python/damask/_result.py index 5bed4347f..a5dcca8db 100644 --- a/python/damask/_result.py +++ b/python/damask/_result.py @@ -18,6 +18,7 @@ from . import Table from . import Orientation from . import grid_filters from . import mechanics +from . import tensor from . import util h5py3 = h5py.__version__[0] == '3' @@ -681,7 +682,7 @@ class Result: label,p = 'Minimum',0 return { - 'data': mechanics.eigenvalues(T_sym['data'])[:,p], + 'data': tensor.eigenvalues(T_sym['data'])[:,p], 'label': f"lambda_{eigenvalue}({T_sym['label']})", 'meta' : { 'Unit': T_sym['meta']['Unit'], @@ -713,7 +714,7 @@ class Result: elif eigenvalue == 'min': label,p = 'minimum',0 return { - 'data': mechanics.eigenvectors(T_sym['data'])[:,p], + 'data': tensor.eigenvectors(T_sym['data'])[:,p], 'label': f"v_{eigenvalue}({T_sym['label']})", 'meta' : { 'Unit': '1', diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index deae8b6ac..e60b14ecb 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -1,6 +1,6 @@ import numpy as np -from . import mechanics +from . import tensor from . import util from . import grid_filters @@ -549,7 +549,7 @@ class Rotation: raise ValueError('Invalid shape.') if reciprocal: - om = np.linalg.inv(mechanics.transpose(om)/np.pi) # transform reciprocal basis set + om = np.linalg.inv(tensor.transpose(om)/np.pi) # transform reciprocal basis set orthonormal = False # contains stretch if not orthonormal: (U,S,Vh) = np.linalg.svd(om) # singular value decomposition diff --git a/python/damask/mechanics.py b/python/damask/mechanics.py index e94c22705..dca26b9de 100644 --- a/python/damask/mechanics.py +++ b/python/damask/mechanics.py @@ -1,5 +1,10 @@ +"""Finite-strain continuum mechanics.""" + +from . import tensor + import numpy as _np + def Cauchy(P,F): """ Return Cauchy stress calculated from first Piola-Kirchhoff stress and deformation gradient. @@ -15,7 +20,7 @@ def Cauchy(P,F): """ sigma = _np.einsum('...,...ij,...kj->...ik',1.0/_np.linalg.det(F),P,F) - return symmetric(sigma) + return tensor.symmetric(sigma) def deviatoric_part(T): @@ -31,43 +36,6 @@ def deviatoric_part(T): return T - _np.einsum('...ij,...->...ij',_np.eye(3),spherical_part(T)) -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. - - Parameters - ---------- - T_sym : numpy.ndarray of shape (...,3,3) - Symmetric tensor of which the eigenvalues are computed. - - """ - 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. - - 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. - - """ - (u,v) = _np.linalg.eigh(symmetric(T_sym)) - - if RHS: - v[_np.linalg.det(v) < 0.0,:,2] *= -1.0 - return v - - def left_stretch(T): """ Return the left stretch of a tensor. @@ -91,7 +59,7 @@ def maximum_shear(T_sym): Symmetric tensor of which the maximum shear is computed. """ - w = eigenvalues(T_sym) + w = tensor.eigenvalues(T_sym) return (w[...,0] - w[...,2])*0.5 @@ -134,7 +102,7 @@ def PK2(P,F): """ S = _np.einsum('...jk,...kl->...jl',_np.linalg.inv(F),P) - return symmetric(S) + return tensor.symmetric(S) def right_stretch(T): @@ -197,16 +165,16 @@ def strain_tensor(F,t,m): """ if t == 'V': - B = _np.matmul(F,transpose(F)) + B = _np.matmul(F,tensor.transpose(F)) w,n = _np.linalg.eigh(B) elif t == 'U': - C = _np.matmul(transpose(F),F) + C = _np.matmul(tensor.transpose(F),F) w,n = _np.linalg.eigh(C) if m > 0.0: eps = 1.0/(2.0*abs(m)) * (+ _np.matmul(n,_np.einsum('...j,...kj->...jk',w**m,n)) - _np.einsum('...jk->...jk',_np.eye(3))) - + elif m < 0.0: eps = 1.0/(2.0*abs(m)) * (- _np.matmul(n,_np.einsum('...j,...kj->...jk',w**m,n)) + _np.einsum('...jk->...jk',_np.eye(3))) @@ -216,32 +184,6 @@ def strain_tensor(F,t,m): return eps -def symmetric(T): - """ - Return the symmetrized tensor. - - Parameters - ---------- - T : numpy.ndarray of shape (...,3,3) - Tensor of which the symmetrized values are computed. - - """ - return (T+transpose(T))*0.5 - - -def transpose(T): - """ - Return the transpose of a tensor. - - Parameters - ---------- - T : numpy.ndarray of shape (...,3,3) - Tensor of which the transpose is computed. - - """ - return _np.swapaxes(T,axis2=-2,axis1=-1) - - def _polar_decomposition(T,requested): """ Singular value decomposition. diff --git a/python/damask/seeds.py b/python/damask/seeds.py index 9aab953d0..c75dfb5e5 100644 --- a/python/damask/seeds.py +++ b/python/damask/seeds.py @@ -1,3 +1,8 @@ +""" +Functionality for generation of seed points for Voronoi +or Laguerre tessellation. +""" + from scipy import spatial as _spatial import numpy as _np diff --git a/python/damask/tensor.py b/python/damask/tensor.py new file mode 100644 index 000000000..71617e32c --- /dev/null +++ b/python/damask/tensor.py @@ -0,0 +1,74 @@ +""" +Tensor operations. + +Notes +----- +This is not a tensor class, but a collection of routines +to operate on numpy.ndarrays of shape (...,3,3). + +""" + +import numpy as _np + + +def symmetric(T): + """ + Return the symmetrized tensor. + + Parameters + ---------- + T : numpy.ndarray of shape (...,3,3) + Tensor of which the symmetrized values are computed. + + """ + return (T+transpose(T))*0.5 + + +def transpose(T): + """ + Return the transpose of a tensor. + + Parameters + ---------- + T : numpy.ndarray of shape (...,3,3) + Tensor of which the transpose is computed. + + """ + 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. + + Parameters + ---------- + T_sym : numpy.ndarray of shape (...,3,3) + Symmetric tensor of which the eigenvalues are computed. + + """ + 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. + + 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. + + """ + (u,v) = _np.linalg.eigh(symmetric(T_sym)) + + if RHS: + v[_np.linalg.det(v) < 0.0,:,2] *= -1.0 + return v diff --git a/python/tests/test_Result.py b/python/tests/test_Result.py index 929c0ef44..b739bb2ca 100644 --- a/python/tests/test_Result.py +++ b/python/tests/test_Result.py @@ -11,6 +11,7 @@ import h5py from damask import Result from damask import Rotation from damask import Orientation +from damask import tensor from damask import mechanics from damask import grid_filters @@ -152,7 +153,7 @@ class TestResult: default.add_eigenvalue('sigma',eigenvalue) loc = {'sigma' :default.get_dataset_location('sigma'), 'lambda':default.get_dataset_location(f'lambda_{eigenvalue}(sigma)')} - in_memory = function(mechanics.eigenvalues(default.read_dataset(loc['sigma'],0)),axis=1,keepdims=True) + in_memory = function(tensor.eigenvalues(default.read_dataset(loc['sigma'],0)),axis=1,keepdims=True) in_file = default.read_dataset(loc['lambda'],0) assert np.allclose(in_memory,in_file) @@ -162,7 +163,7 @@ class TestResult: default.add_eigenvector('sigma',eigenvalue) loc = {'sigma' :default.get_dataset_location('sigma'), 'v(sigma)':default.get_dataset_location(f'v_{eigenvalue}(sigma)')} - in_memory = mechanics.eigenvectors(default.read_dataset(loc['sigma'],0))[:,idx] + in_memory = tensor.eigenvectors(default.read_dataset(loc['sigma'],0))[:,idx] in_file = default.read_dataset(loc['v(sigma)'],0) assert np.allclose(in_memory,in_file) @@ -210,7 +211,7 @@ class TestResult: in_memory = mechanics.Mises_stress(default.read_dataset(loc['sigma'],0)).reshape(-1,1) in_file = default.read_dataset(loc['sigma_vM'],0) assert np.allclose(in_memory,in_file) - + def test_add_Mises_invalid(self,default): default.add_Cauchy('P','F') default.add_calculation('sigma_y','#sigma#',unit='y') diff --git a/python/tests/test_mechanics.py b/python/tests/test_mechanics.py index 9bcf6c6dc..4c0d68572 100644 --- a/python/tests/test_mechanics.py +++ b/python/tests/test_mechanics.py @@ -1,11 +1,12 @@ import pytest import numpy as np +from damask import tensor from damask import mechanics def Cauchy(P,F): sigma = 1.0/np.linalg.det(F) * np.dot(P,F.T) - return mechanics.symmetric(sigma) + return symmetric(sigma) def deviatoric_part(T): @@ -16,14 +17,6 @@ def eigenvalues(T_sym): return np.linalg.eigvalsh(symmetric(T_sym)) -def eigenvectors(T_sym,RHS=False): - (u,v) = np.linalg.eigh(symmetric(T_sym)) - - if RHS: - if np.linalg.det(v) < 0.0: v[:,2] *= -1.0 - return v - - def left_stretch(T): return polar_decomposition(T,'V')[0] @@ -53,39 +46,35 @@ def right_stretch(T): 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_tensor(F,t,m): - F_ = F.reshape(1,3,3) - + if t == 'V': - B = np.matmul(F_,mechanics.transpose(F_)) + B = np.matmul(F,F.T) w,n = np.linalg.eigh(B) elif t == 'U': - C = np.matmul(mechanics.transpose(F_),F_) + C = np.matmul(F.T,F) w,n = np.linalg.eigh(C) if m > 0.0: - eps = 1.0/(2.0*abs(m)) * (+ np.matmul(n,np.einsum('ij,ikj->ijk',w**m,n)) - - np.broadcast_to(np.eye(3),[F_.shape[0],3,3])) + eps = 1.0/(2.0*abs(m)) * (+ np.matmul(n,np.einsum('j,kj->jk',w**m,n)) + - np.eye(3)) elif m < 0.0: - eps = 1.0/(2.0*abs(m)) * (- np.matmul(n,np.einsum('ij,ikj->ijk',w**m,n)) - + np.broadcast_to(np.eye(3),[F_.shape[0],3,3])) + eps = 1.0/(2.0*abs(m)) * (- np.matmul(n,np.einsum('j,kj->jk',w**m,n)) + + np.eye(3)) else: - eps = np.matmul(n,np.einsum('ij,ikj->ijk',0.5*np.log(w),n)) + eps = np.matmul(n,np.einsum('j,kj->jk',0.5*np.log(w),n)) - return eps.reshape(3,3) + return eps def symmetric(T): - return (T+transpose(T))*0.5 - - -def transpose(T): - return T.T + return (T+T.T)*0.5 def polar_decomposition(T,requested): @@ -96,7 +85,7 @@ def polar_decomposition(T,requested): if 'R' in requested: output.append(R) if 'V' in requested: - output.append(np.dot(T,R.T)) + output.append(np.dot(T,R.T)) if 'U' in requested: output.append(np.dot(R.T,T)) @@ -123,20 +112,15 @@ class TestMechanics: 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.eigenvalues , eigenvalues ), - (mechanics.eigenvectors , eigenvectors ), - (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.symmetric , symmetric ), - (mechanics.transpose , transpose ), - ]) + @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 ), + ]) def test_vectorize_1_arg(self,vectorized,single): epsilon = np.random.rand(self.n,3,3) epsilon_vec = np.reshape(epsilon,(self.n//10,10,3,3)) @@ -171,7 +155,7 @@ class TestMechanics: def test_stress_measures(self,function): """Ensure that all stress measures are equivalent for no deformation.""" P = np.random.rand(self.n,3,3) - assert np.allclose(function(P,np.broadcast_to(np.eye(3),(self.n,3,3))),mechanics.symmetric(P)) + 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)) @@ -237,9 +221,9 @@ class TestMechanics: def test_spherical_mapping(self): """Ensure that mapping to tensor is correct.""" x = np.random.rand(self.n,3,3) - tensor = mechanics.spherical_part(x,True) + tnsr = mechanics.spherical_part(x,True) scalar = mechanics.spherical_part(x) - assert np.allclose(np.linalg.det(tensor), + assert np.allclose(np.linalg.det(tnsr), scalar**3.0) def test_spherical_Mises(self): @@ -249,17 +233,6 @@ class TestMechanics: 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.""" - x = np.random.rand(self.n,3,3) - assert np.allclose(mechanics.symmetric(x)*2.0, - mechanics.transpose(x)+x) - - def test_transpose(self): - """Ensure that a symmetric tensor equals its transpose.""" - x = mechanics.symmetric(np.random.rand(self.n,3,3)) - assert np.allclose(mechanics.transpose(x), - x) def test_Mises(self): """Ensure that equivalent stress is 3/2 of equivalent strain.""" @@ -267,31 +240,7 @@ class TestMechanics: assert np.allclose(mechanics.Mises_stress(x)/mechanics.Mises_strain(x), 1.5) - def test_eigenvalues(self): - """Ensure that the characteristic polynomial can be solved.""" - A = mechanics.symmetric(np.random.rand(self.n,3,3)) - lambd = mechanics.eigenvalues(A) - s = np.random.randint(self.n) - for i in range(3): - assert np.allclose(np.linalg.det(A[s]-lambd[s,i]*np.eye(3)),.0) - - def test_eigenvalues_and_vectors(self): - """Ensure that eigenvalues and -vectors are the solution to the characteristic polynomial.""" - A = mechanics.symmetric(np.random.rand(self.n,3,3)) - lambd = mechanics.eigenvalues(A) - x = mechanics.eigenvectors(A) - s = np.random.randint(self.n) - for i in range(3): - assert np.allclose(np.dot(A[s]-lambd[s,i]*np.eye(3),x[s,:,i]),.0) - - def test_eigenvectors_RHS(self): - """Ensure that RHS coordinate system does only change sign of determinant.""" - A = mechanics.symmetric(np.random.rand(self.n,3,3)) - LRHS = np.linalg.det(mechanics.eigenvectors(A,RHS=False)) - RHS = np.linalg.det(mechanics.eigenvectors(A,RHS=True)) - assert np.allclose(np.abs(LRHS),RHS) - def test_spherical_no_shear(self): """Ensure that sherical stress has max shear of 0.0.""" - A = mechanics.spherical_part(mechanics.symmetric(np.random.rand(self.n,3,3)),True) + A = mechanics.spherical_part(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 new file mode 100644 index 000000000..24fbe6126 --- /dev/null +++ b/python/tests/test_tensor.py @@ -0,0 +1,73 @@ +import pytest +import numpy as np + +from damask import tensor + + +def eigenvalues(T_sym): + return np.linalg.eigvalsh(symmetric(T_sym)) + +def eigenvectors(T_sym,RHS=False): + (u,v) = np.linalg.eigh(symmetric(T_sym)) + + if RHS: + if np.linalg.det(v) < 0.0: v[:,2] *= -1.0 + return v + +def symmetric(T): + return (T+transpose(T))*0.5 + +def transpose(T): + return T.T + + +class TestTensor: + + n = 1000 + c = np.random.randint(n) + + + @pytest.mark.parametrize('vectorized,single',[(tensor.eigenvalues , eigenvalues ), + (tensor.eigenvectors , eigenvectors ), + (tensor.symmetric , symmetric ), + (tensor.transpose , transpose ), + ]) + def test_vectorize_1_arg(self,vectorized,single): + epsilon = np.random.rand(self.n,3,3) + epsilon_vec = np.reshape(epsilon,(self.n//10,10,3,3)) + for i,v in enumerate(np.reshape(vectorized(epsilon_vec),vectorized(epsilon).shape)): + assert np.allclose(single(epsilon[i]),v) + + def test_symmetric(self): + """Ensure that a symmetric tensor is half of the sum of a tensor and its transpose.""" + x = np.random.rand(self.n,3,3) + assert np.allclose(tensor.symmetric(x)*2.0,tensor.transpose(x)+x) + + def test_transpose(self): + """Ensure that a symmetric tensor equals its transpose.""" + x = tensor.symmetric(np.random.rand(self.n,3,3)) + assert np.allclose(tensor.transpose(x),x) + + def test_eigenvalues(self): + """Ensure that the characteristic polynomial can be solved.""" + A = tensor.symmetric(np.random.rand(self.n,3,3)) + lambd = tensor.eigenvalues(A) + s = np.random.randint(self.n) + for i in range(3): + assert np.allclose(np.linalg.det(A[s]-lambd[s,i]*np.eye(3)),.0) + + def test_eigenvalues_and_vectors(self): + """Ensure that eigenvalues and -vectors are the solution to the characteristic polynomial.""" + A = tensor.symmetric(np.random.rand(self.n,3,3)) + lambd = tensor.eigenvalues(A) + x = tensor.eigenvectors(A) + s = np.random.randint(self.n) + for i in range(3): + assert np.allclose(np.dot(A[s]-lambd[s,i]*np.eye(3),x[s,:,i]),.0) + + def test_eigenvectors_RHS(self): + """Ensure that RHS coordinate system does only change sign of determinant.""" + A = tensor.symmetric(np.random.rand(self.n,3,3)) + 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)