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

need to transpose eigenvectors to find the correct one

This commit is contained in:
Martin Diehl 2020-04-09 15:01:01 +02:00
parent e502573e05
commit bab3581b11
2 changed files with 13 additions and 11 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
python/damask/_rotation.py
View File

@ -631,7 +631,6 @@ class Rotation:
# 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-6:
ax = np.array([ 0.0, 0.0, 1.0, 0.0]) ax = np.array([ 0.0, 0.0, 1.0, 0.0])
else: else:
@ -639,22 +638,25 @@ class Rotation:
# next, find the eigenvalue (1,0j) # next, find the eigenvalue (1,0j)
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 = 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]])
ax[0:3] = np.where(np.abs(diagDelta)<0, ax[0:3],np.abs(ax[0:3])*np.sign(-P*diagDelta)) 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: else:
diag_delta = 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)
ax = np.where(np.abs(diag_delta)<0, 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.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(diag_delta)*np.sign(-P*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-6] = [ 0.0, 0.0, 1.0, 0.0]
return ax return ax

2
python/tests/test_Rotation.py
View File

@ -166,7 +166,7 @@ class TestRotation:
assert np.allclose(conversion(q),c) assert np.allclose(conversion(q),c)
@pytest.mark.parametrize('conversion',[Rotation.om2eu, @pytest.mark.parametrize('conversion',[Rotation.om2eu,
#Rotation.om2ax, Rotation.om2ax,
]) ])
def test_matrix_vectorization(self,default,conversion): def test_matrix_vectorization(self,default,conversion):
om = np.array([rot.asMatrix() for rot in default]) om = np.array([rot.asMatrix() for rot in default])