Merge branch '379-logic-of-to_frame-and-to_lattice-in-damask-orientat…

…ion' into 'development'

Lab-frame for Orientation.to_lattice/to_frame

Closes #379

See merge request damask/DAMASK!939

python/damask/_crystal.py

@@ -542,8 +542,10 @@ def __repr__(self):
  return family if self.lattice is None else \
  util.srepr([family,
  f'Bravais lattice: {self.lattice}',
- 'a={:.5g} m, b={:.5g} m, c={:.5g} m'.format(*self.parameters[:3]),
- 'α={:.5g}°, β={:.5g}°, γ={:.5g}°'.format(*np.degrees(self.parameters[3:]))])
+ 'a={a:.5g} m, b={b:.5g} m, c={c:.5g} m'.format(**self.parameters),
+ 'α={alpha:.5g}°, β={beta:.5g}°, γ={gamma:.5g}°'
+ .format(**dict(map(lambda kv: (kv[0], np.degrees(kv[1])), self.parameters.items()))),
+ ])
 
 
  def __eq__(self,
@@ -565,9 +567,10 @@ def __eq__(self,
  self.family == other.family) # type: ignore
 
  @property
- def parameters(self) -> Optional[Tuple]:
+ def parameters(self) -> Optional[Dict]:
  """Return lattice parameters a, b, c, alpha, beta, gamma."""
- return (self.a,self.b,self.c,self.alpha,self.beta,self.gamma) if hasattr(self,'a') else None
+ return dict(a=self.a,b=self.b,c=self.c,
+ alpha=self.alpha,beta=self.beta,gamma=self.gamma) if hasattr(self,'a') else None
 
  @property
  def immutable(self) -> Dict[str, float]:
@@ -849,9 +852,9 @@ def to_lattice(self, *,
  """
  if (direction is not None) ^ (plane is None):
  raise KeyError('specify either "direction" or "plane"')
- basis,axis = (self.basis_reciprocal,np.array(direction)) \
+ basis,axis = (self.basis_reciprocal,np.asarray(direction)) \
  if plane is None else \
- (self.basis_real,np.array(plane))
+ (self.basis_real,np.asarray(plane))
  return np.einsum('li,...l',basis,axis)
 
 
@@ -891,9 +894,9 @@ def to_frame(self, *,
  """
  if (uvw is not None) ^ (hkl is None):
  raise KeyError('specify either "uvw" or "hkl"')
- basis,axis = (self.basis_real,np.array(uvw)) \
+ basis,axis = (self.basis_real,np.asarray(uvw)) \
  if hkl is None else \
- (self.basis_reciprocal,np.array(hkl))
+ (self.basis_reciprocal,np.asarray(hkl))
  return np.einsum('il,...l',basis,axis)
 
 
@@ -1208,12 +1211,12 @@ def relation_operations(self,
 
  m_l,o_l = transform[0].split(sep) # type: ignore
  m_p,o_p = orientation_relationships[model][m_l+sep+o_l]
- other = Crystal(lattice=o_l) if target is None else target
- m_p = np.stack((self.to_frame(uvw=m_p[:,0] if len(m_p[0,0])==3 else util.Bravais_to_Miller(uvtw=m_p[:,0])),
-  self.to_frame(hkl=m_p[:,1] if len(m_p[0,1])==3 else util.Bravais_to_Miller(hkil=m_p[:,1]))),
- axis=1)
- o_p = np.stack((other.to_frame(uvw=o_p[:,0] if len(o_p[0,0])==3 else util.Bravais_to_Miller(uvtw=o_p[:,0])),
-  other.to_frame(hkl=o_p[:,1] if len(o_p[0,1])==3 else util.Bravais_to_Miller(hkil=o_p[:,1]))),
- axis=1)
-
+ m = Crystal(lattice=m_l) if self.parameters is None else Crystal(lattice=m_l,**self.parameters) # type: ignore
+ o = Crystal(lattice=o_l) if target is None else target
+ m_p = np.stack((m.to_frame(uvw=m_p[:,0] if len(m_p[0,0])==3 else util.Bravais_to_Miller(uvtw=m_p[:,0])),
+ m.to_frame(hkl=m_p[:,1] if len(m_p[0,1])==3 else util.Bravais_to_Miller(hkil=m_p[:,1]))),
+  axis=-2)
+ o_p = np.stack((o.to_frame(uvw=o_p[:,0] if len(o_p[0,0])==3 else util.Bravais_to_Miller(uvtw=o_p[:,0])),
+ o.to_frame(hkl=o_p[:,1] if len(o_p[0,1])==3 else util.Bravais_to_Miller(hkil=o_p[:,1]))),
+ axis=-2)
  return (o_l,Rotation.from_parallel(a=m_p,b=o_p))

python/damask/_orientation.py

@@ -57,7 +57,7 @@ class Orientation(Rotation,Crystal):
  and inherits the corresponding crystal family.
  Specifying a Bravais lattice, compared to just the crystal family,
  extends the functionality of Orientation objects to include operations such as
- "Schmid", "related", or "to_pole" that require a lattice type and its parameters.
+ "Schmid", "related", or "to_frame" that require a lattice type and its parameters.
 
  Examples
  --------
@@ -353,7 +353,7 @@ def from_directions(cls,
  new : damask.Orientation
 
  """
- o = cls(**kwargs)
+ o = cls(**kwargs,rotation=[1,0,0,0])
  x = o.to_frame(uvw=uvw)
  z = o.to_frame(hkl=hkl)
  om = np.stack([x,np.cross(z,x),z],axis=-2)
@@ -607,7 +607,7 @@ def to_SST(self,
  Parameters
  ----------
  vector : numpy.ndarray, shape (...,3)
- Lab frame vector to align with crystal frame direction.
+ Lab frame vector to align with an SST crystal frame direction.
  Shape of vector blends with shape of own rotation array.
  For example, a rotation array of shape (3,2) and a vector array of shape (2,4) result in (3,2,4) outputs.
  proper : bool, optional
@@ -629,8 +629,8 @@ def to_SST(self,
  if vector_.shape[-1] != 3:
  raise ValueError('input is not a field of three-dimensional vectors')
 
- blend = util.shapeblender( self.shape,vector_.shape[:-1])
- eq = self.broadcast_to(util.shapeshifter( self.shape,blend,mode='right')).equivalent
+ blend = util.shapeblender(self.shape,vector_.shape[:-1])
+ eq = self.broadcast_to(util.shapeshifter(self.shape,blend,mode='right')).equivalent
  poles = np.atleast_2d(eq @ np.broadcast_to(vector_,(1,)+blend+(3,)))
  ok = self.in_SST(poles,proper=proper)
  ok &= np.cumsum(ok,axis=0) == 1
@@ -775,11 +775,58 @@ def IPF_color(self,
 ####################################################################################################
  # functions that require lattice, not just family
 
- def to_pole(self, *,
- uvw: Optional[FloatSequence] = None,
- hkl: Optional[FloatSequence] = None,
- with_symmetry: bool = False,
- normalize: bool = True) -> np.ndarray:
+ def to_lattice(self, *,
+ direction: Optional[FloatSequence] = None,
+ plane: Optional[FloatSequence] = None) -> np.ndarray:
+ """
+ Calculate lattice vector corresponding to lab frame direction or plane normal.
+
+ Parameters
+ ----------
+ direction|plane : numpy.ndarray, shape (...,3)
+ Real space vector along direction or
+ reciprocal space vector along plane normal.
+ Shape of vector blends with shape of own rotation array.
+ For example, a rotation array of shape (3,2) and a vector
+ array of shape (2,4) result in (3,2,4) outputs.
+
+ Returns
+ -------
+ Miller : numpy.ndarray, shape (...,3)
+ Lattice vector of direction or plane.
+ Use util.scale_to_coprime to convert to (integer) Miller indices.
+
+ Examples
+ --------
+ >>> import np
+ >>> import damask
+ >>> np.set_printoptions(precision=2,suppress=True,floatmode='maxprec')
+ >>> cubic = damask.Orientation.from_axis_angle(n_omega=[1,0,0,90],degrees=True,lattice='cI')
+ >>> cubic.to_lattice(direction=[1, 0, 0])
+ array([1., 0., 0.])
+ >>> cubic.to_lattice(direction=[0, 1, 0])
+ array([ 0., 0., -1.])
+ >>> cubic.to_lattice(direction=[0, 0, 1])
+ array([-0., 1., 0.])
+ >>> tetragonal = damask.Orientation(lattice='tI',c=0.5)
+ >>> damask.util.scale_to_coprime(tetragonal.to_lattice(direction=[1,1,1]))
+ array([1, 1, 2])
+ >>> damask.util.scale_to_coprime(tetragonal.to_lattice(plane=[1,1,1]))
+ array([2, 2, 1])
+
+ """
+ if (direction is not None) ^ (plane is None):
+ raise KeyError('specify either "direction" or "plane"')
+ return (super().to_lattice(direction=self@np.asarray(direction)) if plane is None else
+ super().to_lattice(plane=self@np.asarray(plane)))
+
+
+ def to_frame(self, *,
+ uvw: Optional[FloatSequence] = None,
+ hkl: Optional[FloatSequence] = None,
+ with_symmetry: bool = False,
+ normalize: bool = True,
+ ) -> np.ndarray:
  """
  Calculate lab frame vector along lattice direction [uvw] or plane normal (hkl).
 
@@ -799,12 +846,12 @@ def to_pole(self, *,
 
  Returns
  -------
- vector : numpy.ndarray, shape (...,3) or (...,N,3)
- Lab frame vector (or vectors if with_symmetry) along
+ vector : numpy.ndarray, shape (...,3) or (N,...,3)
+ Lab frame vector (or N vectors if with_symmetry) along
  [uvw] direction or (hkl) plane normal.
 
  """
- v = self.to_frame(uvw=uvw,hkl=hkl)
+ v = super().to_frame(uvw=uvw,hkl=hkl)
  s_v = v.shape[:-1]
  blend = util.shapeblender(self.shape,s_v)
  if normalize:
@@ -815,7 +862,9 @@ def to_pole(self, *,
  blend += sym_ops.shape
  v = sym_ops.broadcast_to(s_v) @ v[...,np.newaxis,:]
 
- return ~(self.broadcast_to(blend)) @ np.broadcast_to(v,blend+(3,))
+ return np.moveaxis(~(self.broadcast_to(blend)) @ np.broadcast_to(v,blend+(3,)),
+ -2 if with_symmetry else 0,
+ 0)
 
 
  def Schmid(self, *,
@@ -859,8 +908,8 @@ def Schmid(self, *,
 
  if not active:
  raise ValueError('Schmid matrix not defined')
- d = self.to_frame(uvw=np.vstack([kinematics['direction'][i][:n] for i,n in enumerate(active)]))
- p = self.to_frame(hkl=np.vstack([kinematics['plane'][i][:n] for i,n in enumerate(active)]))
+ d = super().to_frame(uvw=np.vstack([kinematics['direction'][i][:n] for i,n in enumerate(active)]))
+ p = super().to_frame(hkl=np.vstack([kinematics['plane'][i][:n] for i,n in enumerate(active)]))
  P = np.einsum('...i,...j',d/np.linalg.norm(d,axis=1,keepdims=True),
  p/np.linalg.norm(p,axis=1,keepdims=True))
 
@@ -913,7 +962,6 @@ def related(self: MyType,
  """
  lattice,o = self.relation_operations(model,target)
  target = Crystal(lattice=lattice) if target is None else target
-
  return Orientation(rotation=o*Rotation(self.quaternion)[np.newaxis,...], # type: ignore
  lattice=lattice,
  b = self.b if target.ratio['b'] is None else self.a*target.ratio['b'],

python/damask/_result.py

@@ -1132,7 +1132,7 @@ def add_pole(self,
  Miller indices of crystallographic direction or plane normal.
  with_symmetry : bool, optional
  Calculate all N symmetrically equivalent vectors.
- Defaults to True.
+ Defaults to False.
  normalize : bool, optional
  Normalize output vector.
  Defaults to True.
@@ -1145,17 +1145,19 @@ def pole(q: DADF5Dataset,
  c = q['meta']['c/a'] if 'c/a' in q['meta'] else 1.0
  brackets = ['[]','()','⟨⟩','{}'][(uvw is None)*1+with_symmetry*2]
  label = 'p^' + '{}{} {} {}{}'.format(brackets[0],
-  *(uvw if uvw else hkl),
-  brackets[-1],)
+ *(uvw if uvw else hkl),
+ brackets[-1],)
  ori = Orientation(q['data'],lattice=q['meta']['lattice'],a=1,c=c)
 
  return {
- 'data': ori.to_pole(uvw=uvw,hkl=hkl,with_symmetry=with_symmetry,normalize=normalize),
+ 'data': np.moveaxis(ori.to_frame(uvw=uvw,hkl=hkl,
+ with_symmetry=with_symmetry,
+ normalize=normalize),0,-2 if with_symmetry else 0),
  'label': label,
  'meta' : {
  'unit': '1',
  'description': f'{"normalized " if normalize else ""}lab frame vector along lattice ' \
- + ('direction' if uvw is not None else 'plane') \
+ + ('plane' if uvw is None else 'direction') \
  + ('s' if with_symmetry else ''),
  'creator': 'add_pole'
  }

python/damask/_rotation.py

@@ -1002,7 +1002,7 @@ def from_parallel(a: np.ndarray,
  """
  a_ = np.array(a,dtype=float)
  b_ = np.array(b,dtype=float)
- if a_.shape[-2:] != (2,3) or b_.shape[-2:] != (2,3) or a_.shape != b_.shape:
+ if a_.shape[-2:] != (2,3) or b_.shape[-2:] != (2,3):
  raise ValueError('invalid shape')
  am = np.stack([ a_[...,0,:],
  a_[...,1,:],

python/tests/test_Crystal.py

@@ -46,6 +46,19 @@ def test_basis_invalid(self):
  with pytest.raises(KeyError):
  Crystal(family='cubic').basis_real
 
+ def test_basis_real(self):
+ for gamma in np.random.random(2**8)*np.pi:
+ basis = np.tril(np.random.random((3,3))+1e-6)
+ basis[1,:2] = basis[1,1]*np.array([np.cos(gamma),np.sin(gamma)])
+ basis[2,:2] = basis[2,:2]*2-1
+ lengths = np.linalg.norm(basis,axis=-1)
+ cosines = np.roll(np.einsum('ij,ij->i',basis,np.roll(basis,1,axis=0))/lengths/np.roll(lengths,1),1)
+ o = Crystal(lattice='aP',
+ **dict(zip(['a','b','c'],lengths)),
+ **dict(zip(['alpha','beta','gamma'],np.arccos(cosines))),
+ )
+ assert np.allclose(o.to_frame(uvw=np.eye(3)),basis,rtol=1e-4), 'Lattice basis disagrees with initialization'
+
  @pytest.mark.parametrize('keyFrame,keyLattice',[('uvw','direction'),('hkl','plane'),])
  @pytest.mark.parametrize('vector',np.array([
  [1.,1.,1.],

python/tests/test_Orientation.py

@@ -180,8 +180,8 @@ def test_from_directions(self,kwargs):
  c = np.cross(b,a)
  if np.allclose(c,0): continue
  o = Orientation.from_directions(uvw=a,hkl=c,**kwargs)
- x = o.to_pole(uvw=a)
- z = o.to_pole(hkl=c)
+ x = o.to_frame(uvw=a)
+ z = o.to_frame(hkl=c)
  assert np.isclose(np.dot(x,np.array([1,0,0])),1) \
  and np.isclose(np.dot(z,np.array([0,0,1])),1)
 
@@ -317,19 +317,6 @@ def test_relationship_reference(self,update,res_path,model,lattice):
  Table({'Eulers':(3,)},eu).set('pos',coords).save(reference)
  assert np.allclose(eu,Table.load(reference).get('Eulers'))
 
- def test_basis_real(self):
- for gamma in np.random.random(2**8)*np.pi:
- basis = np.tril(np.random.random((3,3))+1e-6)
- basis[1,:2] = basis[1,1]*np.array([np.cos(gamma),np.sin(gamma)])
- basis[2,:2] = basis[2,:2]*2-1
- lengths = np.linalg.norm(basis,axis=-1)
- cosines = np.roll(np.einsum('ij,ij->i',basis,np.roll(basis,1,axis=0))/lengths/np.roll(lengths,1),1)
- o = Orientation.from_random(lattice='aP',
- **dict(zip(['a','b','c'],lengths)),
- **dict(zip(['alpha','beta','gamma'],np.arccos(cosines))),
- )
- assert np.allclose(o.to_frame(uvw=np.eye(3)),basis,rtol=1e-4), 'Lattice basis disagrees with initialization'
-
  @pytest.mark.parametrize('lattice,a,b,c,alpha,beta,gamma',
  [
  ('aP',0.5,2.0,3.0,0.8,0.5,1.2),
@@ -347,14 +334,14 @@ def test_basis_real(self):
  np.random.random( (4,3)),
  np.random.random((4,8,3)),
  ])
- def test_to_pole(self,shape,lattice,a,b,c,alpha,beta,gamma,vector,kw,with_symmetry):
+ def test_to_frame(self,shape,lattice,a,b,c,alpha,beta,gamma,vector,kw,with_symmetry):
  o = Orientation.from_random(shape=shape,
  lattice=lattice,
  a=a,b=b,c=c,
  alpha=alpha,beta=beta,gamma=gamma)
- assert o.to_pole(**{kw:vector,'with_symmetry':with_symmetry}).shape \
- == util.shapeblender(o.shape,vector.shape[:-1]) \
- + (o.symmetry_operations.shape if with_symmetry else ()) \
+ assert o.to_frame(**{kw:vector,'with_symmetry':with_symmetry}).shape \
+ == (o.symmetry_operations.shape if with_symmetry else ()) \
+ + util.shapeblender(o.shape,vector.shape[:-1]) \
  + vector.shape[-1:]
 
  @pytest.mark.parametrize('lattice',['hP','cI','cF']) #tI not included yet
@@ -506,7 +493,7 @@ def test_to_SST_blending(self,family,left,right):
  ((3,1),(1,3)),
  (None,(3,)),
  ])
- def test_to_pole_blending(self,lattice,a,b,c,alpha,beta,gamma,left,right):
+ def test_to_frame_blending(self,lattice,a,b,c,alpha,beta,gamma,left,right):
  o = Orientation.from_random(shape=left,
  lattice=lattice,
  a=a,b=b,c=c,
@@ -516,8 +503,8 @@ def test_to_pole_blending(self,lattice,a,b,c,alpha,beta,gamma,left,right):
  for loc in np.random.randint(0,blend,(10,len(blend))):
  l = () if left is None else tuple(np.minimum(np.array(left )-1,loc[:len(left)]))
  r = () if right is None else tuple(np.minimum(np.array(right)-1,loc[-len(right):]))
- assert np.allclose(o[l].to_pole(uvw=v[r]),
- o.to_pole(uvw=v)[tuple(loc)])
+ assert np.allclose(o[l].to_frame(uvw=v[r]),
+ o.to_frame(uvw=v)[tuple(loc)])
 
  def test_mul_invalid(self):
  with pytest.raises(TypeError):
@@ -528,12 +515,11 @@ def test_mul_invalid(self):
  def test_OR_plot(self,update,res_path,tmp_path,OR,pole):
  # https://doi.org/10.3390/cryst13040663 for comparison
  O = Orientation(lattice='cF')
- poles = O.related(OR).to_pole(uvw=pole,with_symmetry=True).reshape(-1,3)
+ poles = O.related(OR).to_frame(uvw=pole,with_symmetry=True).reshape(-1,3)
  points = util.project_equal_area(poles,'z')
 
  fig, ax = plt.subplots()
- c = plt.Circle((0,0),1, color='k',fill=False)
- ax.add_patch(c)
+ ax.add_patch(plt.Circle((0,0),1, color='k',fill=False))
  ax.scatter(points[:,0],points[:,1])
  ax.set_aspect('equal', 'box')
  fname=f'{OR}-{"".join(map(str,pole))}.png'
@@ -543,4 +529,3 @@ def test_OR_plot(self,update,res_path,tmp_path,OR,pole):
  current = np.array(Image.open(tmp_path/fname))
  reference = np.array(Image.open(res_path/fname))
  assert np.allclose(current,reference)
-

python/tests/test_Result.py

@@ -280,7 +280,8 @@ def test_add_stress_second_Piola_Kirchhoff(self,default):
  def test_add_pole(self,default,options):
  default.add_pole(**options)
  rot = default.place('O')
- in_memory = Orientation(rot,lattice=rot.dtype.metadata['lattice']).to_pole(**options)
+ in_memory = np.moveaxis(Orientation(rot,lattice=rot.dtype.metadata['lattice']).to_frame(**options),
+ 0,-2 if options['with_symmetry'] else 0)
  brackets = [['[[]','[]]'],'()','⟨⟩','{}'][('hkl' in options)*1+(options['with_symmetry'])*2] # escape fnmatch
  label = 'p^{}{} {} {}{}'.format(brackets[0],
  *(list(options.values())[0]),