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

bugfix and further test

This commit is contained in:
Martin Diehl 2020-05-18 18:41:48 +02:00
parent 3dc90ddb94
commit b200894a40
3 changed files with 49 additions and 12 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

41
python/damask/_rotation.py
View File

@ -534,13 +534,13 @@ class Rotation:
def qu2om(qu):
qq = qu[...,0:1]**2-(qu[...,1:2]**2 + qu[...,2:3]**2 + qu[...,3:4]**2)
om = np.block([qq + 2.0*qu[...,1:2]**2,
2.0*(qu[...,2:3]*qu[...,1:2]+_P*qu[...,0:1]*qu[...,3:4]),
2.0*(qu[...,3:4]*qu[...,1:2]-_P*qu[...,0:1]*qu[...,2:3]),
2.0*(qu[...,1:2]*qu[...,2:3]-_P*qu[...,0:1]*qu[...,3:4]),
2.0*(qu[...,2:3]*qu[...,1:2]-_P*qu[...,0:1]*qu[...,3:4]),
2.0*(qu[...,3:4]*qu[...,1:2]+_P*qu[...,0:1]*qu[...,2:3]),
2.0*(qu[...,1:2]*qu[...,2:3]+_P*qu[...,0:1]*qu[...,3:4]),
qq + 2.0*qu[...,2:3]**2,
2.0*(qu[...,3:4]*qu[...,2:3]+_P*qu[...,0:1]*qu[...,1:2]),
2.0*(qu[...,1:2]*qu[...,3:4]+_P*qu[...,0:1]*qu[...,2:3]),
2.0*(qu[...,2:3]*qu[...,3:4]-_P*qu[...,0:1]*qu[...,1:2]),
2.0*(qu[...,3:4]*qu[...,2:3]-_P*qu[...,0:1]*qu[...,1:2]),
2.0*(qu[...,1:2]*qu[...,3:4]-_P*qu[...,0:1]*qu[...,2:3]),
2.0*(qu[...,2:3]*qu[...,3:4]+_P*qu[...,0:1]*qu[...,1:2]),
qq + 2.0*qu[...,3:4]**2,
]).reshape(qu.shape[:-1]+(3,3))
return om
@ -626,9 +626,31 @@ class Rotation:
"""
Rotation matrix to quaternion.
The original formulation (direct conversion) had (numerical?) issues
This formulation is from www.euclideanspace.com/maths/geometry/rotations/conversions/matrixToQuaternion.
The original formulation had issues.
"""
return Rotation.eu2qu(Rotation.om2eu(om))
def x(a):
trace = a[0,0] + a[1,1] + a[2,2]
if trace > 0:
s = 0.5 / np.sqrt(trace+ 1.0)
qu = np.array([0.25 / s,( a[2,1] - a[1,2] ) * s,( a[0,2] - a[2,0] ) * s,( a[1,0] - a[0,1] ) * s])
else:
if ( a[0,0] > a[1,1] and a[0,0] > a[2,2] ):
s = 2.0 * np.sqrt( 1.0 + a[0,0] - a[1,1] - a[2,2])
qu = np.array([ (a[2,1] - a[1,2]) / s,0.25 * s,(a[0,1] + a[1,0]) / s,(a[0,2] + a[2,0]) / s])
elif (a[1,1] > a[2,2]):
s = 2.0 * np.sqrt( 1.0 + a[1,1] - a[0,0] - a[2,2])
qu = np.array([ (a[0,2] - a[2,0]) / s,(a[0,1] + a[1,0]) / s,0.25 * s,(a[1,2] + a[2,1]) / s])
else:
s = 2.0 * np.sqrt( 1.0 + a[2,2] - a[0,0] - a[1,1] )
qu = np.array([ (a[1,0] - a[0,1]) / s,(a[0,2] + a[2,0]) / s,(a[1,2] + a[2,1]) / s,0.25 * s])
return qu
if len(om.shape) > 2:
om_ = om.reshape(-1,3,3)
qu = np.array([x(o) for o in om_]).reshape(om.shape[:-2]+(4,))
else:
qu = x(om)
return qu*np.array([1,_P,_P,_P])
@staticmethod
def om2eu(om):
@ -649,7 +671,6 @@ class Rotation:
eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu)
return eu
@staticmethod
def om2ax(om):
"""Rotation matrix to axis angle pair."""
@ -672,7 +693,6 @@ class Rotation:
ax[np.abs(ax[...,3])<1.e-6] = [ 0.0, 0.0, 1.0, 0.0]
return ax
@staticmethod
def om2ro(om):
"""Rotation matrix to Rodrigues-Frank vector."""
@ -703,7 +723,6 @@ class Rotation:
qu[qu[...,0]<0.0]*=-1
return qu
@staticmethod
def eu2om(eu):
"""Bunge-Euler angles to rotation matrix."""

17
python/tests/rotation_conversion.py
View File

@ -112,6 +112,23 @@ def qu2ho(qu):
#---------- Rotation matrix ----------
def om2qu(a):
trace = a[0,0] + a[1,1] + a[2,2]
if trace > 0:
s = 0.5 / np.sqrt(trace+ 1.0)
qu = np.array([0.25 / s,( a[2,1] - a[1,2] ) * s,( a[0,2] - a[2,0] ) * s,( a[1,0] - a[0,1] ) * s])
else:
if ( a[0,0] > a[1,1] and a[0,0] > a[2,2] ):
s = 2.0 * np.sqrt( 1.0 + a[0,0] - a[1,1] - a[2,2])
qu = np.array([ (a[2,1] - a[1,2]) / s,0.25 * s,(a[0,1] + a[1,0]) / s,(a[0,2] + a[2,0]) / s])
elif (a[1,1] > a[2,2]):
s = 2.0 * np.sqrt( 1.0 + a[1,1] - a[0,0] - a[2,2])
qu = np.array([ (a[0,2] - a[2,0]) / s,(a[0,1] + a[1,0]) / s,0.25 * s,(a[1,2] + a[2,1]) / s])
else:
s = 2.0 * np.sqrt( 1.0 + a[2,2] - a[0,0] - a[1,1] )
qu = np.array([ (a[1,0] - a[0,1]) / s,(a[0,2] + a[2,0]) / s,(a[1,2] + a[2,1]) / s,0.25 * s])
return qu
def om2eu(om):
"""Rotation matrix to Bunge-Euler angles."""
if not np.isclose(np.abs(om[2,2]),1.0,1.e-4):

3
python/tests/test_Rotation.py
View File

@ -196,7 +196,8 @@ class TestRotation:
assert np.allclose(single(q),c) and np.allclose(single(q),vectorized(q))
@pytest.mark.parametrize('vectorized, single',[(Rotation.om2eu,rotation_conversion.om2eu),
@pytest.mark.parametrize('vectorized, single',[(Rotation.om2qu,rotation_conversion.om2qu),
(Rotation.om2eu,rotation_conversion.om2eu),
(Rotation.om2ax,rotation_conversion.om2ax)])
def test_matrix_vectorization(self,default,vectorized,single):
om = np.array([rot.as_matrix() for rot in default])