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

Merge branch '326-faster-euler-to-quaternion' into 'development' 

adopt faster quaternion-to-Euler algorithm

Closes #326

See merge request damask/DAMASK!812
This commit is contained in:
Martin Diehl 2023-09-01 14:04:54 +00:00
parent 8728491988 9b0a051737
commit be87d9eed8
1 changed files with 34 additions and 21 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

55
python/damask/_rotation.py
View File

@ -1324,28 +1324,41 @@ class Rotation:
@staticmethod @staticmethod
def _qu2eu(qu: np.ndarray) -> np.ndarray: def _qu2eu(qu: np.ndarray) -> np.ndarray:
"""Quaternion to Bunge Euler angles.""" """
q02 = qu[...,0:1]*qu[...,2:3] Quaternion to Bunge Euler angles.
q13 = qu[...,1:2]*qu[...,3:4]
q01 = qu[...,0:1]*qu[...,1:2]
q23 = qu[...,2:3]*qu[...,3:4]
q03_s = qu[...,0:1]**2+qu[...,3:4]**2
q12_s = qu[...,1:2]**2+qu[...,2:3]**2
chi = np.sqrt(q03_s*q12_s)
eu = np.where(np.abs(q12_s) < 1.e-8, References
np.block([np.arctan2(-_P*2.*qu[...,0:1]*qu[...,3:4],qu[...,0:1]**2-qu[...,3:4]**2), ----------
np.zeros(qu.shape[:-1]+(2,))]), E. Bernardes and S. Viollet, PLoS ONE 17(11):e0276302, 2022
np.where(np.abs(q03_s) < 1.e-8, https://doi.org/10.1371/journal.pone.0276302
np.block([np.arctan2( 2.*qu[...,1:2]*qu[...,2:3],qu[...,1:2]**2-qu[...,2:3]**2),
np.broadcast_to(np.pi,qu[...,0:1].shape), """
np.zeros(qu.shape[:-1]+(1,))]), a = qu[...,0:1]
np.block([np.arctan2((-_P*q02+q13)*chi, (-_P*q01-q23)*chi), b = -_P*qu[...,3:4]
np.arctan2( 2.*chi, q03_s-q12_s ), c = -_P*qu[...,1:2]
np.arctan2(( _P*q02+q13)*chi, (-_P*q01+q23)*chi)]) d = -_P*qu[...,2:3]
)
) eu = np.block([
eu[np.abs(eu) < 1.e-6] = 0. np.arctan2(b,a),
np.arccos(2*(a**2+b**2)/(a**2+b**2+c**2+d**2)-1),
np.arctan2(-d,c),
])
eu_sum = eu[...,0] + eu[...,2]
eu_diff = eu[...,0] - eu[...,2]
is_zero = np.isclose(eu[...,1],0.0)
is_pi = np.isclose(eu[...,1],np.pi)
is_ok = ~np.logical_or(is_zero,is_pi)
eu[...,0][is_zero] = 2*eu[...,0][is_zero]
eu[...,0][is_pi] = -2*eu[...,2][is_pi]
eu[...,2][~is_ok] = 0.0
eu[...,0][is_ok] = eu_diff[is_ok]
eu[...,2][is_ok] = eu_sum [is_ok]
eu[np.logical_or(np.abs(eu) < 1.e-6,
np.abs(eu-2*np.pi) < 1.e-6)] = 0.
return np.where(eu < 0., eu%(np.pi*np.array([2.,1.,2.])),eu) return np.where(eu < 0., eu%(np.pi*np.array([2.,1.,2.])),eu)
@staticmethod @staticmethod