Merge branch 'smart-blend' into 'development'

robust Rotation/Orientation blending

See merge request damask/DAMASK!831

python/damask/_orientation.py

@@ -434,30 +434,20 @@ class Orientation(Rotation,Crystal):
         https://doi.org/10.1107/S0108767391006864
 
         """
-        rho = self.as_Rodrigues_vector(compact=True)*(1.0-1.0e-9)
+        def larger_or_equal(v,c):
+            return ((np.isclose(c[0],v[...,0]) | (v[...,0] > c[0])) &
+                    (np.isclose(c[1],v[...,1]) | (v[...,1] > c[1])) &
+                    (np.isclose(c[2],v[...,2]) | (v[...,2] > c[2]))).astype(bool)
 
-        with np.errstate(invalid='ignore'):
+        rho  = self.as_Rodrigues_vector(compact=True)
-            if   self.family == 'cubic':
+        return larger_or_equal(rho,
-                return ((rho[...,0] >= rho[...,1]) &
+                                     [rho[...,1],           rho[...,2],0] if self.family == 'cubic'
-                        (rho[...,1] >= rho[...,2]) &
+                                else [rho[...,1]*np.sqrt(3),0,         0] if self.family == 'hexagonal'
-                        (rho[...,2] >= 0)).astype(bool)
+                                else [rho[...,1],           0,         0] if self.family == 'tetragonal'
-            if self.family == 'hexagonal':
+                                else [0,                    0,         0] if self.family == 'orthorhombic'
-                return ((rho[...,0] >= rho[...,1]*np.sqrt(3)) &
+                                else [-np.inf,              0,         0] if self.family == 'monoclinic'
-                        (rho[...,1] >= 0) &
+                                else [-np.inf,        -np.inf,   -np.inf]) & self.in_FZ
-                        (rho[...,2] >= 0)).astype(bool)
-            if self.family == 'tetragonal':
-                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)
-            if self.family == 'monoclinic':
-                return ((rho[...,1] >= 0) &
-                        (rho[...,2] >= 0)).astype(bool)
 
-            return np.ones_like(rho[...,0],dtype=bool)
 
     def disorientation(self,
                        other: 'Orientation',
@@ -521,20 +511,21 @@ class Orientation(Rotation,Crystal):
         if self.family != other.family:
             raise NotImplementedError('disorientation between different crystal families')
 
-        blend = util.shapeblender(self.shape,other.shape)
+        blend = util.shapeblender( self.shape,other.shape)
-        s = self.equivalent
+        s_m   = util.shapeshifter( self.shape,blend,mode='right')
-        o = other.equivalent
+        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')
+        s =  self.broadcast_to(s_m).equivalent
-        o_ = o.reshape((1,o.shape[0])+other.shape).broadcast_to((s.shape[0],o.shape[0])+blend,mode='right')
+        o = other.broadcast_to(s_o).equivalent
-        r_ = s_.misorientation(o_)
+
+        r_ = s[:,np.newaxis,...].misorientation(o[np.newaxis,:,...]) # type: ignore[index]
         _r = ~r_
 
-        forward = r_.in_FZ & r_.in_disorientation_FZ
+        forward = r_.in_disorientation_FZ
-        reverse = _r.in_FZ & _r.in_disorientation_FZ
+        reverse = _r.in_disorientation_FZ
         ok  = forward | reverse
-        ok &= (np.cumsum(ok.reshape((-1,)+ok.shape[2:]),axis=0) == 1).reshape(ok.shape)
+        ok &= (np.cumsum(ok.reshape((-1,)+blend),axis=0) == 1).reshape(ok.shape)
-        r = np.where(np.any(forward[...,np.newaxis],axis=(0,1),keepdims=True),
+        r = np.where(np.any((ok&forward)[...,np.newaxis],axis=(0,1),keepdims=True),
                      r_.quaternion,
                      _r.quaternion)
         loc  = np.where(ok)
@@ -584,6 +575,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 +612,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])
+        blend = util.shapeblender( self.shape,vector_.shape[:-1])
-        poles = eq.broadcast_to(blend,mode='right') @ np.broadcast_to(vector_,blend+(3,))
+        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 +789,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) \
+            v = sym_ops.broadcast_to(s_v) @ v[...,np.newaxis,:]
-              @ np.broadcast_to(v.reshape(util.shapeshifter(v.shape,shape+(3,))),shape+(3,))
+
-        return ~(self.broadcast_to(blend))@np.broadcast_to(v,blend+(3,))
+        return ~(self.broadcast_to(blend)) @ np.broadcast_to(v,blend+(3,))
 
 
     def Schmid(self, *,
@@ -854,6 +849,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 +893,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'],

python/damask/_rotation.py

@@ -112,7 +112,7 @@ class Rotation:
 
 
     def __getitem__(self,
-                    item: Union[Tuple[int], int, bool, np.bool_, np.ndarray]):
+                    item: Union[Tuple[Union[None, int, slice]], int, bool, np.bool_, np.ndarray]):
         """
         Return self[item].
 
@@ -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]
+            blend = util.shapeblender( self.shape,other.shape)
-            p_m = self.quaternion[...,1:]
+            s_m   = util.shapeshifter( self.shape,blend,mode='right')
-            q_o = other.quaternion[...,0:1]
+            s_o   = util.shapeshifter(other.shape,blend,mode='left')
-            p_o = other.quaternion[...,1:]
+
+            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,7 +442,7 @@ 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,):
                     bs = util.shapeblender(self.shape,other.shape[:-l],False)

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):

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):

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