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

more tests 

Hotfix needed for axis angle to matrix (not used in DAMASK, needs
further investigation)
This commit is contained in:
Martin Diehl 2020-05-19 19:41:50 +02:00
parent de79a5af43
commit a6b0aaffba
3 changed files with 39 additions and 21 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
python/damask/_rotation.py
View File

@ -669,7 +669,7 @@ class Rotation:
"""Rotation matrix to Bunge-Euler angles.""" """Rotation matrix to Bunge-Euler angles."""
with np.errstate(invalid='ignore',divide='ignore'): with np.errstate(invalid='ignore',divide='ignore'):
zeta = 1.0/np.sqrt(1.0-om[...,2,2:3]**2) zeta = 1.0/np.sqrt(1.0-om[...,2,2:3]**2)
eu = np.where(np.isclose(np.abs(om[...,2,2:3]),1.0,1e-4), eu = np.where(np.isclose(np.abs(om[...,2,2:3]),1.0,1e-9),
np.block([np.arctan2(om[...,0,1:2],om[...,0,0:1]), np.block([np.arctan2(om[...,0,1:2],om[...,0,0:1]),
np.pi*0.5*(1-om[...,2,2:3]), np.pi*0.5*(1-om[...,2,2:3]),
np.zeros(om.shape[:-2]+(1,)), np.zeros(om.shape[:-2]+(1,)),
@ -679,30 +679,30 @@ class Rotation:
np.arctan2(om[...,0,2:3]*zeta,+om[...,1,2:3]*zeta) np.arctan2(om[...,0,2:3]*zeta,+om[...,1,2:3]*zeta)
]) ])
) )
eu[np.abs(eu)<1.e-6] = 0.0 eu[np.abs(eu)<1.e-8] = 0.0
eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu) 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 return eu
@staticmethod @staticmethod
def om2ax(om): def om2ax(om):
"""Rotation matrix to axis angle pair.""" """Rotation matrix to axis angle pair."""
return Rotation.qu2ax(Rotation.om2qu(om)) # HOTFIX
diag_delta = -_P*np.block([om[...,1,2:3]-om[...,2,1:2], diag_delta = -_P*np.block([om[...,1,2:3]-om[...,2,1:2],
om[...,2,0:1]-om[...,0,2:3], om[...,2,0:1]-om[...,0,2:3],
om[...,0,1:2]-om[...,1,0:1] om[...,0,1:2]-om[...,1,0:1]
]) ])
diag_delta[np.abs(diag_delta)<1.e-6] = 1.0
t = 0.5*(om.trace(axis2=-2,axis1=-1) -1.0).reshape(om.shape[:-2]+(1,)) t = 0.5*(om.trace(axis2=-2,axis1=-1) -1.0).reshape(om.shape[:-2]+(1,))
w,vr = np.linalg.eig(om) w,vr = np.linalg.eig(om)
# mask duplicated real eigenvalues # mask duplicated real eigenvalues
w[np.isclose(w[...,0],1.0+0.0j),1:] = 0. w[np.isclose(w[...,0],1.0+0.0j),1:] = 0.
w[np.isclose(w[...,1],1.0+0.0j),2:] = 0. w[np.isclose(w[...,1],1.0+0.0j),2:] = 0.
vr = np.swapaxes(vr,-1,-2) vr = np.swapaxes(vr,-1,-2)
ax = np.where(np.abs(diag_delta)<0, ax = np.where(np.abs(diag_delta)<1e-12,
np.real(vr[np.isclose(w,1.0+0.0j)]).reshape(om.shape[:-2]+(3,)), np.real(vr[np.isclose(w,1.0+0.0j)]).reshape(om.shape[:-2]+(3,)),
np.abs(np.real(vr[np.isclose(w,1.0+0.0j)]).reshape(om.shape[:-2]+(3,))) \ np.abs(np.real(vr[np.isclose(w,1.0+0.0j)]).reshape(om.shape[:-2]+(3,))) \
*np.sign(diag_delta)) *np.sign(diag_delta))
ax = np.block([ax,np.arccos(np.clip(t,-1.0,1.0))]) ax = np.block([ax,np.arccos(np.clip(t,-1.0,1.0))])
ax[np.abs(ax[...,3])<1.e-6] = [ 0.0, 0.0, 1.0, 0.0] ax[np.abs(ax[...,3])<1.e-8] = [ 0.0, 0.0, 1.0, 0.0]
return ax return ax
@staticmethod @staticmethod
@ -804,6 +804,7 @@ 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

