diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index ee497e72b..0f0a081bc 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -584,10 +584,10 @@ class Rotation: with np.errstate(invalid='ignore',divide='ignore'): s = np.sign(qu[...,0:1])/np.sqrt(qu[...,1:2]**2+qu[...,2:3]**2+qu[...,3:4]**2) omega = 2.0 * np.arccos(np.clip(qu[...,0:1],-1.0,1.0)) - ax = np.where(np.broadcast_to(qu[...,0:1] < 1.0e-6,qu.shape), + ax = np.where(np.broadcast_to(qu[...,0:1] < 1.0e-8,qu.shape), np.block([qu[...,1:4],np.broadcast_to(np.pi,qu.shape[:-1]+(1,))]), np.block([qu[...,1:4]*s,omega])) - ax[np.sum(np.abs(qu[...,1:4])**2,axis=-1) < 1.0e-6,] = [0.0, 0.0, 1.0, 0.0] + ax[np.isclose(qu[...,0],1.,rtol=0.0)] = [0.0, 0.0, 1.0, 0.0] return ax @staticmethod @@ -804,7 +804,6 @@ class Rotation: @staticmethod def ax2om(ax): """Axis angle pair to rotation matrix.""" - return Rotation.qu2om(Rotation.ax2qu(ax)) # HOTFIX c = np.cos(ax[...,3:4]) s = np.sin(ax[...,3:4]) omc = 1. -c diff --git a/python/tests/rotation_conversion.py b/python/tests/rotation_conversion.py index ab503f402..3aca1a212 100644 --- a/python/tests/rotation_conversion.py +++ b/python/tests/rotation_conversion.py @@ -81,9 +81,9 @@ def qu2ax(qu): Modified version of the original formulation, should be numerically more stable """ - if np.abs(np.sum(qu[1:4]**2)) < 1.e-6: # set axis to [001] if the angle is 0/360 + if np.isclose(qu[0],1.,rtol=0.0): # set axis to [001] if the angle is 0/360 ax = np.array([ 0.0, 0.0, 1.0, 0.0 ]) - elif qu[0] > 1.e-6: + elif qu[0] > 1.e-8: s = np.sign(qu[0])/np.sqrt(qu[1]**2+qu[2]**2+qu[3]**2) omega = 2.0 * np.arccos(np.clip(qu[0],-1.0,1.0)) ax = ax = np.array([ qu[1]*s, qu[2]*s, qu[3]*s, omega ]) @@ -229,7 +229,6 @@ def ax2qu(ax): def ax2om(ax): """Axis angle pair to rotation matrix.""" - return qu2om(ax2qu(ax)) # HOTFIX c = np.cos(ax[3]) s = np.sin(ax[3]) omc = 1.0-c @@ -239,7 +238,7 @@ def ax2om(ax): q = omc*ax[idx[0]] * ax[idx[1]] om[idx[0],idx[1]] = q + s*ax[idx[2]] om[idx[1],idx[0]] = q - s*ax[idx[2]] - return om if _P < 0.0 else np.swapaxes(om,(-1,-2)) + return om if _P < 0.0 else om.T def ax2ro(ax): """Axis angle pair to Rodrigues-Frank vector.""" diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index adf507e49..f9565a5e2 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -110,9 +110,9 @@ class TestRotation: @pytest.mark.parametrize('forward,backward',[(Rotation.om2qu,Rotation.qu2om), (Rotation.om2eu,Rotation.eu2om), (Rotation.om2ax,Rotation.ax2om), - (Rotation.om2ro,Rotation.ro2om)]) - #(Rotation.om2ho,Rotation.ho2om), - #(Rotation.om2cu,Rotation.cu2om)]) + (Rotation.om2ro,Rotation.ro2om), + (Rotation.om2ho,Rotation.ho2om), + (Rotation.om2cu,Rotation.cu2om)]) def test_matrix_internal(self,default,forward,backward): for rot in default: m = rot.as_matrix()