From 550c757cdc87d913f7be9747a14d0efab1321e85 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Wed, 4 Oct 2023 09:34:23 -0400 Subject: [PATCH] overall smart blending; more robust disFZ --- python/damask/_orientation.py | 84 +++++++++++++++++--------------- python/damask/_rotation.py | 19 ++++++-- python/damask/util.py | 24 +++++---- python/tests/test_Orientation.py | 4 +- python/tests/test_util.py | 4 +- 5 files changed, 80 insertions(+), 55 deletions(-) diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index 15f8464fd..ef71366c8 100644 --- a/python/damask/_orientation.py +++ b/python/damask/_orientation.py @@ -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'], diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index c03bff4b2..21197f105 100644 --- a/python/damask/_rotation.py +++ b/python/damask/_rotation.py @@ -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: diff --git a/python/damask/util.py b/python/damask/util.py index 242c259b3..513cab87f 100644 --- a/python/damask/util.py +++ b/python/damask/util.py @@ -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): diff --git a/python/tests/test_Orientation.py b/python/tests/test_Orientation.py index 9836ee16e..6de2c99b1 100644 --- a/python/tests/test_Orientation.py +++ b/python/tests/test_Orientation.py @@ -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): diff --git a/python/tests/test_util.py b/python/tests/test_util.py index d5ec36d75..ce5d1b3a9 100644 --- a/python/tests/test_util.py +++ b/python/tests/test_util.py @@ -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