overall smart blending; more robust disFZ
This commit is contained in:
Philip Eisenlohr
parent
05e675bbdb
commit
550c757cdc
|
|
@ -434,28 +434,29 @@ class Orientation(Rotation,Crystal):
|
|||
https://doi.org/10.1107/S0108767391006864
|
||||
|
||||
"""
|
||||
rho = self.as_Rodrigues_vector(compact=True)*(1.0-1.0e-9)
|
||||
|
||||
with np.errstate(invalid='ignore'):
|
||||
rho = self.as_Rodrigues_vector(compact=True)
|
||||
rho_ = rho + np.linalg.norm(rho,axis=-1,keepdims=True)*1e-9
|
||||
|
||||
if self.family == 'cubic':
|
||||
return ((rho[...,0] >= rho[...,1]) &
|
||||
(rho[...,1] >= rho[...,2]) &
|
||||
(rho[...,2] >= 0)).astype(bool)
|
||||
return ((rho_[...,0] >= rho[...,1]) &
|
||||
(rho_[...,1] >= rho[...,2]) &
|
||||
(rho_[...,2] >= 0)).astype(bool)
|
||||
if self.family == 'hexagonal':
|
||||
return ((rho[...,0] >= rho[...,1]*np.sqrt(3)) &
|
||||
(rho[...,1] >= 0) &
|
||||
(rho[...,2] >= 0)).astype(bool)
|
||||
return ((rho_[...,0] >= rho[...,1]*np.sqrt(3)) &
|
||||
(rho_[...,1] >= 0) &
|
||||
(rho_[...,2] >= 0)).astype(bool)
|
||||
if self.family == 'tetragonal':
|
||||
return ((rho[...,0] >= rho[...,1]) &
|
||||
(rho[...,1] >= 0) &
|
||||
(rho[...,2] >= 0)).astype(bool)
|
||||
return ((rho_[...,0] >= rho[...,1]) &
|
||||
(rho_[...,1] >= 0) &
|
||||
(rho_[...,2] >= 0)).astype(bool)
|
||||
if self.family == 'orthorhombic':
|
||||
return ((rho[...,0] >= 0) &
|
||||
(rho[...,1] >= 0) &
|
||||
(rho[...,2] >= 0)).astype(bool)
|
||||
return ((rho_[...,0] >= 0) &
|
||||
(rho_[...,1] >= 0) &
|
||||
(rho_[...,2] >= 0)).astype(bool)
|
||||
if self.family == 'monoclinic':
|
||||
return ((rho[...,1] >= 0) &
|
||||
(rho[...,2] >= 0)).astype(bool)
|
||||
return ((rho_[...,1] >= 0) &
|
||||
(rho_[...,2] >= 0)).astype(bool)
|
||||
|
||||
return np.ones_like(rho[...,0],dtype=bool)
|
||||
|
||||
|
|
@ -521,20 +522,21 @@ class Orientation(Rotation,Crystal):
|
|||
if self.family != other.family:
|
||||
raise NotImplementedError('disorientation between different crystal families')
|
||||
|
||||
blend = util.shapeblender(self.shape,other.shape)
|
||||
s = self.equivalent
|
||||
o = other.equivalent
|
||||
blend = util.shapeblender( self.shape,other.shape)
|
||||
s_m = util.shapeshifter( self.shape,blend,mode='right')
|
||||
s_o = util.shapeshifter(other.shape,blend,mode='left')
|
||||
|
||||
s_ = s.reshape((s.shape[0],1)+ self.shape).broadcast_to((s.shape[0],o.shape[0])+blend,mode='right')
|
||||
o_ = o.reshape((1,o.shape[0])+other.shape).broadcast_to((s.shape[0],o.shape[0])+blend,mode='right')
|
||||
r_ = s_.misorientation(o_)
|
||||
s = self.broadcast_to(s_m).equivalent
|
||||
o = other.broadcast_to(s_o).equivalent
|
||||
|
||||
r_ = s[:,np.newaxis,...].misorientation(o[np.newaxis,:,...])
|
||||
_r = ~r_
|
||||
|
||||
forward = r_.in_FZ & r_.in_disorientation_FZ
|
||||
reverse = _r.in_FZ & _r.in_disorientation_FZ
|
||||
ok = forward | reverse
|
||||
ok &= (np.cumsum(ok.reshape((-1,)+ok.shape[2:]),axis=0) == 1).reshape(ok.shape)
|
||||
r = np.where(np.any(forward[...,np.newaxis],axis=(0,1),keepdims=True),
|
||||
forward = r_.in_disorientation_FZ
|
||||
reverse = _r.in_disorientation_FZ
|
||||
ok = (forward | reverse) & r_.in_FZ
|
||||
ok &= (np.cumsum(ok.reshape((-1,)+blend),axis=0) == 1).reshape(ok.shape)
|
||||
r = np.where(np.any((ok&forward)[...,np.newaxis],axis=(0,1),keepdims=True),
|
||||
r_.quaternion,
|
||||
_r.quaternion)
|
||||
loc = np.where(ok)
|
||||
|
|
@ -584,6 +586,7 @@ class Orientation(Rotation,Crystal):
|
|||
np.argmin(m,axis=0)[np.newaxis,...,np.newaxis],
|
||||
axis=0),
|
||||
axis=0))
|
||||
|
||||
return ((self.copy(Rotation(r).average(weights)),self.copy(Rotation(r))) if return_cloud else
|
||||
self.copy(Rotation(r).average(weights))
|
||||
)
|
||||
|
|
@ -620,17 +623,19 @@ class Orientation(Rotation,Crystal):
|
|||
vector_ = np.array(vector,float)
|
||||
if vector_.shape[-1] != 3:
|
||||
raise ValueError('input is not a field of three-dimensional vectors')
|
||||
eq = self.equivalent
|
||||
blend = util.shapeblender(eq.shape,vector_.shape[:-1])
|
||||
poles = eq.broadcast_to(blend,mode='right') @ np.broadcast_to(vector_,blend+(3,))
|
||||
|
||||
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
|
||||
loc = np.where(ok)
|
||||
sort = 0 if len(loc) == 1 else np.lexsort(loc[:0:-1])
|
||||
|
||||
return (
|
||||
(poles[ok][sort].reshape(blend[1:]+(3,)), (np.vstack(loc[:1]).T)[sort].reshape(blend[1:]))
|
||||
(poles[ok][sort].reshape(blend+(3,)), (np.vstack(loc[:1]).T)[sort].reshape(blend))
|
||||
if return_operators else
|
||||
poles[ok][sort].reshape(blend[1:]+(3,))
|
||||
poles[ok][sort].reshape(blend+(3,))
|
||||
)
|
||||
|
||||
|
||||
|
|
@ -795,16 +800,17 @@ class Orientation(Rotation,Crystal):
|
|||
|
||||
"""
|
||||
v = self.to_frame(uvw=uvw,hkl=hkl)
|
||||
blend = util.shapeblender(self.shape,v.shape[:-1])
|
||||
s_v = v.shape[:-1]
|
||||
blend = util.shapeblender(self.shape,s_v)
|
||||
if normalize:
|
||||
v /= np.linalg.norm(v,axis=-1,keepdims=len(v.shape)>1)
|
||||
v /= np.linalg.norm(v,axis=-1,keepdims=len(s_v)>0)
|
||||
if with_symmetry:
|
||||
sym_ops = self.symmetry_operations
|
||||
shape = v.shape[:-1]+sym_ops.shape
|
||||
s_v += sym_ops.shape
|
||||
blend += sym_ops.shape
|
||||
v = sym_ops.broadcast_to(shape) \
|
||||
@ np.broadcast_to(v.reshape(util.shapeshifter(v.shape,shape+(3,))),shape+(3,))
|
||||
return ~(self.broadcast_to(blend))@np.broadcast_to(v,blend+(3,))
|
||||
v = sym_ops.broadcast_to(s_v) @ v[...,np.newaxis,:]
|
||||
|
||||
return ~(self.broadcast_to(blend)) @ np.broadcast_to(v,blend+(3,))
|
||||
|
||||
|
||||
def Schmid(self, *,
|
||||
|
|
@ -854,6 +860,7 @@ class Orientation(Rotation,Crystal):
|
|||
p/np.linalg.norm(p,axis=1,keepdims=True))
|
||||
|
||||
shape = P.shape[0:1]+self.shape+(3,3)
|
||||
|
||||
return ~self.broadcast_to(shape[:-2]) \
|
||||
@ np.broadcast_to(P.reshape(util.shapeshifter(P.shape,shape)),shape)
|
||||
|
||||
|
|
@ -897,6 +904,7 @@ class Orientation(Rotation,Crystal):
|
|||
lattice,o = self.relation_operations(model)
|
||||
target = Crystal(lattice=lattice)
|
||||
o = o.broadcast_to(o.shape+self.shape,mode='right')
|
||||
|
||||
return Orientation(rotation=o*Rotation(self.quaternion).broadcast_to(o.shape,mode='left'),
|
||||
lattice=lattice,
|
||||
b = self.b if target.ratio['b'] is None else self.a*target.ratio['b'],
|
||||
|
|
|
|||
|
|
@ -295,6 +295,7 @@ class Rotation:
|
|||
----------
|
||||
other : Rotation, shape (self.shape)
|
||||
Rotation for composition.
|
||||
Compatible innermost dimensions will blend.
|
||||
|
||||
Returns
|
||||
-------
|
||||
|
|
@ -303,10 +304,15 @@ class Rotation:
|
|||
|
||||
"""
|
||||
if isinstance(other,Rotation):
|
||||
q_m = self.quaternion[...,0:1]
|
||||
p_m = self.quaternion[...,1:]
|
||||
q_o = other.quaternion[...,0:1]
|
||||
p_o = other.quaternion[...,1:]
|
||||
blend = util.shapeblender( self.shape,other.shape)
|
||||
s_m = util.shapeshifter( self.shape,blend,mode='right')
|
||||
s_o = util.shapeshifter(other.shape,blend,mode='left')
|
||||
|
||||
q_m = self.broadcast_to(s_m).quaternion[...,0:1]
|
||||
p_m = self.broadcast_to(s_m).quaternion[...,1:]
|
||||
q_o = other.broadcast_to(s_o).quaternion[...,0:1]
|
||||
p_o = other.broadcast_to(s_o).quaternion[...,1:]
|
||||
|
||||
qmo = q_m*q_o
|
||||
q = (qmo - np.einsum('...i,...i',p_m,p_o).reshape(qmo.shape))
|
||||
p = q_m*p_o + q_o*p_m + _P * np.cross(p_m,p_o)
|
||||
|
|
@ -325,6 +331,7 @@ class Rotation:
|
|||
----------
|
||||
other : Rotation, shape (self.shape)
|
||||
Rotation for composition.
|
||||
Compatible innermost dimensions will blend.
|
||||
|
||||
"""
|
||||
return self*other
|
||||
|
|
@ -341,6 +348,7 @@ class Rotation:
|
|||
----------
|
||||
other : damask.Rotation, shape (self.shape)
|
||||
Rotation to invert for composition.
|
||||
Compatible innermost dimensions will blend.
|
||||
|
||||
Returns
|
||||
-------
|
||||
|
|
@ -434,9 +442,10 @@ class Rotation:
|
|||
|
||||
"""
|
||||
if isinstance(other, np.ndarray):
|
||||
obs = util.shapeblender(self.shape,other.shape,keep_ones=False)[len(self.shape):]
|
||||
obs = util.shapeblender(self.shape,other.shape)[len(self.shape):]
|
||||
for l in [4,2,1]:
|
||||
if obs[-l:] == l*(3,):
|
||||
print(f'rotate {l}')
|
||||
bs = util.shapeblender(self.shape,other.shape[:-l],False)
|
||||
self_ = self.broadcast_to(bs) if self.shape != bs else self
|
||||
if l==1:
|
||||
|
|
|
|||
|
|
@ -513,7 +513,7 @@ def shapeshifter(fro: _Tuple[int, ...],
|
|||
|
||||
def shapeblender(a: _Tuple[int, ...],
|
||||
b: _Tuple[int, ...],
|
||||
keep_ones: bool = True) -> _Tuple[int, ...]:
|
||||
keep_ones: bool = False) -> _Tuple[int, ...]:
|
||||
"""
|
||||
Return a shape that overlaps the rightmost entries of 'a' with the leftmost of 'b'.
|
||||
|
||||
|
|
@ -525,20 +525,24 @@ def shapeblender(a: _Tuple[int, ...],
|
|||
Shape of second array.
|
||||
keep_ones : bool, optional
|
||||
Treat innermost '1's as literal value instead of dimensional placeholder.
|
||||
Defaults to True.
|
||||
Defaults to False.
|
||||
|
||||
Examples
|
||||
--------
|
||||
>>> shapeblender((4,4,3),(3,2,1))
|
||||
(4,4,3,2,1)
|
||||
>>> shapeblender((1,2),(1,2,3))
|
||||
(1,2,3)
|
||||
>>> shapeblender((1,),(2,2,1))
|
||||
(1,2,2,1)
|
||||
>>> shapeblender((1,),(2,2,1),False)
|
||||
(2,2,1)
|
||||
>>> shapeblender((3,2),(3,2))
|
||||
(3,2)
|
||||
>>> shapeblender((4,3),(3,2))
|
||||
(4,3,2)
|
||||
>>> shapeblender((4,4),(3,2))
|
||||
(4,4,3,2)
|
||||
>>> shapeblender((1,2),(1,2,3))
|
||||
(1,2,3)
|
||||
>>> shapeblender((),(2,2,1))
|
||||
(2,2,1)
|
||||
>>> shapeblender((1,),(2,2,1))
|
||||
(2,2,1)
|
||||
>>> shapeblender((1,),(2,2,1),True)
|
||||
(1,2,2,1)
|
||||
|
||||
"""
|
||||
def is_broadcastable(a,b):
|
||||
|
|
|
|||
|
|
@ -358,7 +358,9 @@ class TestOrientation:
|
|||
a=a,b=b,c=c,
|
||||
alpha=alpha,beta=beta,gamma=gamma)
|
||||
assert o.to_pole(**{kw:vector,'with_symmetry':with_symmetry}).shape \
|
||||
== o.shape + vector.shape[:-1] + (o.symmetry_operations.shape if with_symmetry else ()) + vector.shape[-1:]
|
||||
== util.shapeblender(o.shape,vector.shape[:-1]) \
|
||||
+ (o.symmetry_operations.shape if with_symmetry else ()) \
|
||||
+ vector.shape[-1:]
|
||||
|
||||
@pytest.mark.parametrize('lattice',['hP','cI','cF']) #tI not included yet
|
||||
def test_Schmid(self,update,res_path,lattice):
|
||||
|
|
|
|||
|
|
@ -136,11 +136,13 @@ class TestUtil:
|
|||
((1,),(7,),False,(7,)),
|
||||
((1,),(7,),True,(1,7)),
|
||||
((2,),(2,2),False,(2,2)),
|
||||
((1,2),(2,2),False,(2,2)),
|
||||
((1,3),(2,3),False,(2,3)),
|
||||
((1,1,2),(2,2),False,(1,2,2)),
|
||||
((1,1,2),(2,2),True,(1,1,2,2)),
|
||||
((1,2,3),(2,3,4),False,(1,2,3,4)),
|
||||
((1,2,3),(1,2,3),False,(1,2,3)),
|
||||
((2,3,1,1),(2,3),False,(2,3,2,3)),
|
||||
((2,3,1,1),(2,3),True,(2,3,1,1,2,3)),
|
||||
])
|
||||
def test_shapeblender(self,a,b,ones,answer):
|
||||
assert util.shapeblender(a,b,ones) == answer
|
||||
|
|
|
|||
Loading…
Reference in New Issue