From 26251a0d1383c3e1f928091a151a1823d7b26399 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 27 Jul 2024 08:10:44 +0200 Subject: [PATCH] fast calculation of disorientation angle MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit code is from Lionel Germain, Université de Lorraine. Publication is under preparation --- python/damask/_orientation.py | 89 +++++++++++++++++++++++++++++++- python/damask/_rotation.py | 21 ++++++++ python/tests/test_Orientation.py | 25 ++++++++- python/tests/test_Rotation.py | 19 ++++++- 4 files changed, 149 insertions(+), 5 deletions(-) diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index 9deb57b3b..56878c621 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -550,6 +550,91 @@ def disorientation(self: MyType, ) + def disorientation_angle(self: MyType, + other: MyType) -> np.ndarray: + """ + Calculate disorientation angle between self and given other orientation. + + Parameters + ---------- + other : Orientation + Orientation to which the disorientation angle is computed. + Compatible innermost dimensions will blend. + + Returns + ------- + omega : np.ndarray + Disorientation angle. + + Notes + ----- + Requires same crystal family for both orientations. + + References + ---------- + Lionel Germain, personal communication. + + """ + q_abs = np.abs((self*~other).quaternion) + + if 'triclinic' == other.family == self.family: + trace_max = q_abs[...,0:1] + + elif 'monoclinic' == other.family == self.family: + trace_max = np.maximum(q_abs[...,0:1], + q_abs[...,2:3]) + + elif 'orthorhombic' == other.family == self.family: + trace_max = np.maximum.reduce([q_abs[...,0:1], + q_abs[...,1:2], + q_abs[...,2:3], + q_abs[...,3:4]]) + + elif 'tetragonal' == other.family == self.family: + m1,m2,m3,m4 = np.split(q_abs,4,axis=-1) + + trace_max = np.maximum.reduce([m1,m2,m3,m4, + (m1+m4)*np.sqrt(2.)/2., + (m2+m3)*np.sqrt(2.)/2.]) + + elif 'hexagonal' == other.family == self.family: + m1,m2,m3,m4 = np.split(q_abs,4,axis=-1) + + mask = m1 < m4 + m1[mask],m4[mask] = m4[mask],m1[mask] + mask = m2 < m3 + m2[mask],m3[mask] = m3[mask],m2[mask] + + trace_max = np.maximum.reduce([m1,m2, + m1*np.sqrt(3.)/2.+m4*.5, + m2*np.sqrt(3.)/2.+m3*.5]) + + elif 'cubic' == other.family == self.family: + m1,m2,m3,m4 = np.split(q_abs,4,axis=-1) + + trace_max = np.sum(q_abs,axis=-1,keepdims=True)*.5 + + mask = m1 < m2 + m1[mask],m2[mask] = m2[mask],m1[mask] + mask = m3 < m4 + m3[mask],m4[mask] = m4[mask],m3[mask] + + mask1 = m1 > m3 + mask2 = np.logical_and(mask1,m2 Union[Tuple[MyType, MyType], MyType]: @@ -579,8 +664,8 @@ def average(self: MyType, """ eq = self.equivalent - m = eq.misorientation(self[...,0].reshape((1,)+self.shape[:-1]+(1,)) # type: ignore - .broadcast_to(eq.shape)).as_axis_angle()[...,3] # type: ignore + m = eq.misorientation_angle(self[...,0].reshape((1,)+self.shape[:-1]+(1,)) + .broadcast_to(eq.shape)) r = Rotation(np.squeeze(np.take_along_axis(eq.quaternion, np.argmin(m,axis=0)[np.newaxis,...,np.newaxis], axis=0), diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index f16bf2706..835b7c4cb 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -637,6 +637,27 @@ def misorientation(self: MyType, return ~(self*~other) + def misorientation_angle(self: MyType, + other: MyType) -> np.ndarray: + """ + Calculate misorientation angle to other Rotation. + + Parameters + ---------- + other : damask.Rotation + Rotation to which the misorientation angle is computed. + Compatible innermost dimensions will blend. + + Returns + ------- + omega : np.ndarray + Misorientation angle. + + """ + trace_max = np.abs((self*~other).quaternion[...,0]) + return 2.*np.arccos(np.clip(np.round(trace_max,15),None,1.)) + + ################################################################################################ # convert to different orientation representations (numpy arrays) diff --git a/python/tests/test_Orientation.py b/python/tests/test_Orientation.py index 088a90281..932a6664b 100644 --- a/python/tests/test_Orientation.py +++ b/python/tests/test_Orientation.py @@ -257,6 +257,18 @@ def test_disorientation360(self,family): o_2 = Orientation.from_Euler_angles(family=family,phi=[360,0,0],degrees=True) assert np.allclose((o_1.disorientation(o_2)).as_matrix(),np.eye(3)) + @pytest.mark.parametrize('family',crystal_families) + @pytest.mark.parametrize('shape',[[None,None], + [[2,3,4],[2,3,4]], + [[3,4],[4,3]], + [1000,1000]]) + def test_disorientation_angle(self,family,shape): + o_1 = Orientation.from_random(shape=shape[0],family=family) + o_2 = Orientation.from_random(shape=shape[1],family=family) + angle = o_1.disorientation_angle(o_2) + full = o_1.disorientation(o_2).as_axis_angle(pair=True)[1] + assert np.allclose(angle,full,atol=1e-13,rtol=0) + @pytest.mark.parametrize('shape',[[None,None,()], [[2,3,4],[2,3,4],(2,3,4)], [[3,4],[4,5],(3,4,5)], @@ -266,9 +278,10 @@ def test_disorientation360(self,family): def test_shape_blending(self,shape): o_1 = Orientation.from_random(shape=shape[0],family='triclinic') o_2 = Orientation.from_random(shape=shape[1],family='triclinic') + angle = o_1.misorientation_angle(o_2) full = o_1.misorientation(o_2).as_axis_angle(pair=True)[1] composition = o_1*o_2 - assert full.shape == composition.shape == shape[2] + assert angle.shape == full.shape == composition.shape == shape[2] def test_disorientation_invalid(self): a,b = np.random.choice(list(crystal_families),2,False) @@ -276,6 +289,14 @@ def test_disorientation_invalid(self): o_2 = Orientation.from_random(family=b) with pytest.raises(NotImplementedError): o_1.disorientation(o_2) + with pytest.raises(NotImplementedError): + o_1.disorientation_angle(o_2) + + @pytest.mark.parametrize('family',crystal_families) + def test_disorientation_zero(self,set_of_quaternions,family): + o = Orientation.from_quaternion(q=set_of_quaternions,family=family) + assert np.allclose(o.disorientation_angle(o),0.0,atol=1e-15,rtol=0.) + assert np.allclose(o.disorientation(o).as_axis_angle(pair=True)[1],0.,atol=1e-15,rtol=0.) @pytest.mark.parametrize('color',[{'label':'red', 'RGB':[1,0,0],'direction':[0,0,1]}, {'label':'green','RGB':[0,1,0],'direction':[0,1,1]}, @@ -385,7 +406,7 @@ def test_to_frame_symmetries(self,lattice,mode,vector,N_sym): o = Orientation.from_random(lattice=lattice) frame = o.to_frame(**{keyword:vector,'with_symmetry':True}) shape_full = frame.shape[0] - shape_reduced = np.unique(np.around(frame,12),axis=0).shape[0] + shape_reduced = np.unique(np.around(frame,11),axis=0).shape[0] assert shape_full//N_sym == shape_reduced diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index f37929a8b..f04375a64 100644 --- a/python/tests/test_Rotation.py +++ b/python/tests/test_Rotation.py @@ -1101,6 +1101,22 @@ def test_misorientation_360deg(self): R_2 = Rotation.from_Euler_angles([360,0,0],degrees=True) assert np.allclose(R_1.misorientation(R_2).as_matrix(),np.eye(3)) + def test_misorientation_zero(self,set_of_quaternions): + r = Rotation.from_quaternion(set_of_quaternions) + assert np.allclose(r.misorientation_angle(r),0.0,atol=1e-15,rtol=0.) + assert np.allclose(r.misorientation(r).as_axis_angle(pair=True)[1],0.,atol=1e-15,rtol=0.) + + @pytest.mark.parametrize('shape',[[None,None], + [[2,3,4],[2,3,4]], + [[3,4],[4,3]], + [1000,1000]]) + def test_misorientation_angle(self,shape): + r_1 = Rotation.from_random(shape=shape[0]) + r_2 = Rotation.from_random(shape=shape[1]) + angle = r_1.misorientation_angle(r_2) + full = r_1.misorientation(r_2).as_axis_angle(pair=True)[1] + assert np.allclose(angle,full,atol=1e-13,rtol=0) + def test_composition(self): a,b = (Rotation.from_random(),Rotation.from_random()) c = a * b @@ -1120,9 +1136,10 @@ def test_composition_invalid(self): def test_shape_blending(self,shape): r_1 = Rotation.from_random(shape=shape[0]) r_2 = Rotation.from_random(shape=shape[1]) + angle = r_1.misorientation_angle(r_2) full = r_1.misorientation(r_2).as_axis_angle(pair=True)[1] composition = r_1*r_2 - assert full.shape == composition.shape == shape[2] + assert angle.shape == full.shape == composition.shape == shape[2] def test_composition_inverse(self): a,b = (Rotation.from_random(),Rotation.from_random())