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

more tests and adjustments to tolerances

This commit is contained in:
Martin Diehl 2020-05-19 21:03:30 +02:00
parent a6b0aaffba
commit 1c53a37de4
3 changed files with 8 additions and 10 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

5
python/damask/_rotation.py
View File

@ -584,10 +584,10 @@ class Rotation:
with np.errstate(invalid='ignore',divide='ignore'): 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) 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)) 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],np.broadcast_to(np.pi,qu.shape[:-1]+(1,))]),
np.block([qu[...,1:4]*s,omega])) 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 return ax
@staticmethod @staticmethod
@ -804,7 +804,6 @@ class Rotation:
@staticmethod @staticmethod
def ax2om(ax): def ax2om(ax):
"""Axis angle pair to rotation matrix.""" """Axis angle pair to rotation matrix."""
return Rotation.qu2om(Rotation.ax2qu(ax)) # HOTFIX
c = np.cos(ax[...,3:4]) c = np.cos(ax[...,3:4])
s = np.sin(ax[...,3:4]) s = np.sin(ax[...,3:4])
omc = 1. -c omc = 1. -c

7
python/tests/rotation_conversion.py
View File

@ -81,9 +81,9 @@ def qu2ax(qu):
Modified version of the original formulation, should be numerically more stable 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 ]) 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) 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)) 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 ]) ax = ax = np.array([ qu[1]*s, qu[2]*s, qu[3]*s, omega ])
@ -229,7 +229,6 @@ def ax2qu(ax):
def ax2om(ax): def ax2om(ax):
"""Axis angle pair to rotation matrix.""" """Axis angle pair to rotation matrix."""
return qu2om(ax2qu(ax)) # HOTFIX
c = np.cos(ax[3]) c = np.cos(ax[3])
s = np.sin(ax[3]) s = np.sin(ax[3])
omc = 1.0-c omc = 1.0-c
@ -239,7 +238,7 @@ def ax2om(ax):
q = omc*ax[idx[0]] * ax[idx[1]] q = omc*ax[idx[0]] * ax[idx[1]]
om[idx[0],idx[1]] = q + s*ax[idx[2]] om[idx[0],idx[1]] = q + s*ax[idx[2]]
om[idx[1],idx[0]] = 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): def ax2ro(ax):
"""Axis angle pair to Rodrigues-Frank vector.""" """Axis angle pair to Rodrigues-Frank vector."""

6
python/tests/test_Rotation.py
View File

@ -110,9 +110,9 @@ class TestRotation:
@pytest.mark.parametrize('forward,backward',[(Rotation.om2qu,Rotation.qu2om), @pytest.mark.parametrize('forward,backward',[(Rotation.om2qu,Rotation.qu2om),
(Rotation.om2eu,Rotation.eu2om), (Rotation.om2eu,Rotation.eu2om),
(Rotation.om2ax,Rotation.ax2om), (Rotation.om2ax,Rotation.ax2om),
(Rotation.om2ro,Rotation.ro2om)]) (Rotation.om2ro,Rotation.ro2om),
#(Rotation.om2ho,Rotation.ho2om), (Rotation.om2ho,Rotation.ho2om),
#(Rotation.om2cu,Rotation.cu2om)]) (Rotation.om2cu,Rotation.cu2om)])
def test_matrix_internal(self,default,forward,backward): def test_matrix_internal(self,default,forward,backward):
for rot in default: for rot in default:
m = rot.as_matrix() m = rot.as_matrix()