From 23fc58699f28468baacb057e6f1a393cd88aea7b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 15 Apr 2020 23:00:00 +0200 Subject: [PATCH 1/7] vectorize Rotation.fromXXX functions --- python/.coveragerc | 3 + python/damask/_rotation.py | 172 ++++++++++++++++++---------------- python/damask/mechanics.py | 12 +-- python/tests/test_Rotation.py | 13 +++ 4 files changed, 114 insertions(+), 86 deletions(-) diff --git a/python/.coveragerc b/python/.coveragerc index c712d2595..5daa25bb2 100644 --- a/python/.coveragerc +++ b/python/.coveragerc @@ -1,2 +1,5 @@ [run] omit = tests/* + damask/_asciitable.py + damask/_test.py + damask/config/* diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index dcd61822f..8a2441012 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -1,6 +1,7 @@ import numpy as np from ._Lambert import ball_to_cube, cube_to_ball +from . import mechanics _P = -1 @@ -61,6 +62,8 @@ class Rotation: def __repr__(self): """Orientation displayed as unit quaternion, rotation matrix, and Bunge-Euler angles.""" + if self.quaternion.shape != (4,): + raise NotImplementedError return '\n'.join([ 'Quaternion: (real={:.3f}, imag=<{:+.3f}, {:+.3f}, {:+.3f}>)'.format(*(self.quaternion)), 'Matrix:\n{}'.format(self.asMatrix()), @@ -116,7 +119,7 @@ class Rotation: def inverse(self): """In-place inverse rotation/backward rotation.""" - self.quaternion[1:] *= -1 + self.quaternion[...,1:] *= -1 return self def inversed(self): @@ -125,12 +128,12 @@ class Rotation: def standardize(self): - """In-place quaternion representation with positive q.""" - if self.quaternion[0] < 0.0: self.quaternion*=-1 + """In-place quaternion representation with positive real part.""" + self.quaternion[self.quaternion[...,0] < 0.0] *= -1 return self def standardized(self): - """Quaternion representation with positive q.""" + """Quaternion representation with positive real.""" return self.copy().standardize() @@ -165,7 +168,7 @@ class Rotation: 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 ---------- @@ -251,107 +254,106 @@ class Rotation: # static constructors. The input data needs to follow the convention, options allow to # relax these convections @staticmethod - def fromQuaternion(quaternion, - acceptHomomorph = False, - P = -1): + def from_quaternion(quaternion, + acceptHomomorph = False, + P = -1): - qu = quaternion if isinstance(quaternion,np.ndarray) and quaternion.dtype == np.dtype(float) \ - else np.array(quaternion,dtype=float) - if P > 0: qu[1:4] *= -1 # convert from P=1 to P=-1 - if qu[0] < 0.0: - if acceptHomomorph: - qu *= -1. - else: - raise ValueError('Quaternion has negative first component: {}.'.format(qu[0])) - if not np.isclose(np.linalg.norm(qu), 1.0): - raise ValueError('Quaternion is not of unit length: {} {} {} {}.'.format(*qu)) + qu = np.array(quaternion,dtype=float) + + if P > 0: qu[...,1:4] *= -1 # convert from P=1 to P=-1 + if acceptHomomorph: + qu[qu[...,0] < 0.0] *= -1 + else: + if np.any(qu[...,0] < 0.0): + raise ValueError('Quaternions need to have positive first(real) component.') + if not np.all(np.isclose(np.linalg.norm(qu,axis=-1), 1.0)): + raise ValueError('Quaternions need to have unit length.') return Rotation(qu) @staticmethod - def fromEulers(eulers, - degrees = False): + def from_Eulers(eulers, + 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 - if np.any(eu < 0.0) or np.any(eu > 2.0*np.pi) or eu[1] > np.pi: - raise ValueError('Euler angles outside of [0..2π],[0..π],[0..2π]: {} {} {}.'.format(*eu)) + 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 need to be in [0..2π],[0..π],[0..2π].') return Rotation(Rotation.eu2qu(eu)) @staticmethod - def fromAxisAngle(angleAxis, - degrees = False, - normalise = False, - P = -1): + def from_axis_angle(axis_angle, + degrees = False, + normalise = False, + P = -1): - ax = angleAxis if isinstance(angleAxis, np.ndarray) and angleAxis.dtype == np.dtype(float) \ - else np.array(angleAxis,dtype=float) - if P > 0: ax[0:3] *= -1 # convert from P=1 to P=-1 - if degrees: ax[ 3] = np.radians(ax[3]) - if normalise: ax[0:3] /= np.linalg.norm(ax[0:3]) - if ax[3] < 0.0 or ax[3] > np.pi: - raise ValueError('Axis angle rotation angle outside of [0..π]: {}.'.format(ax[3])) - if not np.isclose(np.linalg.norm(ax[0:3]), 1.0): - raise ValueError('Axis angle rotation axis is not of unit length: {} {} {}.'.format(*ax[0:3])) + ax = np.array(axis_angle,dtype=float) + + if P > 0: ax[...,0:3] *= -1 # convert from P=1 to P=-1 + if degrees: ax[..., 3] = np.radians(ax[...,3]) + if normalise: ax[...,0:3] /= np.linalg.norm(ax[...,0:3],axis=-1) + if np.any(ax[...,3] < 0.0) or np.any(ax[...,3] > np.pi): + raise ValueError('Axis angle rotation angle outside of [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.') return Rotation(Rotation.ax2qu(ax)) @staticmethod - def fromBasis(basis, - orthonormal = True, - reciprocal = False, - ): + def from_basis(basis, + orthonormal = True, + reciprocal = False): + + om = np.array(basis,dtype=float) - om = basis if isinstance(basis, np.ndarray) else np.array(basis).reshape(3,3) 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 if not orthonormal: (U,S,Vh) = np.linalg.svd(om) # singular value decomposition - om = np.dot(U,Vh) - if not np.isclose(np.linalg.det(om),1.0): + om = np.einsum('...ij,...jl->...il',U,Vh) + if not np.all(np.isclose(np.linalg.det(om),1.0)): raise ValueError('matrix is not a proper rotation: {}.'.format(om)) - if not np.isclose(np.dot(om[0],om[1]), 0.0) \ - or not np.isclose(np.dot(om[1],om[2]), 0.0) \ - or not np.isclose(np.dot(om[2],om[0]), 0.0): - raise ValueError('matrix is not orthogonal: {}.'.format(om)) + if not np.all(np.isclose(np.einsum('...i,...i',om[...,0],om[...,1]), 0.0)) \ + or not np.all(np.isclose(np.einsum('...i,...i',om[...,1],om[...,2]), 0.0)) \ + or not np.all(np.isclose(np.einsum('...i,...i',om[...,2],om[...,0]), 0.0)): + raise ValueError('matrix is not orthogonal.') return Rotation(Rotation.om2qu(om)) @staticmethod - def fromMatrix(om, - ): + def from_matrix(om): - return Rotation.fromBasis(om) + return Rotation.from_basis(om) @staticmethod - def fromRodrigues(rodrigues, - normalise = False, - P = -1): + def from_Rodrigues(rodrigues, + normalise = False, + P = -1): - ro = rodrigues if isinstance(rodrigues, np.ndarray) and rodrigues.dtype == np.dtype(float) \ - else np.array(rodrigues,dtype=float) - 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 not np.isclose(np.linalg.norm(ro[0:3]), 1.0): - raise ValueError('Rodrigues rotation axis is not of unit length: {} {} {}.'.format(*ro[0:3])) - if ro[3] < 0.0: - raise ValueError('Rodrigues rotation angle not positive: {}.'.format(ro[3])) + ro = np.array(rodrigues,dtype=float) + + 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],axis=-1) + if np.any(ro[...,3] < 0.0): + raise ValueError('Rodrigues rotation angle not positive.') + if not np.all(np.isclose(np.linalg.norm(ro[...,0:3],axis=-1), 1.0)): + raise ValueError('Rodrigues rotation axis is not of unit length.') return Rotation(Rotation.ro2qu(ro)) @staticmethod - def fromHomochoric(homochoric, - P = -1): + def from_homochoric(homochoric, + 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 np.linalg.norm(ho) > (3.*np.pi/4.)**(1./3.)+1e-9: - raise ValueError('Coordinate outside of the sphere: {} {} {}.'.format(ho)) + if np.any(np.linalg.norm(ho,axis=-1) > (3.*np.pi/4.)**(1./3.)+1e-9): + raise ValueError('Coordinate outside of the sphere.') return Rotation(Rotation.ho2qu(ho)) @@ -359,8 +361,7 @@ class Rotation: def fromCubochoric(cubochoric, P = -1): - cu = cubochoric if isinstance(cubochoric, np.ndarray) and cubochoric.dtype == np.dtype(float) \ - else np.array(cubochoric,dtype=float) + cu = np.array(cubochoric,dtype=float) if np.abs(np.max(cu))>np.pi**(2./3.) * 0.5+1e-9: 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) - @staticmethod - def fromRandom(): - r = np.random.random(3) - A = np.sqrt(r[2]) - B = np.sqrt(1.0-r[2]) - return Rotation(np.array([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() + def from_random(shape=None): + if shape is None: + r = np.random.random(3) + else: + r = np.random.random(tuple(shape)+(3,)) + A = np.sqrt(r[...,2]) + B = np.sqrt(1.0-r[...,2]) + 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 diff --git a/python/damask/mechanics.py b/python/damask/mechanics.py index 674ff9c5a..8d19e9b85 100644 --- a/python/damask/mechanics.py +++ b/python/damask/mechanics.py @@ -135,16 +135,16 @@ def PK2(P,F): 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. - F : numpy.ndarray of shape (:,3,3) or (3,3) + F : numpy.ndarray of shape (...,3,3) or (3,3) Deformation gradient. """ if np.shape(F) == np.shape(P) == (3,3): S = np.dot(np.linalg.inv(F),P) 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) @@ -241,7 +241,7 @@ def symmetric(T): 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. """ @@ -254,12 +254,12 @@ def transpose(T): 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. """ 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): diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index 2ac819f4c..8150e35f8 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -157,6 +157,19 @@ class TestRotation: print(m,o,rot.asQuaternion()) 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, Rotation.qu2eu, Rotation.qu2ax, From ae3eca5f98a21c745c3dc3ba780517a5848c8bdc Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 21 Apr 2020 11:59:42 +0200 Subject: [PATCH 2/7] fix for vectorized from_random --- python/damask/_rotation.py | 45 ++++++++++++++++++++------------------ 1 file changed, 24 insertions(+), 21 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 8a2441012..69a725e99 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -265,9 +265,9 @@ class Rotation: qu[qu[...,0] < 0.0] *= -1 else: if np.any(qu[...,0] < 0.0): - raise ValueError('Quaternions need to have positive first(real) component.') + raise ValueError('Quaternion with negative first (real) component.') if not np.all(np.isclose(np.linalg.norm(qu,axis=-1), 1.0)): - raise ValueError('Quaternions need to have unit length.') + raise ValueError('Quaternion is not of unit length.') return Rotation(qu) @@ -279,7 +279,7 @@ class Rotation: eu = np.radians(eu) if degrees else eu 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 need to be in [0..2π],[0..π],[0..2π].') + raise ValueError('Euler angles outside of [0..2π],[0..π],[0..2π].') return Rotation(Rotation.eu2qu(eu)) @@ -315,11 +315,11 @@ class Rotation: (U,S,Vh) = np.linalg.svd(om) # singular value decomposition om = np.einsum('...ij,...jl->...il',U,Vh) if not np.all(np.isclose(np.linalg.det(om),1.0)): - raise ValueError('matrix is not a proper rotation: {}.'.format(om)) + raise ValueError('Orientation matrix has determinant ≠ 1.') if not np.all(np.isclose(np.einsum('...i,...i',om[...,0],om[...,1]), 0.0)) \ or not np.all(np.isclose(np.einsum('...i,...i',om[...,1],om[...,2]), 0.0)) \ or not np.all(np.isclose(np.einsum('...i,...i',om[...,2],om[...,0]), 0.0)): - raise ValueError('matrix is not orthogonal.') + raise ValueError('Orientation matrix is not orthogonal.') return Rotation(Rotation.om2qu(om)) @@ -338,9 +338,9 @@ class Rotation: 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],axis=-1) if np.any(ro[...,3] < 0.0): - raise ValueError('Rodrigues rotation angle not positive.') + raise ValueError('Rodrigues vector rotation angle not positive.') if not np.all(np.isclose(np.linalg.norm(ro[...,0:3],axis=-1), 1.0)): - raise ValueError('Rodrigues rotation axis is not of unit length.') + raise ValueError('Rodrigues vector rotation axis is not of unit length.') return Rotation(Rotation.ro2qu(ro)) @@ -353,7 +353,7 @@ class Rotation: if P > 0: ho *= -1 # convert from P=1 to P=-1 if np.any(np.linalg.norm(ho,axis=-1) > (3.*np.pi/4.)**(1./3.)+1e-9): - raise ValueError('Coordinate outside of the sphere.') + raise ValueError('Homochoric coordinate outside of the sphere.') return Rotation(Rotation.ho2qu(ho)) @@ -364,7 +364,7 @@ class Rotation: cu = np.array(cubochoric,dtype=float) 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('Cubochoric coordinate outside of the cube: {} {} {}.'.format(*cu)) ho = Rotation.cu2ho(cu) if P > 0: ho *= -1 # convert from P=1 to P=-1 @@ -408,14 +408,20 @@ class Rotation: def from_random(shape=None): if shape is None: r = np.random.random(3) - else: + elif hasattr(shape, '__iter__'): r = np.random.random(tuple(shape)+(3,)) + else: + r = np.random.random((shape,3)) + A = np.sqrt(r[...,2]) B = np.sqrt(1.0-r[...,2]) - 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() + q = np.stack([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],axis=-1) + + return Rotation(q.reshape(r.shape[:-1]+(4,)) if shape is not None else q).standardize() + # for compatibility (old names do not follow convention) fromQuaternion = from_quaternion @@ -820,12 +826,11 @@ class Rotation: c = np.cos(ax[3]*0.5) s = np.sin(ax[3]*0.5) qu = np.array([ c, ax[0]*s, ax[1]*s, ax[2]*s ]) - return qu else: c = np.cos(ax[...,3:4]*.5) s = np.sin(ax[...,3:4]*.5) qu = np.where(np.abs(ax[...,3:4])<1.e-6,[1.0, 0.0, 0.0, 0.0],np.block([c, ax[...,:3]*s])) - return qu + return qu @staticmethod def ax2om(ax): @@ -871,7 +876,7 @@ class Rotation: # 180 degree case ro += [np.inf] if np.isclose(ax[3],np.pi,atol=1.0e-15,rtol=0.0) else \ [np.tan(ax[3]*0.5)] - return np.array(ro) + ro = np.array(ro) else: ro = np.block([ax[...,:3], np.where(np.isclose(ax[...,3:4],np.pi,atol=1.e-15,rtol=.0), @@ -879,7 +884,7 @@ class Rotation: np.tan(ax[...,3:4]*0.5)) ]) ro[np.abs(ax[...,3])<1.e-6] = [.0,.0,_P,.0] - return ro + return ro @staticmethod def ax2ho(ax): @@ -887,11 +892,10 @@ class Rotation: if len(ax.shape) == 1: f = (0.75 * ( ax[3] - np.sin(ax[3]) ))**(1.0/3.0) ho = ax[0:3] * f - return ho else: f = (0.75 * ( ax[...,3:4] - np.sin(ax[...,3:4]) ))**(1.0/3.0) ho = ax[...,:3] * f - return ho + return ho @staticmethod def ax2cu(ax): @@ -948,7 +952,6 @@ class Rotation: f = np.where(np.isfinite(ro[...,3:4]),2.0*np.arctan(ro[...,3:4]) -np.sin(2.0*np.arctan(ro[...,3:4])),np.pi) ho = np.where(np.broadcast_to(np.sum(ro[...,0:3]**2.0,axis=-1,keepdims=True) < 1.e-6,ro[...,0:3].shape), np.zeros(3), ro[...,0:3]* (0.75*f)**(1.0/3.0)) - return ho @staticmethod From 75d723837650bdb89b68b5eddfd1b41e583303e9 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 21 Apr 2020 12:27:50 +0200 Subject: [PATCH 3/7] vectorized as_XXX --- python/damask/_rotation.py | 42 +++++++++++++++++++++++--------------- 1 file changed, 25 insertions(+), 17 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 69a725e99..b9be2a606 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -166,7 +166,7 @@ class Rotation: ################################################################################################ # convert to different orientation representations (numpy arrays) - def asQuaternion(self): + def as_quaternion(self): """ Unit quaternion [q, p_1, p_2, p_3]. @@ -178,8 +178,8 @@ class Rotation: """ return self.quaternion - def asEulers(self, - degrees = False): + def as_Eulers(self, + degrees = False): """ Bunge-Euler angles: (φ_1, ϕ, φ_2). @@ -193,9 +193,9 @@ class Rotation: if degrees: eu = np.degrees(eu) return eu - def asAxisAngle(self, - degrees = False, - pair = False): + def as_axis_angle(self, + degrees = False, + pair = False): """ Axis angle representation [n_1, n_2, n_3, ω] unless pair == True: ([n_1, n_2, n_3], ω). @@ -208,15 +208,15 @@ class Rotation: """ ax = Rotation.qu2ax(self.quaternion) - if degrees: ax[3] = np.degrees(ax[3]) - return (ax[:3],np.degrees(ax[3])) if pair else ax + if degrees: ax[...,3] = np.degrees(ax[...,3]) + return (ax[...,:3],ax[...,3]) if pair else ax - def asMatrix(self): + def as_matrix(self): """Rotation matrix.""" return Rotation.qu2om(self.quaternion) - def asRodrigues(self, - vector = False): + def as_Rodrigues(self, + vector = False): """ Rodrigues-Frank vector representation [n_1, n_2, n_3, tan(ω/2)] unless vector == True: [n_1, n_2, n_3] * tan(ω/2). @@ -227,9 +227,9 @@ class Rotation: """ ro = Rotation.qu2ro(self.quaternion) - return ro[:3]*ro[3] if vector else ro + return ro[...,:3]*ro[...,3] if vector else ro - def asHomochoric(self): + def as_homochoric(self): """Homochoric vector: (h_1, h_2, h_3).""" return Rotation.qu2ho(self.quaternion) @@ -237,7 +237,7 @@ class Rotation: """Cubochoric vector: (c_1, c_2, c_3).""" return Rotation.qu2cu(self.quaternion) - def asM(self): + def M(self): # ToDo not sure about the name: as_M or M? we do not have a from_M """ Intermediate representation supporting quaternion averaging. @@ -247,12 +247,20 @@ class Rotation: https://doi.org/10.2514/1.28949 """ - return np.outer(self.quaternion,self.quaternion) + return np.einsum('...i,...j',self.quaternion,self.quaternion) + # for compatibility (old names do not follow convention) + asM = M + asQuaternion = as_quaternion + asEulers = as_Eulers + asAxisAngle = as_axis_angle + asMatrix = as_matrix + asRodrigues = as_Rodrigues + asHomochoric = as_homochoric ################################################################################################ - # static constructors. The input data needs to follow the convention, options allow to - # relax these convections + # Static constructors. The input data needs to follow the conventions, options allow to + # relax the conventions. @staticmethod def from_quaternion(quaternion, acceptHomomorph = False, From b1be4e7ac80a122cde6df4f69970374f180679aa Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 21 Apr 2020 13:22:24 +0200 Subject: [PATCH 4/7] rotation class does not take care of correct shape anymore --- processing/post/addOrientations.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/processing/post/addOrientations.py b/processing/post/addOrientations.py index 137b8121e..e1752b0a4 100755 --- a/processing/post/addOrientations.py +++ b/processing/post/addOrientations.py @@ -172,7 +172,7 @@ for name in filenames: elif inputtype == 'matrix': d = representations['matrix'][1] - o = damask.Rotation.fromMatrix(list(map(float,table.data[column:column+d]))) + o = damask.Rotation.fromMatrix(np.array(list(map(float,table.data[column:column+d]))).reshape(3,3)) elif inputtype == 'frame': M = np.array(list(map(float,table.data[column[0]:column[0]+3] + \ From 97a5880d76311c83e5864f642f24e1d3e51d14cc Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 21 Apr 2020 13:22:55 +0200 Subject: [PATCH 5/7] ensure correct shape --- python/damask/_rotation.py | 14 ++++++++++++++ python/tests/test_Rotation.py | 11 +++++++++++ 2 files changed, 25 insertions(+) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index b9be2a606..6076b4add 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -267,6 +267,8 @@ class Rotation: P = -1): qu = np.array(quaternion,dtype=float) + if qu.shape[:-2:-1] != (4,): + raise ValueError('Invalid shape.') if P > 0: qu[...,1:4] *= -1 # convert from P=1 to P=-1 if acceptHomomorph: @@ -284,6 +286,8 @@ class Rotation: degrees = False): eu = np.array(eulers,dtype=float) + if eu.shape[:-2:-1] != (3,): + raise ValueError('Invalid shape.') eu = np.radians(eu) if degrees else eu 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 @@ -298,6 +302,8 @@ class Rotation: P = -1): ax = np.array(axis_angle,dtype=float) + if ax.shape[:-2:-1] != (4,): + raise ValueError('Invalid shape.') if P > 0: ax[...,0:3] *= -1 # convert from P=1 to P=-1 if degrees: ax[..., 3] = np.radians(ax[...,3]) @@ -315,6 +321,8 @@ class Rotation: reciprocal = False): om = np.array(basis,dtype=float) + if om.shape[:-3:-1] != (3,3): + raise ValueError('Invalid shape.') if reciprocal: om = np.linalg.inv(mechanics.transpose(om)/np.pi) # transform reciprocal basis set @@ -342,6 +350,8 @@ class Rotation: P = -1): ro = np.array(rodrigues,dtype=float) + if ro.shape[:-2:-1] != (4,): + raise ValueError('Invalid shape.') 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],axis=-1) @@ -357,6 +367,8 @@ class Rotation: P = -1): ho = np.array(homochoric,dtype=float) + if ho.shape[:-2:-1] != (3,): + raise ValueError('Invalid shape.') if P > 0: ho *= -1 # convert from P=1 to P=-1 @@ -370,6 +382,8 @@ class Rotation: P = -1): cu = np.array(cubochoric,dtype=float) + if cu.shape[:-2:-1] != (3,): + raise ValueError('Invalid shape.') if np.abs(np.max(cu))>np.pi**(2./3.) * 0.5+1e-9: raise ValueError('Cubochoric coordinate outside of the cube: {} {} {}.'.format(*cu)) diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index 8150e35f8..302e895ee 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -157,6 +157,17 @@ class TestRotation: print(m,o,rot.asQuaternion()) assert ok and o.max() < np.pi**(2./3.)*0.5+1.e-9 + @pytest.mark.parametrize('function',[Rotation.from_quaternion, + Rotation.from_Eulers, + Rotation.from_axis_angle, + Rotation.from_matrix, + Rotation.from_Rodrigues, + Rotation.from_homochoric]) + def test_invalid_shape(self,function): + invalid_shape = np.random.random(np.random.randint(8,32,(3))) + with pytest.raises(ValueError): + function(invalid_shape) + @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])), From 707324887f562c4449ff2091082abdcb1246a71c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 24 Apr 2020 06:52:09 +0200 Subject: [PATCH 6/7] inform the user about planned functionality --- python/damask/_orientation.py | 3 +++ python/damask/_rotation.py | 14 +++++++++----- 2 files changed, 12 insertions(+), 5 deletions(-) diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index 76475e057..e6c732007 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -38,6 +38,9 @@ class Orientation: else: self.rotation = Rotation.fromQuaternion(rotation) # assume quaternion + if self.rotation.quaternion.shape != (4,): + raise NotImplementedError('Support for multiple rotations missing') + def disorientation(self, other, SST = True, diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 6076b4add..716977a79 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -63,7 +63,7 @@ class Rotation: def __repr__(self): """Orientation displayed as unit quaternion, rotation matrix, and Bunge-Euler angles.""" if self.quaternion.shape != (4,): - raise NotImplementedError + raise NotImplementedError('Support for multiple rotations missing') return '\n'.join([ 'Quaternion: (real={:.3f}, imag=<{:+.3f}, {:+.3f}, {:+.3f}>)'.format(*(self.quaternion)), 'Matrix:\n{}'.format(self.asMatrix()), @@ -86,6 +86,8 @@ class Rotation: considere rotation of (3,3,3,3)-matrix """ + if self.quaternion.shape != (4,): + raise NotImplementedError('Support for multiple rotations missing') if isinstance(other, Rotation): # rotate a rotation self_q = self.quaternion[0] self_p = self.quaternion[1:] @@ -110,7 +112,7 @@ class Rotation: elif other.shape == (3,3,): # rotate a single (3x3)-matrix return np.dot(self.asMatrix(),np.dot(other,self.asMatrix().T)) elif other.shape == (3,3,3,3,): - raise NotImplementedError + raise NotImplementedError('Support for rotation of 4th order tensors missing') else: return NotImplemented else: @@ -133,7 +135,7 @@ class Rotation: return self def standardized(self): - """Quaternion representation with positive real.""" + """Quaternion representation with positive real part.""" return self.copy().standardize() @@ -160,6 +162,8 @@ class Rotation: Rotation from which the average is rotated. """ + if self.quaternion.shape != (4,) or other.quaternion.shape != (4,): + raise NotImplementedError('Support for multiple rotations missing') return Rotation.fromAverage([self,other]) @@ -1047,7 +1051,7 @@ class Rotation: if len(ho.shape) == 1: return ball_to_cube(ho) else: - raise NotImplementedError + raise NotImplementedError('Support for multiple rotations missing') #---------- Cubochoric ---------- @@ -1082,4 +1086,4 @@ class Rotation: if len(cu.shape) == 1: return cube_to_ball(cu) else: - raise NotImplementedError + raise NotImplementedError('Support for multiple rotations missing') From 97438713ccbed90d2da6ceedfcb3a80f082b8333 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 24 Apr 2020 12:55:59 +0200 Subject: [PATCH 7/7] require numpy array --- processing/post/addSchmidfactors.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/processing/post/addSchmidfactors.py b/processing/post/addSchmidfactors.py index 94c062766..b3c4c7500 100755 --- a/processing/post/addSchmidfactors.py +++ b/processing/post/addSchmidfactors.py @@ -214,7 +214,7 @@ for name in filenames: outputAlive = True while outputAlive and table.data_read(): # read next data line of ASCII table - o = damask.Rotation(list(map(float,table.data[column:column+4]))) + o = damask.Rotation(np.array(list(map(float,table.data[column:column+4])))) table.data_append( np.abs( np.sum(slip_direction * (o * force) ,axis=1) \ * np.sum(slip_normal * (o * normal),axis=1)))