223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

vectorize Rotation.fromXXX functions

This commit is contained in:
Martin Diehl 2020-04-15 23:00:00 +02:00
parent 3d10266fbc
commit 23fc58699f
4 changed files with 114 additions and 86 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
python/.coveragerc
View File

@ -1,2 +1,5 @@
[run] [run]
omit = tests/* omit = tests/*
damask/_asciitable.py
damask/_test.py
damask/config/*

172
python/damask/_rotation.py
View File

@ -1,6 +1,7 @@
import numpy as np import numpy as np
from ._Lambert import ball_to_cube, cube_to_ball from ._Lambert import ball_to_cube, cube_to_ball
from . import mechanics
_P = -1 _P = -1
@ -61,6 +62,8 @@ class Rotation:
def __repr__(self): def __repr__(self):
"""Orientation displayed as unit quaternion, rotation matrix, and Bunge-Euler angles.""" """Orientation displayed as unit quaternion, rotation matrix, and Bunge-Euler angles."""
if self.quaternion.shape != (4,):
raise NotImplementedError
return '\n'.join([ return '\n'.join([
'Quaternion: (real={:.3f}, imag=<{:+.3f}, {:+.3f}, {:+.3f}>)'.format(*(self.quaternion)), 'Quaternion: (real={:.3f}, imag=<{:+.3f}, {:+.3f}, {:+.3f}>)'.format(*(self.quaternion)),
'Matrix:\n{}'.format(self.asMatrix()), 'Matrix:\n{}'.format(self.asMatrix()),
@ -116,7 +119,7 @@ class Rotation:
def inverse(self): def inverse(self):
"""In-place inverse rotation/backward rotation.""" """In-place inverse rotation/backward rotation."""
self.quaternion[1:] *= -1 self.quaternion[...,1:] *= -1
return self return self
def inversed(self): def inversed(self):
@ -125,12 +128,12 @@ class Rotation:
def standardize(self): def standardize(self):
"""In-place quaternion representation with positive q.""" """In-place quaternion representation with positive real part."""
if self.quaternion[0] < 0.0: self.quaternion*=-1 self.quaternion[self.quaternion[...,0] < 0.0] *= -1
return self return self
def standardized(self): def standardized(self):
"""Quaternion representation with positive q.""" """Quaternion representation with positive real."""
return self.copy().standardize() return self.copy().standardize()
@ -165,7 +168,7 @@ class Rotation:
def asQuaternion(self): def asQuaternion(self):
""" """
Unit quaternion [q, p_1, p_2, p_3] unless quaternion == True: damask.quaternion object. Unit quaternion [q, p_1, p_2, p_3].
Parameters Parameters
---------- ----------
@ -251,107 +254,106 @@ class Rotation:
# static constructors. The input data needs to follow the convention, options allow to # static constructors. The input data needs to follow the convention, options allow to
# relax these convections # relax these convections
@staticmethod @staticmethod
def fromQuaternion(quaternion, def from_quaternion(quaternion,
acceptHomomorph = False, acceptHomomorph = False,
P = -1): P = -1):
qu = quaternion if isinstance(quaternion,np.ndarray) and quaternion.dtype == np.dtype(float) \ qu = np.array(quaternion,dtype=float)
else np.array(quaternion,dtype=float)
if P > 0: qu[1:4] *= -1 # convert from P=1 to P=-1 if P > 0: qu[...,1:4] *= -1 # convert from P=1 to P=-1
if qu[0] < 0.0: if acceptHomomorph:
if acceptHomomorph: qu[qu[...,0] < 0.0] *= -1
qu *= -1. else:
else: if np.any(qu[...,0] < 0.0):
raise ValueError('Quaternion has negative first component: {}.'.format(qu[0])) raise ValueError('Quaternions need to have positive first(real) component.')
if not np.isclose(np.linalg.norm(qu), 1.0): if not np.all(np.isclose(np.linalg.norm(qu,axis=-1), 1.0)):
raise ValueError('Quaternion is not of unit length: {} {} {} {}.'.format(*qu)) raise ValueError('Quaternions need to have unit length.')
return Rotation(qu) return Rotation(qu)
@staticmethod @staticmethod
def fromEulers(eulers, def from_Eulers(eulers,
degrees = False): degrees = False):
eu = np.array(eulers,dtype=float)
eu = eulers if isinstance(eulers, np.ndarray) and eulers.dtype == np.dtype(float) \
else np.array(eulers,dtype=float)
eu = np.radians(eu) if degrees else eu eu = np.radians(eu) if degrees else eu
if np.any(eu < 0.0) or np.any(eu > 2.0*np.pi) or eu[1] > np.pi: if np.any(eu < 0.0) or np.any(eu > 2.0*np.pi) or np.any(eu[...,1] > np.pi): # ToDo: No separate check for PHI
raise ValueError('Euler angles outside of [0..2π],[0..π],[0..2π]: {} {} {}.'.format(*eu)) raise ValueError('Euler angles need to be in [0..2π],[0..π],[0..2π].')
return Rotation(Rotation.eu2qu(eu)) return Rotation(Rotation.eu2qu(eu))
@staticmethod @staticmethod
def fromAxisAngle(angleAxis, def from_axis_angle(axis_angle,
degrees = False, degrees = False,
normalise = False, normalise = False,
P = -1): P = -1):
ax = angleAxis if isinstance(angleAxis, np.ndarray) and angleAxis.dtype == np.dtype(float) \ ax = np.array(axis_angle,dtype=float)
else np.array(angleAxis,dtype=float)
if P > 0: ax[0:3] *= -1 # convert from P=1 to P=-1 if P > 0: ax[...,0:3] *= -1 # convert from P=1 to P=-1
if degrees: ax[ 3] = np.radians(ax[3]) if degrees: ax[..., 3] = np.radians(ax[...,3])
if normalise: ax[0:3] /= np.linalg.norm(ax[0:3]) if normalise: ax[...,0:3] /= np.linalg.norm(ax[...,0:3],axis=-1)
if ax[3] < 0.0 or ax[3] > np.pi: if np.any(ax[...,3] < 0.0) or np.any(ax[...,3] > np.pi):
raise ValueError('Axis angle rotation angle outside of [0..π]: {}.'.format(ax[3])) raise ValueError('Axis angle rotation angle outside of [0..π].')
if not np.isclose(np.linalg.norm(ax[0:3]), 1.0): if not np.all(np.isclose(np.linalg.norm(ax[...,0:3],axis=-1), 1.0)):
raise ValueError('Axis angle rotation axis is not of unit length: {} {} {}.'.format(*ax[0:3])) raise ValueError('Axis angle rotation axis is not of unit length.')
return Rotation(Rotation.ax2qu(ax)) return Rotation(Rotation.ax2qu(ax))
@staticmethod @staticmethod
def fromBasis(basis, def from_basis(basis,
orthonormal = True, orthonormal = True,
reciprocal = False, reciprocal = False):
):
om = np.array(basis,dtype=float)
om = basis if isinstance(basis, np.ndarray) else np.array(basis).reshape(3,3)
if reciprocal: if reciprocal:
om = np.linalg.inv(om.T/np.pi) # transform reciprocal basis set om = np.linalg.inv(mechanics.transpose(om)/np.pi) # transform reciprocal basis set
orthonormal = False # contains stretch orthonormal = False # contains stretch
if not orthonormal: if not orthonormal:
(U,S,Vh) = np.linalg.svd(om) # singular value decomposition (U,S,Vh) = np.linalg.svd(om) # singular value decomposition
om = np.dot(U,Vh) om = np.einsum('...ij,...jl->...il',U,Vh)
if not np.isclose(np.linalg.det(om),1.0): if not np.all(np.isclose(np.linalg.det(om),1.0)):
raise ValueError('matrix is not a proper rotation: {}.'.format(om)) raise ValueError('matrix is not a proper rotation: {}.'.format(om))
if not np.isclose(np.dot(om[0],om[1]), 0.0) \ if not np.all(np.isclose(np.einsum('...i,...i',om[...,0],om[...,1]), 0.0)) \
or not np.isclose(np.dot(om[1],om[2]), 0.0) \ or not np.all(np.isclose(np.einsum('...i,...i',om[...,1],om[...,2]), 0.0)) \
or not np.isclose(np.dot(om[2],om[0]), 0.0): or not np.all(np.isclose(np.einsum('...i,...i',om[...,2],om[...,0]), 0.0)):
raise ValueError('matrix is not orthogonal: {}.'.format(om)) raise ValueError('matrix is not orthogonal.')
return Rotation(Rotation.om2qu(om)) return Rotation(Rotation.om2qu(om))
@staticmethod @staticmethod
def fromMatrix(om, def from_matrix(om):
):
return Rotation.fromBasis(om) return Rotation.from_basis(om)
@staticmethod @staticmethod
def fromRodrigues(rodrigues, def from_Rodrigues(rodrigues,
normalise = False, normalise = False,
P = -1): P = -1):
ro = rodrigues if isinstance(rodrigues, np.ndarray) and rodrigues.dtype == np.dtype(float) \ ro = np.array(rodrigues,dtype=float)
else np.array(rodrigues,dtype=float)
if P > 0: ro[0:3] *= -1 # convert from P=1 to P=-1 if P > 0: ro[...,0:3] *= -1 # convert from P=1 to P=-1
if normalise: ro[0:3] /= np.linalg.norm(ro[0:3]) if normalise: ro[...,0:3] /= np.linalg.norm(ro[...,0:3],axis=-1)
if not np.isclose(np.linalg.norm(ro[0:3]), 1.0): if np.any(ro[...,3] < 0.0):
raise ValueError('Rodrigues rotation axis is not of unit length: {} {} {}.'.format(*ro[0:3])) raise ValueError('Rodrigues rotation angle not positive.')
if ro[3] < 0.0: if not np.all(np.isclose(np.linalg.norm(ro[...,0:3],axis=-1), 1.0)):
raise ValueError('Rodrigues rotation angle not positive: {}.'.format(ro[3])) raise ValueError('Rodrigues rotation axis is not of unit length.')
return Rotation(Rotation.ro2qu(ro)) return Rotation(Rotation.ro2qu(ro))
@staticmethod @staticmethod
def fromHomochoric(homochoric, def from_homochoric(homochoric,
P = -1): P = -1):
ho = np.array(homochoric,dtype=float)
ho = homochoric if isinstance(homochoric, np.ndarray) and homochoric.dtype == np.dtype(float) \
else np.array(homochoric,dtype=float)
if P > 0: ho *= -1 # convert from P=1 to P=-1 if P > 0: ho *= -1 # convert from P=1 to P=-1
if np.linalg.norm(ho) > (3.*np.pi/4.)**(1./3.)+1e-9: if np.any(np.linalg.norm(ho,axis=-1) > (3.*np.pi/4.)**(1./3.)+1e-9):
raise ValueError('Coordinate outside of the sphere: {} {} {}.'.format(ho)) raise ValueError('Coordinate outside of the sphere.')
return Rotation(Rotation.ho2qu(ho)) return Rotation(Rotation.ho2qu(ho))
@ -359,8 +361,7 @@ class Rotation:
def fromCubochoric(cubochoric, def fromCubochoric(cubochoric,
P = -1): P = -1):
cu = cubochoric if isinstance(cubochoric, np.ndarray) and cubochoric.dtype == np.dtype(float) \ cu = np.array(cubochoric,dtype=float)
else np.array(cubochoric,dtype=float)
if np.abs(np.max(cu))>np.pi**(2./3.) * 0.5+1e-9: if np.abs(np.max(cu))>np.pi**(2./3.) * 0.5+1e-9:
raise ValueError('Coordinate outside of the cube: {} {} {}.'.format(*cu)) raise ValueError('Coordinate outside of the cube: {} {} {}.'.format(*cu))
@ -403,17 +404,28 @@ class Rotation:
return Rotation.fromQuaternion(np.real(vec.T[eig.argmax()]),acceptHomomorph = True) return Rotation.fromQuaternion(np.real(vec.T[eig.argmax()]),acceptHomomorph = True)
@staticmethod @staticmethod
def fromRandom(): def from_random(shape=None):
r = np.random.random(3) if shape is None:
A = np.sqrt(r[2]) r = np.random.random(3)
B = np.sqrt(1.0-r[2]) else:
return Rotation(np.array([np.cos(2.0*np.pi*r[0])*A, r = np.random.random(tuple(shape)+(3,))
np.sin(2.0*np.pi*r[1])*B, A = np.sqrt(r[...,2])
np.cos(2.0*np.pi*r[1])*B, B = np.sqrt(1.0-r[...,2])
np.sin(2.0*np.pi*r[0])*A])).standardize() return Rotation(np.block([np.cos(2.0*np.pi*r[...,0])*A,
np.sin(2.0*np.pi*r[...,1])*B,
np.cos(2.0*np.pi*r[...,1])*B,
np.sin(2.0*np.pi*r[...,0])*A])).standardize()
# for compatibility (old names do not follow convention)
fromQuaternion = from_quaternion
fromEulers = from_Eulers
fromAxisAngle = from_axis_angle
fromBasis = from_basis
fromMatrix = from_matrix
fromRodrigues = from_Rodrigues
fromHomochoric = from_homochoric
fromRandom = from_random
#################################################################################################### ####################################################################################################
# Code below available according to the following conditions on https://github.com/MarDiehl/3Drotations # Code below available according to the following conditions on https://github.com/MarDiehl/3Drotations

12
python/damask/mechanics.py
View File

@ -135,16 +135,16 @@ def PK2(P,F):
Parameters Parameters
---------- ----------
P : numpy.ndarray of shape (:,3,3) or (3,3) P : numpy.ndarray of shape (...,3,3) or (3,3)
First Piola-Kirchhoff stress. First Piola-Kirchhoff stress.
F : numpy.ndarray of shape (:,3,3) or (3,3) F : numpy.ndarray of shape (...,3,3) or (3,3)
Deformation gradient. Deformation gradient.
""" """
if np.shape(F) == np.shape(P) == (3,3): if np.shape(F) == np.shape(P) == (3,3):
S = np.dot(np.linalg.inv(F),P) S = np.dot(np.linalg.inv(F),P)
else: else:
S = np.einsum('ijk,ikl->ijl',np.linalg.inv(F),P) S = np.einsum('...jk,...kl->...jl',np.linalg.inv(F),P)
return symmetric(S) return symmetric(S)
@ -241,7 +241,7 @@ def symmetric(T):
Parameters Parameters
---------- ----------
T : numpy.ndarray of shape (:,3,3) or (3,3) T : numpy.ndarray of shape (...,3,3) or (3,3)
Tensor of which the symmetrized values are computed. Tensor of which the symmetrized values are computed.
""" """
@ -254,12 +254,12 @@ def transpose(T):
Parameters Parameters
---------- ----------
T : numpy.ndarray of shape (:,3,3) or (3,3) T : numpy.ndarray of shape (...,3,3) or (3,3)
Tensor of which the transpose is computed. Tensor of which the transpose is computed.
""" """
return T.T if np.shape(T) == (3,3) else \ return T.T if np.shape(T) == (3,3) else \
np.transpose(T,(0,2,1)) np.swapaxes(T,axis2=-2,axis1=-1)
def _polar_decomposition(T,requested): def _polar_decomposition(T,requested):

13
python/tests/test_Rotation.py
View File

@ -157,6 +157,19 @@ class TestRotation:
print(m,o,rot.asQuaternion()) print(m,o,rot.asQuaternion())
assert ok and o.max() < np.pi**(2./3.)*0.5+1.e-9 assert ok and o.max() < np.pi**(2./3.)*0.5+1.e-9
@pytest.mark.parametrize('function,invalid',[(Rotation.from_quaternion, np.array([-1,0,0,0])),
(Rotation.from_quaternion, np.array([1,1,1,0])),
(Rotation.from_Eulers, np.array([1,4,0])),
(Rotation.from_axis_angle, np.array([1,0,0,4])),
(Rotation.from_axis_angle, np.array([1,1,0,1])),
(Rotation.from_matrix, np.random.rand(3,3)),
(Rotation.from_Rodrigues, np.array([1,0,0,-1])),
(Rotation.from_Rodrigues, np.array([1,1,0,1])),
(Rotation.from_homochoric, np.array([2,2,2])) ])
def test_invalid(self,function,invalid):
with pytest.raises(ValueError):
function(invalid)
@pytest.mark.parametrize('conversion',[Rotation.qu2om, @pytest.mark.parametrize('conversion',[Rotation.qu2om,
Rotation.qu2eu, Rotation.qu2eu,
Rotation.qu2ax, Rotation.qu2ax,