Merge branch 'shaped-rotation-application' into 'development'

Shaped rotation application

See merge request damask/DAMASK!819

python/damask/_orientation.py

@@ -804,7 +804,7 @@ class Orientation(Rotation,Crystal):
             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,))
+        return ~(self.broadcast_to(blend))@np.broadcast_to(v,blend+(3,))
 
 
     def Schmid(self, *,
@@ -833,7 +833,7 @@ class Orientation(Rotation,Crystal):
         >>> import damask
         >>> np.set_printoptions(3,suppress=True,floatmode='fixed')
         >>> O = damask.Orientation.from_Euler_angles(phi=[0,45,0],degrees=True,lattice='cF')
-        >>> O.Schmid(N_slip=[1])
+        >>> O.Schmid(N_slip=[12])[0]
         array([[ 0.000,  0.000,  0.000],
                [ 0.577, -0.000,  0.816],
                [ 0.000,  0.000,  0.000]])

python/damask/_rotation.py

@@ -375,6 +375,11 @@ class Rotation:
         Return self@other.
 
         Rotate vector, second-order tensor, or fourth-order tensor.
+        `other` is interpreted as an array of tensor quantities with the highest-possible order
+        considering the shape of `self`. Compatible innermost dimensions will blend.
+        For instance, shapes of (2,) and (3,3) for `self` and `other` prompt interpretation of
+        `other` as a second-rank tensor and result in (2,) rotated tensors, whereas
+        shapes of (2,1) and (3,3) for `self` and `other` result in (2,3) rotated vectors.
 
         Parameters
         ----------
@@ -386,29 +391,73 @@ class Rotation:
         rotated : numpy.ndarray, shape (...,3), (...,3,3), or (...,3,3,3,3)
             Rotated vector or tensor, i.e. transformed to frame defined by rotation.
 
+        Examples
+        --------
+        All below examples rely on imported modules:
+        >>> import numpy as np
+        >>> import damask
+
+        Application of twelve (random) rotations to a set of five vectors.
+
+        >>> r = damask.Rotation.from_random(shape=(12))
+        >>> o = np.ones((5,3))
+        >>> (r@o).shape                                    # (12) @ (5, 3)
+        (12,5, 3)
+
+        Application of a (random) rotation to all twelve second-rank tensors.
+
+        >>> r = damask.Rotation.from_random()
+        >>> o = np.ones((12,3,3))
+        >>> (r@o).shape                                    # (1) @ (12, 3,3)
+        (12,3,3)
+
+        Application of twelve (random) rotations to the corresponding twelve second-rank tensors.
+
+        >>> r = damask.Rotation.from_random(shape=(12))
+        >>> o = np.ones((12,3,3))
+        >>> (r@o).shape                                    # (12) @ (3,3)
+        (12,3,3)
+
+        Application of each of three (random) rotations to all three vectors.
+
+        >>> r = damask.Rotation.from_random(shape=(3))
+        >>> o = np.ones((3,3))
+        >>> (r[...,np.newaxis]@o[np.newaxis,...]).shape    # (3,1) @ (1,3, 3)
+        (3,3,3)
+
+        Application of twelve (random) rotations to all twelve second-rank tensors.
+
+        >>> r = damask.Rotation.from_random(shape=(12))
+        >>> o = np.ones((12,3,3))
+        >>> (r@o[np.newaxis,...]).shape                    # (12) @ (1,12, 3,3)
+        (12,3,3,3)
+
         """
         if isinstance(other, np.ndarray):
-            if self.shape + (3,) == other.shape:
+            obs = util.shapeblender(self.shape,other.shape,keep_ones=False)[len(self.shape):]
-                q_m = self.quaternion[...,0]
+            for l in [4,2,1]:
-                p_m = self.quaternion[...,1:]
+                if obs[-l:] == l*(3,):
-                A = q_m**2 - np.einsum('...i,...i',p_m,p_m)
+                    bs = util.shapeblender(self.shape,other.shape[:-l],False)
-                B = 2. * np.einsum('...i,...i',p_m,other)
+                    self_ = self.broadcast_to(bs) if self.shape != bs else self
-                C = 2. * _P * q_m
+                    if l==1:
-                return np.block([(A * other[...,i]).reshape(self.shape+(1,)) +
+                        q_m = self_.quaternion[...,0]
-                                 (B * p_m[...,i]).reshape(self.shape+(1,)) +
+                        p_m = self_.quaternion[...,1:]
-                                 (C * (  p_m[...,(i+1)%3]*other[...,(i+2)%3]\
+                        A = q_m**2 - np.einsum('...i,...i',p_m,p_m)
-                                       - p_m[...,(i+2)%3]*other[...,(i+1)%3])).reshape(self.shape+(1,))
+                        B = 2. * np.einsum('...i,...i',p_m,other)
-                                 for i in [0,1,2]])
+                        C = 2. * _P * q_m
-            if self.shape + (3,3) == other.shape:
+                        return np.block([(A * other[...,i]) +
-                R = self.as_matrix()
+                                         (B *   p_m[...,i]) +
-                return np.einsum('...im,...jn,...mn',R,R,other)
+                                         (C * ( p_m[...,(i+1)%3]*other[...,(i+2)%3]
-            if self.shape + (3,3,3,3) == other.shape:
+                                              - p_m[...,(i+2)%3]*other[...,(i+1)%3]))
-                R = self.as_matrix()
+                                        for i in [0,1,2]]).reshape(bs+(3,),order='F')
-                return np.einsum('...im,...jn,...ko,...lp,...mnop',R,R,R,R,other)
+                    else:
-            else:
+                        return np.einsum({2: '...im,...jn,...mn',
-                raise ValueError('can only rotate vectors, second-order tensors, and fourth-order tensors')
+                                          4: '...im,...jn,...ko,...lp,...mnop'}[l],
+                                         *l*[self_.as_matrix()],
+                                         other)
+            raise ValueError('can only rotate vectors, second-order tensors, and fourth-order tensors')
         elif isinstance(other, Rotation):
-            raise TypeError('use "R1*R2", i.e. multiplication, to compose rotations "R1" and "R2"')
+            raise TypeError('use "R2*R1", i.e. multiplication, to compose rotations "R1" and "R2"')
         else:
             raise TypeError(f'cannot rotate "{type(other)}"')
 

python/damask/util.py

@@ -512,7 +512,8 @@ def shapeshifter(fro: _Tuple[int, ...],
     return tuple(final_shape[::-1] if mode == 'left' else final_shape)
 
 def shapeblender(a: _Tuple[int, ...],
-                 b: _Tuple[int, ...]) -> _Tuple[int, ...]:
+                 b: _Tuple[int, ...],
+                 keep_ones: bool = True) -> _Tuple[int, ...]:
     """
     Return a shape that overlaps the rightmost entries of 'a' with the leftmost of 'b'.
 
@@ -522,6 +523,9 @@ def shapeblender(a: _Tuple[int, ...],
         Shape of first array.
     b : tuple
         Shape of second array.
+    keep_ones : bool, optional
+        Treat innermost '1's as literal value instead of dimensional placeholder.
+        Defaults to True.
 
     Examples
     --------
@@ -531,13 +535,30 @@ def shapeblender(a: _Tuple[int, ...],
         (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)
 
     """
-    i = min(len(a),len(b))
+    def is_broadcastable(a,b):
-    while i > 0 and a[-i:] != b[:i]: i -= 1
+        try:
-    return a + b[i:]
+            _np.broadcast_shapes(a,b)
+            return True
+        except ValueError:
+            return False
+
+    a_,_b = a,b
+    if keep_ones:
+        i = min(len(a_),len(_b))
+        while i > 0 and a_[-i:] != _b[:i]: i -= 1
+        return a_ + _b[i:]
+    else:
+        a_ += max(0,len(_b)-len(a_))*(1,)
+        while not is_broadcastable(a_,_b):
+            a_ = a_ + ((1,) if len(a_)<=len(_b) else ())
+            _b = ((1,) if len(_b)<len(a_) else ()) + _b
+        return _np.broadcast_shapes(a_,_b)
 
 
 def _docstringer(docstring: _Union[str, _Callable],

python/tests/test_Orientation.py

@@ -162,7 +162,7 @@ class TestOrientation:
                                                                ([np.arccos(3**(-.5)),np.pi/4,0],[0,0],[0,0,1],[0,0,1])])
     def test_fiber_IPF(self,crystal,sample,direction,color):
         fiber = Orientation.from_fiber_component(crystal=crystal,sample=sample,family='cubic',shape=200)
-        print(np.allclose(fiber.IPF_color(direction),color))
+        assert np.allclose(fiber.IPF_color(direction),color)
 
 
     @pytest.mark.parametrize('kwargs',[
@@ -455,11 +455,9 @@ class TestOrientation:
         p = Orientation.from_random(family=family,shape=right)
         blend = util.shapeblender(o.shape,p.shape)
         for loc in np.random.randint(0,blend,(10,len(blend))):
-            # print(f'{a}/{b} @ {loc}')
+            l = () if  left is None else tuple(np.minimum(np.array(left )-1,loc[:len(left)]))
-            # print(o[tuple(loc[:len(o.shape)])].disorientation(p[tuple(loc[-len(p.shape):])]))
+            r = () if right is None else tuple(np.minimum(np.array(right)-1,loc[-len(right):]))
-            # print(o.disorientation(p)[tuple(loc)])
+            assert o[l].disorientation(p[r]).isclose(o.disorientation(p)[tuple(loc)])
-            assert o[tuple(loc[:len(o.shape)])].disorientation(p[tuple(loc[-len(p.shape):])]) \
-                   .isclose(o.disorientation(p)[tuple(loc)])
 
     @pytest.mark.parametrize('family',crystal_families)
     @pytest.mark.parametrize('left,right',[
@@ -467,13 +465,16 @@ class TestOrientation:
                                     ((2,2),(4,4)),
                                     ((3,1),(1,3)),
                                     (None,(3,)),
+                                    (None,()),
                                    ])
     def test_IPF_color_blending(self,family,left,right):
         o = Orientation.from_random(family=family,shape=left)
         v = np.random.random(right+(3,))
         blend = util.shapeblender(o.shape,v.shape[:-1])
         for loc in np.random.randint(0,blend,(10,len(blend))):
-            assert np.allclose(o[tuple(loc[:len(o.shape)])].IPF_color(v[tuple(loc[-len(v.shape[:-1]):])]),
+            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].IPF_color(v[r]),
                                o.IPF_color(v)[tuple(loc)])
 
     @pytest.mark.parametrize('family',crystal_families)
@@ -488,7 +489,9 @@ class TestOrientation:
         v = np.random.random(right+(3,))
         blend = util.shapeblender(o.shape,v.shape[:-1])
         for loc in np.random.randint(0,blend,(10,len(blend))):
-            assert np.allclose(o[tuple(loc[:len(o.shape)])].to_SST(v[tuple(loc[-len(v.shape[:-1]):])]),
+            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_SST(v[r]),
                                o.to_SST(v)[tuple(loc)])
 
     @pytest.mark.parametrize('lattice,a,b,c,alpha,beta,gamma',
@@ -514,8 +517,10 @@ class TestOrientation:
         v = np.random.random(right+(3,))
         blend = util.shapeblender(o.shape,v.shape[:-1])
         for loc in np.random.randint(0,blend,(10,len(blend))):
-            assert np.allclose(o[tuple(loc[:len(o.shape)])].to_pole(uvw=v[tuple(loc[-len(v.shape[:-1]):])]),
+            l = () if  left is None else tuple(np.minimum(np.array(left )-1,loc[:len(left)]))
-                               o.to_pole(uvw=v)[tuple(loc)])
+            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)])
 
     def test_mul_invalid(self):
         with pytest.raises(TypeError):

python/tests/test_Rotation.py

@@ -1065,7 +1065,7 @@ class TestRotation:
 
     @pytest.mark.parametrize('data',[np.random.rand(4),
                                      np.random.rand(3,2),
-                                     np.random.rand(3,2,3,3)])
+                                     np.random.rand(3,3,3,1)])
     def test_rotate_invalid_shape(self,data):
         R = Rotation.from_random()
         with pytest.raises(ValueError):

python/tests/test_util.py

@@ -128,18 +128,22 @@ class TestUtil:
         with pytest.raises(ValueError):
             util.shapeshifter(fro,to,mode)
 
-    @pytest.mark.parametrize('a,b,answer',
+    @pytest.mark.parametrize('a,b,ones,answer',
                              [
-                              ((),(1,),(1,)),
+                              ((),(1,),True,(1,)),
-                              ((1,),(),(1,)),
+                              ((1,),(),False,(1,)),
-                              ((1,),(7,),(1,7)),
+                              ((1,1),(7,),False,(1,7)),
-                              ((2,),(2,2),(2,2)),
+                              ((1,),(7,),False,(7,)),
-                              ((1,2),(2,2),(1,2,2)),
+                              ((1,),(7,),True,(1,7)),
-                              ((1,2,3),(2,3,4),(1,2,3,4)),
+                              ((2,),(2,2),False,(2,2)),
-                              ((1,2,3),(1,2,3),(1,2,3)),
+                              ((1,2),(2,2),False,(2,2)),
+                              ((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)),
                              ])
-    def test_shapeblender(self,a,b,answer):
+    def test_shapeblender(self,a,b,ones,answer):
-        assert util.shapeblender(a,b) == answer
+        assert util.shapeblender(a,b,ones) == answer
 
     @pytest.mark.parametrize('style',[util.emph,util.deemph,util.warn,util.strikeout])
     def test_decorate(self,style):