11
python/tests/rotation_conversion.py
View File

@ -134,25 +134,26 @@ def om2qu(a):
def om2eu(om): def om2eu(om):
"""Rotation matrix to Bunge-Euler angles.""" """Rotation matrix to Bunge-Euler angles."""
if not np.isclose(np.abs(om[2,2]),1.0,1.e-4): if not np.isclose(np.abs(om[2,2]),1.0,1.e-9):
zeta = 1.0/np.sqrt(1.0-om[2,2]**2) zeta = 1.0/np.sqrt(1.0-om[2,2]**2)
eu = np.array([np.arctan2(om[2,0]*zeta,-om[2,1]*zeta), eu = np.array([np.arctan2(om[2,0]*zeta,-om[2,1]*zeta),
np.arccos(om[2,2]), np.arccos(om[2,2]),
np.arctan2(om[0,2]*zeta, om[1,2]*zeta)]) np.arctan2(om[0,2]*zeta, om[1,2]*zeta)])
else: else:
eu = np.array([np.arctan2( om[0,1],om[0,0]), np.pi*0.5*(1-om[2,2]),0.0]) # following the paper, not the reference implementation eu = np.array([np.arctan2( om[0,1],om[0,0]), np.pi*0.5*(1-om[2,2]),0.0]) # following the paper, not the reference implementation
eu[np.abs(eu)<1.e-6] = 0.0 eu[np.abs(eu)<1.e-8] = 0.0
eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu) 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 return eu
def om2ax(om): def om2ax(om):
"""Rotation matrix to axis angle pair.""" """Rotation matrix to axis angle pair."""
return qu2ax(om2qu(om)) # HOTFIX
ax=np.empty(4) ax=np.empty(4)
# first get the rotation angle # first get the rotation angle
t = 0.5*(om.trace() -1.0) t = 0.5*(om.trace() -1.0)
ax[3] = np.arccos(np.clip(t,-1.0,1.0)) ax[3] = np.arccos(np.clip(t,-1.0,1.0))
if np.abs(ax[3])<1.e-6: if np.abs(ax[3])<1.e-8:
ax = np.array([ 0.0, 0.0, 1.0, 0.0]) ax = np.array([ 0.0, 0.0, 1.0, 0.0])
else: else:
w,vr = np.linalg.eig(om) w,vr = np.linalg.eig(om)
@ -160,8 +161,7 @@ def om2ax(om):
i = np.where(np.isclose(w,1.0+0.0j))[0][0] i = np.where(np.isclose(w,1.0+0.0j))[0][0]
ax[0:3] = np.real(vr[0:3,i]) ax[0:3] = np.real(vr[0:3,i])
diagDelta = -_P*np.array([om[1,2]-om[2,1],om[2,0]-om[0,2],om[0,1]-om[1,0]]) diagDelta = -_P*np.array([om[1,2]-om[2,1],om[2,0]-om[0,2],om[0,1]-om[1,0]])
diagDelta[np.abs(diagDelta)<1.e-6] = 1.0 ax[0:3] = np.where(np.abs(diagDelta)<1e-12, ax[0:3],np.abs(ax[0:3])*np.sign(diagDelta))
ax[0:3] = np.where(np.abs(diagDelta)<0, ax[0:3],np.abs(ax[0:3])*np.sign(diagDelta))
return ax return ax
#---------- Bunge-Euler angles ---------- #---------- Bunge-Euler angles ----------
@ -229,6 +229,7 @@ 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

