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

Merge branch 'polish-rotation-error-determination' into 'development' 

cleaner and more efficient determination of ValueErrors

See merge request damask/DAMASK!656
This commit is contained in:
Martin Diehl 2022-11-17 17:52:14 +00:00
parent 70248c44cb d887f0a9d8
commit 05974fb300
2 changed files with 34 additions and 40 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

70
python/damask/_rotation.py
View File

@ -767,10 +767,8 @@ class Rotation:
""" """
qu = np.array(q,dtype=float) qu = np.array(q,dtype=float)
if qu.shape[:-2:-1] != (4,): if qu.shape[:-2:-1] != (4,): raise ValueError('invalid shape')
raise ValueError('invalid shape') if abs(P) != 1: raise ValueError('P ∉ {-1,1}')
if abs(P) != 1:
raise ValueError('P ∉ {-1,1}')
qu[...,1:4] *= -P qu[...,1:4] *= -P
@ -778,7 +776,8 @@ class Rotation:
qu[qu[...,0] < 0.0] *= -1 qu[qu[...,0] < 0.0] *= -1
elif np.any(qu[...,0] < 0.0): elif np.any(qu[...,0] < 0.0):
raise ValueError('quaternion with negative 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,rtol=1e-8)):
if not np.allclose(np.linalg.norm(qu,axis=-1), 1.0,rtol=1e-8):
raise ValueError('quaternion is not of unit length') raise ValueError('quaternion is not of unit length')
return Rotation(qu) return Rotation(qu)
@ -803,11 +802,10 @@ class Rotation:
""" """
eu = np.array(phi,dtype=float) eu = np.array(phi,dtype=float)
if eu.shape[:-2:-1] != (3,): if eu.shape[:-2:-1] != (3,): raise ValueError('invalid shape')
raise ValueError('invalid shape')
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 np.any(eu[...,1] > np.pi): # ToDo: No separate check for PHI if np.any(eu < 0.0) or np.any(eu > np.pi*np.array([2.0,1.0,2.0])):
raise ValueError('Euler angles outside of [0..2π],[0..π],[0..2π]') raise ValueError('Euler angles outside of [0..2π],[0..π],[0..2π]')
return Rotation(Rotation._eu2qu(eu)) return Rotation(Rotation._eu2qu(eu))
@ -834,18 +832,17 @@ class Rotation:
""" """
ax = np.array(axis_angle,dtype=float) ax = np.array(axis_angle,dtype=float)
if ax.shape[:-2:-1] != (4,): if ax.shape[:-2:-1] != (4,): raise ValueError('invalid shape')
raise ValueError('invalid shape') if abs(P) != 1: raise ValueError('P ∉ {-1,1}')
if abs(P) != 1:
raise ValueError('P ∉ {-1,1}')
ax[...,0:3] *= -P ax[...,0:3] *= -P
if degrees: ax[..., 3] = np.radians(ax[...,3]) if degrees: ax[..., 3] = np.radians(ax[...,3])
if normalize: ax[...,0:3] /= np.linalg.norm(ax[...,0:3],axis=-1,keepdims=True)
if np.any(ax[...,3] < 0.0) or np.any(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..π]') 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)):
print(np.linalg.norm(ax[...,0:3],axis=-1)) if normalize:
ax[...,0:3] /= np.linalg.norm(ax[...,0:3],axis=-1,keepdims=True)
elif not np.allclose(np.linalg.norm(ax[...,0:3],axis=-1), 1.0):
raise ValueError('axis–angle rotation axis is not of unit length') raise ValueError('axis–angle rotation axis is not of unit length')
return Rotation(Rotation._ax2qu(ax)) return Rotation(Rotation._ax2qu(ax))
@ -868,22 +865,23 @@ class Rotation:
""" """
om = np.array(basis,dtype=float) om = np.array(basis,dtype=float)
if om.shape[-2:] != (3,3): if om.shape[-2:] != (3,3): raise ValueError('invalid shape')
raise ValueError('invalid shape')
if reciprocal: if reciprocal:
om = np.linalg.inv(tensor.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 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.einsum('...ij,...jl',U,Vh) om = np.einsum('...ij,...jl',U,Vh)
if not np.all(np.isclose(np.linalg.det(om),1.0)): elif not np.allclose(np.einsum('...i,...i',om[...,0],om[...,1]), 0.0) \
raise ValueError('orientation matrix has determinant ≠ 1') or not np.allclose(np.einsum('...i,...i',om[...,1],om[...,2]), 0.0) \
if not np.all(np.isclose(np.einsum('...i,...i',om[...,0],om[...,1]), 0.0)) \ or not np.allclose(np.einsum('...i,...i',om[...,2],om[...,0]), 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('orientation matrix is not orthogonal') raise ValueError('orientation matrix is not orthogonal')
if not np.allclose(np.linalg.det(om),1.0):
raise ValueError('orientation matrix has determinant ≠ 1')
return Rotation(Rotation._om2qu(om)) return Rotation(Rotation._om2qu(om))
@staticmethod @staticmethod
@ -946,16 +944,15 @@ class Rotation:
""" """
ro = np.array(rho,dtype=float) ro = np.array(rho,dtype=float)
if ro.shape[:-2:-1] != (4,): if ro.shape[:-2:-1] != (4,): raise ValueError('invalid shape')
raise ValueError('invalid shape') if abs(P) != 1: raise ValueError('P ∉ {-1,1}')
if abs(P) != 1:
raise ValueError('P ∉ {-1,1}')
ro[...,0:3] *= -P ro[...,0:3] *= -P
if normalize: ro[...,0:3] /= np.linalg.norm(ro[...,0:3],axis=-1,keepdims=True) if np.any(ro[...,3] < 0.0): raise ValueError('Rodrigues vector rotation angle is negative')
if np.any(ro[...,3] < 0.0):
raise ValueError('Rodrigues vector rotation angle is negative') if normalize:
if not np.all(np.isclose(np.linalg.norm(ro[...,0:3],axis=-1), 1.0)): ro[...,0:3] /= np.linalg.norm(ro[...,0:3],axis=-1,keepdims=True)
elif not np.allclose(np.linalg.norm(ro[...,0:3],axis=-1), 1.0):
raise ValueError('Rodrigues vector rotation axis is not of unit length') raise ValueError('Rodrigues vector rotation axis is not of unit length')
return Rotation(Rotation._ro2qu(ro)) return Rotation(Rotation._ro2qu(ro))
@ -975,10 +972,8 @@ class Rotation:
""" """
ho = np.array(h,dtype=float) ho = np.array(h,dtype=float)
if ho.shape[:-2:-1] != (3,): if ho.shape[:-2:-1] != (3,): raise ValueError('invalid shape')
raise ValueError('invalid shape') if abs(P) != 1: raise ValueError('P ∉ {-1,1}')
if abs(P) != 1:
raise ValueError('P ∉ {-1,1}')
ho *= -P ho *= -P
@ -1002,11 +997,8 @@ class Rotation:
""" """
cu = np.array(x,dtype=float) cu = np.array(x,dtype=float)
if cu.shape[:-2:-1] != (3,): if cu.shape[:-2:-1] != (3,): raise ValueError('invalid shape')
raise ValueError('invalid shape') if abs(P) != 1: raise ValueError('P ∉ {-1,1}')
if abs(P) != 1:
raise ValueError('P ∉ {-1,1}')
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('cubochoric coordinate outside of the cube') raise ValueError('cubochoric coordinate outside of the cube')

4
python/tests/test_Rotation.py
View File

@ -927,6 +927,7 @@ class TestRotation:
@pytest.mark.parametrize('function,invalid',[(Rotation.from_quaternion, np.array([-1,0,0,0])), @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_quaternion, np.array([1,1,1,0])),
(Rotation.from_Euler_angles, np.array([1,4,0])), (Rotation.from_Euler_angles, np.array([1,4,0])),
(Rotation.from_Euler_angles, np.array([-1,0,0])),
(Rotation.from_axis_angle, np.array([1,0,0,4])), (Rotation.from_axis_angle, np.array([1,0,0,4])),
(Rotation.from_axis_angle, np.array([1,1,0,1])), (Rotation.from_axis_angle, np.array([1,1,0,1])),
(Rotation.from_matrix, np.random.rand(3,3)), (Rotation.from_matrix, np.random.rand(3,3)),
@ -934,7 +935,8 @@ class TestRotation:
(Rotation.from_Rodrigues_vector, np.array([1,0,0,-1])), (Rotation.from_Rodrigues_vector, np.array([1,0,0,-1])),
(Rotation.from_Rodrigues_vector, np.array([1,1,0,1])), (Rotation.from_Rodrigues_vector, np.array([1,1,0,1])),
(Rotation.from_homochoric, np.array([2,2,2])), (Rotation.from_homochoric, np.array([2,2,2])),
(Rotation.from_cubochoric, np.array([1.1,0,0])) ]) (Rotation.from_cubochoric, np.array([1.1,0,0])),
])
def test_invalid_value(self,function,invalid): def test_invalid_value(self,function,invalid):
with pytest.raises(ValueError): with pytest.raises(ValueError):
function(invalid) function(invalid)