Merge branch '456-schmid-matrix-in-crystal' into 'development'

"Schmid" can be defined in Crystal

Closes #456

See merge request damask/DAMASK!1018

python/damask/_crystal.py

@@ -1302,3 +1302,51 @@ def relation_operations(self,
                         o.to_frame(hkl=o_p[:,1] if o_l != 'hP' else util.Bravais_to_Miller(hkil=o_p[:,1]))),
                         axis=-2)
         return (o_l,Rotation.from_parallel(a=m_p,b=o_p,active=True))
+
+
+    def Schmid(self, *,
+               N_slip: Optional[Union[IntSequence, Literal['*']]] = None,
+               N_twin: Optional[Union[IntSequence, Literal['*']]] = None) -> np.ndarray:
+        u"""
+        Calculate Schmid matrix P = d ⨂ n for selected deformation systems.
+
+        Parameters
+        ----------
+        N_slip|N_twin : '*' or sequence of int
+            Number of deformation systems per family of the deformation system.
+            Use '*' to select all.
+
+        Returns
+        -------
+        P : numpy.ndarray, shape (N,3,3)
+            Schmid matrix for each of the N deformation systems.
+
+        Examples
+        --------
+        Schmid matrix of first octahedral slip system of a face-centered
+        cubic crystal.
+
+        >>> import numpy as np
+        >>> import damask
+        >>> C = damask.Crystal(lattice='cF')
+        >>> C.Schmid(N_slip=[12])[0]
+        array([[ 0.    ,  0.    ,  0.    ],
+               [ 0.4082,  0.4082,  0.4082],
+               [-0.4082, -0.4082, -0.4082]])
+
+        """
+        if (N_slip is not None) ^ (N_twin is None):
+            raise KeyError('specify either "N_slip" or "N_twin"')
+
+        kinematics,active = (self.kinematics('slip'),N_slip) if N_twin is None else \
+                            (self.kinematics('twin'),N_twin)
+        everylen = list(map(len,kinematics['direction']))
+
+        if active == '*': active = everylen
+        if not active or (np.array(active) > everylen[:len(active)]).any():
+            raise ValueError('Invalid number of slip/twin systems')
+
+        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)]))
+        return np.einsum('...i,...j',d/np.linalg.norm(d,axis=-1,keepdims=True),
+                                     p/np.linalg.norm(p,axis=-1,keepdims=True))

python/damask/_orientation.py

@@ -987,24 +987,10 @@ def Schmid(self, *,
                [ 0.000,  0.000,  0.000]])
 
         """
-        if (N_slip is not None) ^ (N_twin is None):
-            raise KeyError('specify either "N_slip" or "N_twin"')
-
-        kinematics,active = (self.kinematics('slip'),N_slip) if N_twin is None else \
-                            (self.kinematics('twin'),N_twin)
-        if active == '*': active = [len(a) for a in kinematics['direction']]
-
-        if not active:
-            raise ValueError('Schmid matrix not defined')
-        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))
-
-        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)
+        if len(self.shape) == 0:
+            return self @ super().Schmid(N_slip=N_slip, N_twin=N_twin)
+        P = np.moveaxis(self @ super().Schmid(N_slip=N_slip, N_twin=N_twin),-3,0)
+        return P.reshape((P.shape[0],)+self.shape+(3,3))
 
 
     def related(self: MyType,

python/damask/_rotation.py

@@ -459,7 +459,7 @@ def __matmul__(self,
         >>> import damask
         >>> r = damask.Rotation.from_random(shape=(12))
         >>> o = np.ones((12,3,3))
-        >>> (r@o[np.newaxis,...]).shape                    # (12) @ (1,12, 3,3)
+        >>> (r@o[np.newaxis,...]).shape                                                             # (12) @ (1,12, 3,3)
         (12, 12, 3, 3)
 
         """

python/tests/test_Crystal.py

@@ -3,9 +3,16 @@
 from pathlib import Path
 
 import damask
+from damask import Table
 from damask import Crystal
 from damask import util
 
+
+@pytest.fixture
+def res_path(res_path_base):
+    """Directory containing testing resources."""
+    return res_path_base/'Crystal'
+
 class TestCrystal:
 
     @pytest.mark.parametrize('lattice,family',[('aP','cubic'),('xI','cubic')])
@@ -117,6 +124,21 @@ def test_N_twin(self,crystal,length):
         assert [len(s) for s in crystal.kinematics('twin')['direction']] == length
         assert [len(s) for s in crystal.kinematics('twin')['plane']] == length
 
+    @pytest.mark.parametrize('lattice',['hP','cI','cF','tI'])
+    def test_Schmid(self,update,res_path,lattice):
+        C = Crystal(lattice=lattice,c=(1.2 if lattice == 'tI' else None))                           # noqa
+        for mode in ['slip']+([] if lattice == 'tI' else ['twin']):
+            reference = res_path/f'{lattice}_{mode}.txt'
+            P = C.Schmid(N_slip='*') if mode == 'slip' else C.Schmid(N_twin='*')
+            if update:
+                Table({'Schmid':(3,3,)},P.reshape(-1,9)).save(reference)
+            assert np.allclose(P,Table.load(reference).get('Schmid'))
+
+    def test_Schmid_invalid(self):
+        with pytest.raises(KeyError):
+            Crystal(lattice='fcc').Schmid()
+
+    # https://doi.org/10.1016/0079-6425(94)00007-7, Fig. 22
     @pytest.mark.parametrize('c_a,mode',
                             [(np.sqrt(2)*0.99,['c','c','c','c']),
                              (np.sqrt(2)*1.01,['c','c','c','t']),

python/tests/test_Orientation.py

@@ -408,20 +408,6 @@ def test_to_frame_symmetries(self,np_rng,lattice,mode,vector,N_sym):
         assert shape_full//N_sym == shape_reduced
 
 
-    @pytest.mark.parametrize('lattice',['hP','cI','cF','tI'])
-    def test_Schmid(self,update,res_path,lattice):
-        O = Orientation(lattice=lattice,c=(1.2 if lattice == 'tI' else None))                       # noqa
-        for mode in ['slip']+([] if lattice == 'tI' else ['twin']):
-            reference = res_path/f'{lattice}_{mode}.txt'
-            P = O.Schmid(N_slip='*') if mode == 'slip' else O.Schmid(N_twin='*')
-            if update:
-                Table({'Schmid':(3,3,)},P.reshape(-1,9)).save(reference)
-            assert np.allclose(P,Table.load(reference).get('Schmid'))
-
-    def test_Schmid_invalid(self):
-        with pytest.raises(KeyError):
-            Orientation(lattice='fcc').Schmid()
-
     # https://doi.org/10.1016/0079-6425(94)00007-7, Fig. 22
     @pytest.mark.parametrize('c_a,mode',
                             [(np.sqrt(2)*0.99,['c','c','c','c']),