38
python/tests/test_Rotation.py
View File

@ -107,10 +107,10 @@ class TestRotation:
print(m,o,rot.as_quaternion()) print(m,o,rot.as_quaternion())
assert ok and np.isclose(np.linalg.norm(o),1.0) assert ok and np.isclose(np.linalg.norm(o),1.0)
@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):
@ -121,12 +121,12 @@ class TestRotation:
print(m,o,rot.as_quaternion()) print(m,o,rot.as_quaternion())
assert ok and np.isclose(np.linalg.det(o),1.0) assert ok and np.isclose(np.linalg.det(o),1.0)
@pytest.mark.parametrize('forward,backward',[(Rotation.eu2qu,Rotation.qu2eu)]) @pytest.mark.parametrize('forward,backward',[(Rotation.eu2qu,Rotation.qu2eu),
#(Rotation.eu2om,Rotation.om2eu), (Rotation.eu2om,Rotation.om2eu),
#(Rotation.eu2ax,Rotation.ax2eu), (Rotation.eu2ax,Rotation.ax2eu),
#(Rotation.eu2ro,Rotation.ro2eu), (Rotation.eu2ro,Rotation.ro2eu),
#(Rotation.eu2ho,Rotation.ho2eu), (Rotation.eu2ho,Rotation.ho2eu),
#(Rotation.eu2cu,Rotation.cu2eu)]) (Rotation.eu2cu,Rotation.cu2eu)])
def test_Eulers_internal(self,default,forward,backward): def test_Eulers_internal(self,default,forward,backward):
for rot in default: for rot in default:
m = rot.as_Eulers() m = rot.as_Eulers()
@ -140,6 +140,22 @@ class TestRotation:
print(m,o,rot.as_quaternion()) print(m,o,rot.as_quaternion())
assert ok and (np.zeros(3)-1.e-9 <= o).all() and (o <= np.array([np.pi*2.,np.pi,np.pi*2.])+1.e-9).all() assert ok and (np.zeros(3)-1.e-9 <= o).all() and (o <= np.array([np.pi*2.,np.pi,np.pi*2.])+1.e-9).all()
@pytest.mark.parametrize('forward,backward',[(Rotation.ax2qu,Rotation.qu2ax),
(Rotation.ax2om,Rotation.om2ax),
(Rotation.ax2eu,Rotation.eu2ax),
(Rotation.ax2ro,Rotation.ro2ax),
(Rotation.ax2ho,Rotation.ho2ax),
(Rotation.ax2cu,Rotation.cu2ax)])
def test_axis_angle_internal(self,default,forward,backward):
for rot in default:
m = rot.as_axis_angle()
o = backward(forward(m))
ok = np.allclose(m,o,atol=atol)
if np.isclose(m[3],np.pi,atol=atol):
ok = ok or np.allclose(m*np.array([-1.,-1.,-1.,1.]),o,atol=atol)
print(m,o,rot.as_quaternion())
assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) and o[3]<=np.pi+1.e-7
@pytest.mark.parametrize('degrees',[True,False]) @pytest.mark.parametrize('degrees',[True,False])
def test_Eulers(self,default,degrees): def test_Eulers(self,default,degrees):
@ -177,7 +193,7 @@ class TestRotation:
if np.isclose(m[3],np.pi,atol=atol): if np.isclose(m[3],np.pi,atol=atol):
ok = ok or np.allclose(m*np.array([-1.,-1.,-1.,1.]),o,atol=atol) ok = ok or np.allclose(m*np.array([-1.,-1.,-1.,1.]),o,atol=atol)
print(m,o,rot.as_quaternion()) print(m,o,rot.as_quaternion())
assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) and o[3]<=np.pi++1.e-9 assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) and o[3]<=np.pi+1.e-7
@pytest.mark.parametrize('P',[1,-1]) @pytest.mark.parametrize('P',[1,-1])
@pytest.mark.parametrize('normalise',[True,False]) @pytest.mark.parametrize('normalise',[True,False])