Merge branch '434-fast-disorientation-calculation' into 'development'

fast calculation of disorientation angle

Closes #434

See merge request damask/DAMASK!973

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<m3)
+ mask3 = np.logical_not(mask1)
+
+ m2[mask2] = m3[mask2]
+ m2[mask3] = np.where(m4[mask3]<m1[mask3],m1[mask3],m4[mask3])
+ m1[mask3] = m3[mask3]
+
+ trace_max = np.maximum.reduce([trace_max,m1,(m1+m2)*np.sqrt(2.)/2.])
+
+ else:
+ return self.disorientation(other).as_axis_angle(pair=True)[1] # type: ignore
+
+ return 2.*np.arccos(np.clip(np.round(trace_max[...,0],15),None,1.))
+
+
  def average(self: MyType, # type: ignore[override]
  weights: Optional[FloatSequence] = None,
  return_cloud: bool = False) -> 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),

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)
 

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('shapes',[[None,None],
+ [[2,3,4],[2,3,4]],
+ [[3,4],[4,3]],
+ [1000,1000]])
+ def test_disorientation_angle(self,family,shapes):
+ o_1 = Orientation.from_random(shape=shapes[0],family=family)
+ o_2 = Orientation.from_random(shape=shapes[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('shapes',[[None,None,()],
  [[2,3,4],[2,3,4],(2,3,4)],
  [[3,4],[4,5],(3,4,5)],
@@ -267,16 +279,25 @@ def test_shape_blending(self,shapes):
  me,other,blend = shapes
  o_1 = Orientation.from_random(shape=me,family='triclinic')
  o_2 = Orientation.from_random(shape=other,family='triclinic')
+ angle = o_1.misorientation_angle(o_2)
  full = o_1.misorientation(o_2)
  composition = o_1*o_2
- assert full.shape == composition.shape == blend
+ assert angle.shape == full.shape == composition.shape == blend
 
  def test_disorientation_invalid(self):
  a,b = np.random.choice(list(crystal_families),2,False)
  o_1 = Orientation.from_random(family=a)
  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]},

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('shapes',[[None,None],
+ [[2,3,4],[2,3,4]],
+ [[3,4],[4,3]],
+ [1000,1000]])
+ def test_misorientation_angle(self,shapes):
+ r_1 = Rotation.from_random(shape=shapes[0])
+ r_2 = Rotation.from_random(shape=shapes[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
@@ -1121,9 +1137,10 @@ def test_shape_blending(self,shapes):
  me,other,blend = shapes
  r_1 = Rotation.from_random(shape=me)
  r_2 = Rotation.from_random(shape=other)
+ angle = r_1.misorientation_angle(r_2)
  full = r_1.misorientation(r_2)
  composition = r_1*r_2
- assert full.shape == composition.shape == blend
+ assert angle.shape == full.shape == composition.shape == blend
 
  def test_composition_inverse(self):
  a,b = (Rotation.from_random(),Rotation.from_random())