From bab3581b1194c74e5ca3c217db8d2ad1e78b0531 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 9 Apr 2020 15:01:01 +0200 Subject: [PATCH] need to transpose eigenvectors to find the correct one --- python/damask/_rotation.py | 22 ++++++++++++---------- python/tests/test_Rotation.py | 2 +- 2 files changed, 13 insertions(+), 11 deletions(-) diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 4756e94c2..3336eb6c0 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -631,7 +631,6 @@ class Rotation: # first get the rotation angle t = 0.5*(om.trace() -1.0) ax[3] = np.arccos(np.clip(t,-1.0,1.0)) - if np.abs(ax[3])<1.e-6: ax = np.array([ 0.0, 0.0, 1.0, 0.0]) else: @@ -639,22 +638,25 @@ class Rotation: # next, find the eigenvalue (1,0j) i = np.where(np.isclose(w,1.0+0.0j))[0][0] ax[0:3] = np.real(vr[0:3,i]) - diagDelta = np.array([om[1,2]-om[2,1],om[2,0]-om[0,2],om[0,1]-om[1,0]]) - ax[0:3] = np.where(np.abs(diagDelta)<0, ax[0:3],np.abs(ax[0:3])*np.sign(-P*diagDelta)) + 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)<0, ax[0:3],np.abs(ax[0:3])*np.sign(diagDelta)) else: - diag_delta = np.block([om[...,1,2:3]-om[...,2,1:2], - om[...,2,0:1]-om[...,0,2:3], - om[...,0,1:2]-om[...,1,0:1] - ]) + diag_delta = -P*np.block([om[...,1,2:3]-om[...,2,1:2], + om[...,2,0:1]-om[...,0,2:3], + 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,)) w,vr = np.linalg.eig(om) # mask duplicated real eigenvalues w[np.isclose(w[...,0],1.0+0.0j),1:] = 0. w[np.isclose(w[...,1],1.0+0.0j),2:] = 0. + vr = np.swapaxes(vr,-1,-2) ax = np.where(np.abs(diag_delta)<0, - 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(diag_delta)*np.sign(-P*diag_delta)) + 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)) 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] return ax diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index 6bb0f5160..02d561b66 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -166,7 +166,7 @@ class TestRotation: assert np.allclose(conversion(q),c) @pytest.mark.parametrize('conversion',[Rotation.om2eu, - #Rotation.om2ax, + Rotation.om2ax, ]) def test_matrix_vectorization(self,default,conversion): om = np.array([rot.asMatrix() for rot in default])