Merge branch 'misc-improvements' into 'development'

Misc improvements

See merge request damask/DAMASK!181

examples/ConfigFiles/Phase_DisloUCLA_Tungsten.config

@@ -27,7 +27,7 @@ SolidSolutionStrength 0.0           # Strength due to elements in solid solution
 
 ### Dislocation glide parameters ###
 #per family
-Nslip               12  0        
+Nslip               12
 slipburgers         2.72e-10                # Burgers vector of slip system [m]
 rhoedge0            1.0e12                  # Initial edge dislocation density [m/m**3]
 rhoedgedip0         1.0                     # Initial edged dipole dislocation density [m/m**3]

python/damask/_rotation.py

@@ -20,8 +20,8 @@ class Rotation:
       when viewing from the end point of the rotation axis towards the origin.
     - rotations will be interpreted in the passive sense.
     - Euler angle triplets are implemented using the Bunge convention,
-      with the angular ranges as [0, 2π],[0, π],[0, 2π].
+      with the angular ranges as [0,2π], [0,π], [0,2π].
-    - the rotation angle ω is limited to the interval [0, π].
+    - the rotation angle ω is limited to the interval [0,π].
     - the real part of a quaternion is positive, Re(q) > 0
     - P = -1 (as default).
 
@@ -49,7 +49,8 @@ class Rotation:
         Parameters
         ----------
         quaternion : numpy.ndarray, optional
-            Unit quaternion that follows the conventions. Use .from_quaternion to perform a sanity check.
+            Unit quaternion in positive real hemisphere.
+            Use .from_quaternion to perform a sanity check.
 
         """
         self.quaternion = quaternion.copy()
@@ -73,7 +74,7 @@ class Rotation:
             raise NotImplementedError('Support for multiple rotations missing')
         return '\n'.join([
                'Quaternion: (real={:.3f}, imag=<{:+.3f}, {:+.3f}, {:+.3f}>)'.format(*(self.quaternion)),
-               'Matrix:\n{}'.format(self.as_matrix()),
+               'Matrix:\n{}'.format(np.round(self.as_matrix(),8)),
                'Bunge Eulers / deg: ({:3.2f}, {:3.2f}, {:3.2f})'.format(*self.as_Eulers(degrees=True)),
                 ])
 
@@ -87,10 +88,6 @@ class Rotation:
         other : numpy.ndarray or Rotation
             Vector, second or fourth order tensor, or rotation object that is rotated.
 
-        Todo
-        ----
-        Check rotation of 4th order tensor
-
         """
         if isinstance(other, Rotation):
             q_m = self.quaternion[...,0:1]
@@ -99,7 +96,7 @@ class Rotation:
             p_o = other.quaternion[...,1:]
             q = (q_m*q_o - np.einsum('...i,...i',p_m,p_o).reshape(self.shape+(1,)))
             p = q_m*p_o + q_o*p_m + _P * np.cross(p_m,p_o)
-            return self.__class__(np.block([q,p])).standardize()
+            return self.__class__(np.block([q,p]))._standardize()
 
         elif isinstance(other,np.ndarray):
             if self.shape + (3,) == other.shape:
@@ -124,24 +121,24 @@ class Rotation:
         else:
             raise TypeError('Cannot rotate {}'.format(type(other)))
 
-    def inverse(self):
-        """In-place inverse rotation/backward rotation."""
-        self.quaternion[...,1:] *= -1
-        return self
 
-    def inversed(self):
+    def _standardize(self):
-        """Inverse rotation/backward rotation."""
+        """Standardize (ensure positive real hemisphere)."""
-        return self.copy().inverse()
-
-
-    def standardize(self):
-        """In-place quaternion representation with positive real part."""
         self.quaternion[self.quaternion[...,0] < 0.0] *= -1
         return self
 
-    def standardized(self):
+    def inverse(self):
-        """Quaternion representation with positive real part."""
+        """In-place inverse rotation (backward rotation)."""
-        return self.copy().standardize()
+        self.quaternion[...,1:] *= -1
+        return self
+
+    def __invert__(self):
+        """Inverse rotation (backward rotation)."""
+        return self.copy().inverse()
+
+    def inversed(self):
+        """Inverse rotation (backward rotation)."""
+        return ~ self
 
 
     def misorientation(self,other):
@@ -154,7 +151,7 @@ class Rotation:
             Rotation to which the misorientation is computed.
 
         """
-        return other*self.inversed()
+        return other@~self
 
 
     def broadcast_to(self,shape):
@@ -169,7 +166,7 @@ class Rotation:
         return self.__class__(q)
 
 
-    def average(self,other):
+    def average(self,other): #ToDo: discuss calling for vectors
         """
         Calculate the average rotation.
 
@@ -189,25 +186,31 @@ class Rotation:
 
     def as_quaternion(self):
         """
-        Unit quaternion [q, p_1, p_2, p_3].
+        Represent as unit quaternion.
 
-        Parameters
+        Returns
-        ----------
+        -------
-        quaternion : bool, optional
+        q : numpy.ndarray of shape (...,4)
-            return quaternion as DAMASK object.
+            Unit quaternion in positive real hemisphere: (q_0, q_1, q_2, q_3), |q|=1, q_0 ≥ 0.
 
         """
-        return self.quaternion
+        return self.quaternion.copy()
 
     def as_Eulers(self,
                   degrees = False):
         """
-        Bunge-Euler angles: (φ_1, ϕ, φ_2).
+        Represent as Bunge-Euler angles.
 
         Parameters
         ----------
         degrees : bool, optional
-            return angles in degrees.
+            Return angles in degrees.
+
+        Returns
+        -------
+        phi : numpy.ndarray of shape (...,3)
+            Bunge-Euler angles: (φ_1, ϕ, φ_2), φ_1 ∈ [0,2π], ϕ ∈ [0,π], φ_2 ∈ [0,2π]
+            unless degrees == True: φ_1 ∈ [0,360], ϕ ∈ [0,180], φ_2 ∈ [0,360]
 
         """
         eu = Rotation._qu2eu(self.quaternion)
@@ -218,14 +221,21 @@ class Rotation:
                       degrees = False,
                       pair = False):
         """
-        Axis angle representation [n_1, n_2, n_3, ω] unless pair == True: ([n_1, n_2, n_3], ω).
+        Represent as axis angle pair.
 
         Parameters
         ----------
         degrees : bool, optional
-            return rotation angle in degrees.
+            Return rotation angle in degrees. Defaults to False.
         pair : bool, optional
-            return tuple of axis and angle.
+            Return tuple of axis and angle. Defaults to False.
+
+        Returns
+        -------
+        axis_angle : numpy.ndarray of shape (...,4) unless pair == True:
+            tuple containing numpy.ndarray of shapes (...,3) and (...)
+            Axis angle pair: (n_1, n_2, n_3, ω), |n| = 1 and ω ∈ [0,π]
+            unless degrees = True: ω ∈ [0,180].
 
         """
         ax = Rotation._qu2ax(self.quaternion)
@@ -233,29 +243,60 @@ class Rotation:
         return (ax[...,:3],ax[...,3]) if pair else ax
 
     def as_matrix(self):
-        """Rotation matrix."""
+        """
+        Represent as rotation matrix.
+
+        Returns
+        -------
+        R : numpy.ndarray of shape (...,3,3)
+            Rotation matrix R, det(R) = 1, R.T∙R=I.
+
+        """
         return Rotation._qu2om(self.quaternion)
 
     def as_Rodrigues(self,
                      vector = False):
         """
-        Rodrigues-Frank vector representation [n_1, n_2, n_3, tan(ω/2)] unless vector == True: [n_1, n_2, n_3] * tan(ω/2).
+        Represent as Rodrigues-Frank vector with separated axis and angle argument.
 
         Parameters
         ----------
         vector : bool, optional
-            return as actual Rodrigues--Frank vector, i.e. rotation axis scaled by tan(ω/2).
+            Return as actual Rodrigues-Frank vector, i.e. axis
+            and angle argument are not separated.
+
+        Returns
+        -------
+        rho : numpy.ndarray of shape (...,4) unless vector == True:
+            numpy.ndarray of shape (...,3)
+            Rodrigues-Frank vector: [n_1, n_2, n_3, tan(ω/2)], |n| = 1 and ω ∈ [0,π].
 
         """
         ro = Rotation._qu2ro(self.quaternion)
         return ro[...,:3]*ro[...,3] if vector else ro
 
     def as_homochoric(self):
-        """Homochoric vector: (h_1, h_2, h_3)."""
+        """
+        Represent as homochoric vector.
+
+        Returns
+        -------
+        h : numpy.ndarray of shape (...,3)
+            Homochoric vector: (h_1, h_2, h_3), |h| < 1/2*π^(2/3).
+
+        """
         return Rotation._qu2ho(self.quaternion)
 
     def as_cubochoric(self):
-        """Cubochoric vector: (c_1, c_2, c_3)."""
+        """
+        Represent as cubochoric vector.
+
+        Returns
+        -------
+        c : numpy.ndarray of shape (...,3)
+              Cubochoric vector: (c_1, c_2, c_3), max(c_i) < 1/2*π^(2/3).
+
+        """
         return Rotation._qu2cu(self.quaternion)
 
     def M(self): # ToDo not sure about the name: as_M or M? we do not have a from_M
@@ -275,18 +316,34 @@ class Rotation:
     # Static constructors. The input data needs to follow the conventions, options allow to
     # relax the conventions.
     @staticmethod
-    def from_quaternion(quaternion,
+    def from_quaternion(q,
                         accept_homomorph = False,
                         P = -1,
-                        acceptHomomorph = None):
+                        acceptHomomorph = None): # old name (for compatibility)
+        """
+        Initialize from quaternion.
 
+        Parameters
+        ----------
+        q : numpy.ndarray of shape (...,4)
+            Unit quaternion in positive real hemisphere: (q_0, q_1, q_2, q_3),
+            |q|=1, q_0 ≥ 0.
+        accept_homomorph : boolean, optional
+            Allow homomorphic variants, i.e. q_0 < 0 (negative real hemisphere).
+            Defaults to False.
+        P : integer ∈ {-1,1}, optional
+            Convention used. Defaults to -1.
+
+        """
         if acceptHomomorph is not None:
-            accept_homomorph = acceptHomomorph
+            accept_homomorph = acceptHomomorph # for compatibility
-        qu = np.array(quaternion,dtype=float)
+        qu = np.array(q,dtype=float)
         if qu.shape[:-2:-1] != (4,):
             raise ValueError('Invalid shape.')
+        if abs(P) != 1:
+            raise ValueError('P ∉ {-1,1}')
 
-        if P > 0: qu[...,1:4] *= -1                                                                 # convert from P=1 to P=-1
+        if P == 1: qu[...,1:4] *= -1
         if accept_homomorph:
             qu[qu[...,0] < 0.0] *= -1
         else:
@@ -298,10 +355,21 @@ class Rotation:
         return Rotation(qu)
 
     @staticmethod
-    def from_Eulers(eulers,
+    def from_Eulers(phi,
                     degrees = False):
+        """
+        Initialize from Bunge-Euler angles.
 
-        eu = np.array(eulers,dtype=float)
+        Parameters
+        ----------
+        phi : numpy.ndarray of shape (...,3)
+            Bunge-Euler angles: (φ_1, ϕ, φ_2), φ_1 ∈ [0,2π], ϕ ∈ [0,π], φ_2 ∈ [0,2π]
+            unless degrees == True: φ_1 ∈ [0,360], ϕ ∈ [0,180], φ_2 ∈ [0,360].
+        degrees : boolean, optional
+            Bunge-Euler angles are given in degrees. Defaults to False.
+
+        """
+        eu = np.array(phi,dtype=float)
         if eu.shape[:-2:-1] != (3,):
             raise ValueError('Invalid shape.')
 
@@ -316,12 +384,29 @@ class Rotation:
                         degrees = False,
                         normalise = False,
                         P = -1):
+        """
+        Initialize from Axis angle pair.
 
+        Parameters
+        ----------
+        axis_angle : numpy.ndarray of shape (...,4)
+            Axis angle pair: [n_1, n_2, n_3, ω], |n| = 1 and ω ∈ [0,π]
+            unless degrees = True: ω ∈ [0,180].
+        degrees : boolean, optional
+            Angle ω is given in degrees. Defaults to False.
+        normalize: boolean, optional
+            Allow |n| ≠ 1. Defaults to False.
+        P : integer ∈ {-1,1}, optional
+            Convention used. Defaults to -1.
+
+        """
         ax = np.array(axis_angle,dtype=float)
         if ax.shape[:-2:-1] != (4,):
             raise ValueError('Invalid shape.')
+        if abs(P) != 1:
+            raise ValueError('P ∉ {-1,1}')
 
-        if P > 0:     ax[...,0:3] *= -1                                                             # convert from P=1 to P=-1
+        if P == 1:    ax[...,0:3] *= -1
         if degrees:   ax[...,  3]  = np.radians(ax[...,3])
         if normalise: ax[...,0:3] /= np.linalg.norm(ax[...,0:3],axis=-1)
         if np.any(ax[...,3] < 0.0) or np.any(ax[...,3] > np.pi):
@@ -335,7 +420,19 @@ class Rotation:
     def from_basis(basis,
                    orthonormal = True,
                    reciprocal = False):
+        """
+        Initialize from lattice basis vectors.
 
+        Parameters
+        ----------
+        basis : numpy.ndarray of shape (...,3,3)
+            Three lattice basis vectors in three dimensions.
+        orthonormal : boolean, optional
+            Basis is strictly orthonormal, i.e. is free of stretch components. Defaults to True.
+        reciprocal : boolean, optional
+            Basis vectors are given in reciprocal (instead of real) space. Defaults to False.
+
+        """
         om = np.array(basis,dtype=float)
         if om.shape[:-3:-1] != (3,3):
             raise ValueError('Invalid shape.')
@@ -356,20 +453,43 @@ class Rotation:
         return Rotation(Rotation._om2qu(om))
 
     @staticmethod
-    def from_matrix(om):
+    def from_matrix(R):
+        """
+        Initialize from rotation matrix.
 
-        return Rotation.from_basis(om)
+        Parameters
+        ----------
+        R : numpy.ndarray of shape (...,3,3)
+            Rotation matrix: det(R) = 1, R.T∙R=I.
+
+        """
+        return Rotation.from_basis(R)
 
     @staticmethod
-    def from_Rodrigues(rodrigues,
+    def from_Rodrigues(rho,
                        normalise = False,
                        P = -1):
+        """
+        Initialize from Rodrigues-Frank vector.
 
-        ro = np.array(rodrigues,dtype=float)
+        Parameters
+        ----------
+        rho : numpy.ndarray of shape (...,4)
+            Rodrigues-Frank vector (angle separated from axis).
+            (n_1, n_2, n_3, tan(ω/2)), |n| = 1  and ω ∈ [0,π].
+        normalize : boolean, optional
+            Allow |n| ≠ 1. Defaults to False.
+        P : integer ∈ {-1,1}, optional
+            Convention used. Defaults to -1.
+
+        """
+        ro = np.array(rho,dtype=float)
         if ro.shape[:-2:-1] != (4,):
             raise ValueError('Invalid shape.')
+        if abs(P) != 1:
+            raise ValueError('P ∉ {-1,1}')
 
-        if P > 0:     ro[...,0:3] *= -1                                                             # convert from P=1 to P=-1
+        if P == 1:    ro[...,0:3] *= -1
         if normalise: ro[...,0:3] /= np.linalg.norm(ro[...,0:3],axis=-1)
         if np.any(ro[...,3] < 0.0):
             raise ValueError('Rodrigues vector rotation angle not positive.')
@@ -379,14 +499,26 @@ class Rotation:
         return Rotation(Rotation._ro2qu(ro))
 
     @staticmethod
-    def from_homochoric(homochoric,
+    def from_homochoric(h,
                         P = -1):
+        """
+        Initialize from homochoric vector.
 
-        ho = np.array(homochoric,dtype=float)
+        Parameters
+        ----------
+        h : numpy.ndarray of shape (...,3)
+            Homochoric vector: (h_1, h_2, h_3), |h| < (3/4*π)^(1/3).
+        P : integer ∈ {-1,1}, optional
+            Convention used. Defaults to -1.
+
+        """
+        ho = np.array(h,dtype=float)
         if ho.shape[:-2:-1] != (3,):
             raise ValueError('Invalid shape.')
+        if abs(P) != 1:
+            raise ValueError('P ∉ {-1,1}')
 
-        if P > 0: ho *= -1                                                                          # convert from P=1 to P=-1
+        if P == 1: ho *= -1
 
         if np.any(np.linalg.norm(ho,axis=-1) >_R1+1e-9):
             raise ValueError('Homochoric coordinate outside of the sphere.')
@@ -394,18 +526,30 @@ class Rotation:
         return Rotation(Rotation._ho2qu(ho))
 
     @staticmethod
-    def from_cubochoric(cubochoric,
+    def from_cubochoric(c,
-                       P = -1):
+                        P = -1):
+        """
+        Initialize from cubochoric vector.
 
-        cu = np.array(cubochoric,dtype=float)
+        Parameters
+        ----------
+        c : numpy.ndarray of shape (...,3)
+            Cubochoric vector: (c_1, c_2, c_3), max(c_i) < 1/2*π^(2/3).
+        P : integer ∈ {-1,1}, optional
+            Convention used. Defaults to -1.
+
+        """
+        cu = np.array(c,dtype=float)
         if cu.shape[:-2:-1] != (3,):
             raise ValueError('Invalid shape.')
+        if abs(P) != 1:
+            raise ValueError('P ∉ {-1,1}')
 
-        if np.abs(np.max(cu))>np.pi**(2./3.) * 0.5+1e-9:
+        if np.abs(np.max(cu)) > np.pi**(2./3.) * 0.5+1e-9:
-            raise ValueError('Cubochoric coordinate outside of the cube: {} {} {}.'.format(*cu))
+            raise ValueError('Cubochoric coordinate outside of the cube.')
 
         ho = Rotation._cu2ho(cu)
-        if P > 0: ho *= -1                                                                          # convert from P=1 to P=-1
+        if P == 1: ho *= -1
 
         return Rotation(Rotation._ho2qu(ho))
 
@@ -458,7 +602,7 @@ class Rotation:
                       np.cos(2.0*np.pi*r[...,1])*B,
                       np.sin(2.0*np.pi*r[...,0])*A],axis=-1)
 
-        return Rotation(q.reshape(r.shape[:-1]+(4,)) if shape is not None else q).standardize()
+        return Rotation(q.reshape(r.shape[:-1]+(4,)) if shape is not None else q)._standardize()
 
 
     # for compatibility (old names do not follow convention)
@@ -471,7 +615,7 @@ class Rotation:
 ####################################################################################################
 # Code below available according to the following conditions on https://github.com/MarDiehl/3Drotations
 ####################################################################################################
-# Copyright (c) 2017-2019, Martin Diehl/Max-Planck-Institut für Eisenforschung GmbH
+# Copyright (c) 2017-2020, Martin Diehl/Max-Planck-Institut für Eisenforschung GmbH
 # Copyright (c) 2013-2014, Marc De Graef/Carnegie Mellon University
 # All rights reserved.
 #
@@ -530,7 +674,7 @@ class Rotation:
                           np.zeros(qu.shape[:-1]+(2,))]),
                       np.where(np.abs(q03_s) < 1.0e-8,
                           np.block([np.arctan2(   2.0*qu[...,1:2]*qu[...,2:3],qu[...,1:2]**2-qu[...,2:3]**2),
-                                    np.broadcast_to(np.pi,qu.shape[:-1]+(1,)),
+                                    np.broadcast_to(np.pi,qu[...,0:1].shape),
                                     np.zeros(qu.shape[:-1]+(1,))]),
                           np.block([np.arctan2((-_P*q02+q13)*chi, (-_P*q01-q23)*chi),
                                     np.arctan2( 2.0*chi,          q03_s-q12_s    ),
@@ -553,7 +697,7 @@ class Rotation:
             s = np.sign(qu[...,0:1])/np.sqrt(qu[...,1:2]**2+qu[...,2:3]**2+qu[...,3:4]**2)
             omega = 2.0 * np.arccos(np.clip(qu[...,0:1],-1.0,1.0))
             ax = np.where(np.broadcast_to(qu[...,0:1] < 1.0e-8,qu.shape),
-                          np.block([qu[...,1:4],np.broadcast_to(np.pi,qu.shape[:-1]+(1,))]),
+                          np.block([qu[...,1:4],np.broadcast_to(np.pi,qu[...,0:1].shape)]),
                           np.block([qu[...,1:4]*s,omega]))
         ax[np.isclose(qu[...,0],1.,rtol=0.0)] = [0.0, 0.0, 1.0, 0.0]
         return ax
@@ -564,7 +708,7 @@ class Rotation:
         with np.errstate(invalid='ignore',divide='ignore'):
             s  = np.linalg.norm(qu[...,1:4],axis=-1,keepdims=True)
             ro = np.where(np.broadcast_to(np.abs(qu[...,0:1]) < 1.0e-12,qu.shape),
-                          np.block([qu[...,1:2], qu[...,2:3], qu[...,3:4], np.broadcast_to(np.inf,qu.shape[:-1]+(1,))]),
+                          np.block([qu[...,1:2], qu[...,2:3], qu[...,3:4], np.broadcast_to(np.inf,qu[...,0:1].shape)]),
                           np.block([qu[...,1:2]/s,qu[...,2:3]/s,qu[...,3:4]/s,
                                     np.tan(np.arccos(np.clip(qu[...,0:1],-1.0,1.0)))
                                    ])

python/damask/util.py

@@ -93,7 +93,7 @@ def strikeout(what):
 
 
 def execute(cmd,
-            streamIn = None,
+            stream_in = None,
             wd = './',
             env = None):
     """
@@ -103,7 +103,7 @@ def execute(cmd,
     ----------
     cmd : str
         Command to be executed.
-    streanIn :, optional
+    stream_in : file object, optional
         Input (via pipe) for executed process.
     wd : str, optional
         Working directory of process. Defaults to ./ .
@@ -119,14 +119,14 @@ def execute(cmd,
                                stderr = subprocess.PIPE,
                                stdin  = subprocess.PIPE,
                                env = myEnv)
-    out,error = [i for i in (process.communicate() if streamIn is None
+    stdout, stderr = [i for i in (process.communicate() if stream_in is None
-                                                   else process.communicate(streamIn.read().encode('utf-8')))]
+                             else process.communicate(stream_in.read().encode('utf-8')))]
     os.chdir(initialPath)
-    out   = out.decode('utf-8').replace('\x08','')
+    stdout = stdout.decode('utf-8').replace('\x08','')
-    error = error.decode('utf-8').replace('\x08','')
+    stderr = stderr.decode('utf-8').replace('\x08','')
     if process.returncode != 0:
         raise RuntimeError('{} failed with returncode {}'.format(cmd,process.returncode))
-    return out,error
+    return stdout, stderr
 
 
 def show_progress(iterable,N_iter=None,prefix='',bar_length=50):

python/tests/test_Result.py

@@ -302,8 +302,8 @@ class TestResult:
 
     @pytest.mark.parametrize('allowed',['off','on'])
     def test_rename(self,default,allowed):
-        F = default.read_dataset(default.get_dataset_location('F'))
         if allowed == 'on':
+            F = default.read_dataset(default.get_dataset_location('F'))
             default.allow_modification()
             default.rename('F','new_name')
             assert np.all(F == default.read_dataset(default.get_dataset_location('new_name')))

python/tests/test_Rotation.py

@@ -14,7 +14,7 @@ scatter=1.e-2
 
 @pytest.fixture
 def default():
-    """A set of n random rotations."""
+    """A set of n rotations (corner cases and random)."""
     specials = np.array([
                          [1.0, 0.0, 0.0, 0.0],
                         #----------------------
@@ -170,7 +170,7 @@ def qu2ax(qu):
 
     Modified version of the original formulation, should be numerically more stable
     """
-    if np.isclose(qu[0],1.,rtol=0.0):                                                              # set axis to [001] if the angle is 0/360
+    if np.isclose(qu[0],1.,rtol=0.0):                                                               # set axis to [001] if the angle is 0/360
         ax = np.array([ 0.0, 0.0, 1.0, 0.0 ])
     elif qu[0] > 1.e-8:
         s = np.sign(qu[0])/np.sqrt(qu[1]**2+qu[2]**2+qu[3]**2)
@@ -534,7 +534,7 @@ def mul(me, other):
         other_p = other.quaternion[1:]
         R = me.__class__(np.append(me_q*other_q - np.dot(me_p,other_p),
                                    me_q*other_p + other_q*me_p + _P * np.cross(me_p,other_p)))
-        return R.standardize()
+        return R._standardize()
     elif isinstance(other, np.ndarray):
         if other.shape == (3,):
             A = me.quaternion[0]**2.0 - np.dot(me.quaternion[1:],me.quaternion[1:])
@@ -554,7 +554,7 @@ def mul(me, other):
             axes = ((0, 2, 4, 6), (0, 1, 2, 3))
             return np.tensordot(RRRR, other, axes)
         else:
-            raise ValueError('Can only rotate vectors, 2nd order ternsors, and 4th order tensors')
+            raise ValueError('Can only rotate vectors, 2nd order tensors, and 4th order tensors')
     else:
         raise TypeError('Cannot rotate {}'.format(type(other)))
 
@@ -854,6 +854,10 @@ class TestRotation:
             rot = Rotation.from_basis(om,False,reciprocal=reciprocal)
             assert np.isclose(np.linalg.det(rot.as_matrix()),1.0)
 
+    @pytest.mark.parametrize('shape',[None,1,(4,4)])
+    def test_random(self,shape):
+        Rotation.from_random(shape)
+
     @pytest.mark.parametrize('function',[Rotation.from_quaternion,
                                          Rotation.from_Eulers,
                                          Rotation.from_axis_angle,
@@ -866,6 +870,16 @@ class TestRotation:
         with pytest.raises(ValueError):
             function(invalid_shape)
 
+    @pytest.mark.parametrize('fr,to',[(Rotation.from_quaternion,'as_quaternion'),
+                                      (Rotation.from_axis_angle,'as_axis_angle'),
+                                      (Rotation.from_Rodrigues, 'as_Rodrigues'),
+                                      (Rotation.from_homochoric,'as_homochoric'),
+                                      (Rotation.from_cubochoric,'as_cubochoric')])
+    def test_invalid_P(self,fr,to):
+        R = Rotation.from_random(np.random.randint(8,32,(3)))                                       # noqa
+        with pytest.raises(ValueError):
+            fr(eval('R.{}()'.format(to)),P=-30)
+
     @pytest.mark.parametrize('shape',[None,(3,),(4,2)])
     def test_broadcast(self,shape):
         rot = Rotation.from_random(shape)
@@ -932,14 +946,18 @@ class TestRotation:
         phi_2 = 2*np.pi - phi_1
         R_1 = Rotation.from_Eulers(np.array([phi_1,0.,0.]))
         R_2 = Rotation.from_Eulers(np.array([0.,0.,phi_2]))
-        assert np.allclose(data,R_2*(R_1*data))
+        assert np.allclose(data,R_2@(R_1@data))
+
+    def test_rotate_inverse(self):
+        R = Rotation.from_random()
+        assert np.allclose(np.eye(3),(~R@R).as_matrix())
 
     @pytest.mark.parametrize('data',[np.random.rand(3),
                                      np.random.rand(3,3),
                                      np.random.rand(3,3,3,3)])
-    def test_rotate_inverse(self,data):
+    def test_rotate_inverse_array(self,data):
         R = Rotation.from_random()
-        assert np.allclose(data,R.inversed()*(R*data))
+        assert np.allclose(data,~R@(R@data))
 
     @pytest.mark.parametrize('data',[np.random.rand(4),
                                      np.random.rand(3,2),
@@ -956,3 +974,7 @@ class TestRotation:
         R = Rotation.from_random()
         with pytest.raises(TypeError):
             R*data
+
+    def test_misorientation(self):
+        R = Rotation.from_random()
+        assert np.allclose(R.misorientation(R).as_matrix(),np.eye(3))

python/tests/test_util.py

@@ -0,0 +1,11 @@
+from damask import util
+
+class TestUtil:
+
+    def test_execute_direct(self):
+        out,err = util.execute('echo test')
+        assert out=='test\n' and err==''
+    
+    def test_execute_env(self):
+        out,err = util.execute('sh -c "echo $test_for_execute"',env={'test_for_execute':'test'})
+        assert out=='test\n' and err==''

src/constitutive_plastic_dislotwin.f90

@@ -830,7 +830,7 @@ module subroutine plastic_dislotwin_results(instance,group)
                                                      'mobile dislocation density','1/m²')
       case('rho_dip')
         if(prm%sum_N_sl>0) call results_writeDataset(group,stt%rho_dip,'rho_dip',&
-                                                     'dislocation dipole density''1/m²')
+                                                     'dislocation dipole density','1/m²')
       case('gamma_sl')
         if(prm%sum_N_sl>0) call results_writeDataset(group,stt%gamma_sl,'gamma_sl',&
                                                      'plastic shear','1')

src/crystallite.f90

@@ -284,11 +284,9 @@ end subroutine crystallite_init
 !--------------------------------------------------------------------------------------------------
 !> @brief calculate stress (P)
 !--------------------------------------------------------------------------------------------------
-function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
+function crystallite_stress()
 
   logical, dimension(discretization_nIP,discretization_nElem) :: crystallite_stress
-  real(pReal), intent(in), optional :: &
-    dummyArgumentToPreventInternalCompilerErrorWithGCC
   real(pReal) :: &
     formerSubStep
   integer :: &

src/discretization.f90

@@ -49,7 +49,7 @@ subroutine discretization_init(homogenizationAt,microstructureAt,&
     IPcoords0, &
     NodeCoords0
   integer, optional,           intent(in) :: &
-    sharedNodesBegin
+    sharedNodesBegin                                                                                !< index of first node shared among different processes (MPI)
 
   write(6,'(/,a)') ' <<<+-  discretization init  -+>>>'; flush(6)
 

src/mesh/discretization_mesh.f90

@@ -8,53 +8,50 @@ module discretization_mesh
 #include <petsc/finclude/petscdmplex.h>
 #include <petsc/finclude/petscis.h>
 #include <petsc/finclude/petscdmda.h>
- use PETScdmplex
+  use PETScdmplex
- use PETScdmda
+  use PETScdmda
- use PETScis
+  use PETScis
 
- use DAMASK_interface
+  use DAMASK_interface
- use IO
+  use IO
- use debug
+  use debug
- use discretization
+  use discretization
- use numerics
+  use numerics
- use FEsolving
+  use FEsolving
- use FEM_quadrature
+  use FEM_quadrature
- use prec
+  use prec
- 
- implicit none
- private
- 
- integer, public, protected :: &
-   mesh_Nboundaries, &
-   mesh_NcpElemsGlobal
 
- integer :: &
+  implicit none
-   mesh_NcpElems                                                                                     !< total number of CP elements in mesh
+  private
+
+  integer, public, protected :: &
+    mesh_Nboundaries, &
+    mesh_NcpElemsGlobal
+
+  integer :: &
+    mesh_NcpElems                                                                                   !< total number of CP elements in mesh
 
 !!!! BEGIN DEPRECATED !!!!!
- integer, public, protected :: &
+  integer, public, protected :: &
-   mesh_maxNips                                                                                     !< max number of IPs in any CP element
+    mesh_maxNips                                                                                    !< max number of IPs in any CP element
 !!!! BEGIN DEPRECATED !!!!!
 
- integer, dimension(:,:), allocatable :: &
+  real(pReal), dimension(:,:), allocatable :: &
-   mesh_element !DEPRECATED
+    mesh_ipVolume, &                                                                                !< volume associated with IP (initially!)
+    mesh_node0                                                                                      !< node x,y,z coordinates (initially!)
 
- real(pReal), dimension(:,:), allocatable :: &
+  real(pReal), dimension(:,:,:), allocatable :: &
-   mesh_ipVolume, &                                                                                 !< volume associated with IP (initially!)
+    mesh_ipCoordinates                                                                              !< IP x,y,z coordinates (after deformation!)
-   mesh_node0                                                                                       !< node x,y,z coordinates (initially!)
- 
- real(pReal), dimension(:,:,:), allocatable :: &
-   mesh_ipCoordinates                                                                               !< IP x,y,z coordinates (after deformation!)
 
- DM, public :: geomMesh
+  DM, public :: geomMesh
- 
- PetscInt, dimension(:), allocatable, public, protected :: &
-   mesh_boundaries
 
- public :: &
+  PetscInt, dimension(:), allocatable, public, protected :: &
-   discretization_mesh_init, &
+    mesh_boundaries
-   mesh_FEM_build_ipVolumes, &
+
-   mesh_FEM_build_ipCoordinates
+  public :: &
+    discretization_mesh_init, &
+    mesh_FEM_build_ipVolumes, &
+    mesh_FEM_build_ipCoordinates
 
 contains
 
@@ -67,39 +64,33 @@ subroutine discretization_mesh_init(restart)
 
   logical, intent(in) :: restart
 
-  integer, dimension(1), parameter:: FE_geomtype = [1]                                              !< geometry type of particular element type
-  integer, dimension(1) :: FE_Nips                                                                  !< number of IPs in a specific type of element
-  
   integer, allocatable, dimension(:) :: chunkPos
   integer :: dimPlex, &
     mesh_Nnodes, &                                                                                  !< total number of nodes in mesh
     j, l
-  integer, parameter :: &
-    mesh_ElemType=1                                                                                 !< Element type of the mesh (only support homogeneous meshes)
   PetscSF :: sf
   DM :: globalMesh
   PetscInt :: nFaceSets
   PetscInt, pointer, dimension(:) :: pFaceSets
   character(len=pStringLen), dimension(:), allocatable :: fileContent
-  IS :: faceSetIS 
+  IS :: faceSetIS
   PetscErrorCode :: ierr
+  integer, dimension(:), allocatable :: &
+    homogenizationAt, &
+    microstructureAt
+
 
-  
   write(6,'(/,a)')   ' <<<+-  mesh init  -+>>>'
 
-  ! read in file
   call DMPlexCreateFromFile(PETSC_COMM_WORLD,geometryFile,PETSC_TRUE,globalMesh,ierr)
   CHKERRQ(ierr)
-  ! get spatial dimension (2 or 3?)
   call DMGetDimension(globalMesh,dimPlex,ierr)
   CHKERRQ(ierr)
-  write(6,*) 'dimension',dimPlex;flush(6)
   call DMGetStratumSize(globalMesh,'depth',dimPlex,mesh_NcpElemsGlobal,ierr)
   CHKERRQ(ierr)
   ! get number of IDs in face sets (for boundary conditions?)
   call DMGetLabelSize(globalMesh,'Face Sets',mesh_Nboundaries,ierr)
   CHKERRQ(ierr)
-  write(6,*) 'number of "Face Sets"',mesh_Nboundaries;flush(6)
   call MPI_Bcast(mesh_Nboundaries,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
   call MPI_Bcast(mesh_NcpElemsGlobal,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
   call MPI_Bcast(dimPlex,1,MPI_INTEGER,0,PETSC_COMM_WORLD,ierr)
@@ -142,48 +133,42 @@ subroutine discretization_mesh_init(restart)
     enddo
     call DMClone(globalMesh,geomMesh,ierr)
     CHKERRQ(ierr)
-  else 
+  else
     call DMPlexDistribute(globalMesh,0,sf,geomMesh,ierr)
     CHKERRQ(ierr)
-  endif  
+  endif
 
   call DMDestroy(globalMesh,ierr); CHKERRQ(ierr)
-  
+
   call DMGetStratumSize(geomMesh,'depth',dimPlex,mesh_NcpElems,ierr)
   CHKERRQ(ierr)
   call DMGetStratumSize(geomMesh,'depth',0,mesh_Nnodes,ierr)
   CHKERRQ(ierr)
 
-  FE_Nips(FE_geomtype(1)) = FEM_nQuadrature(dimPlex,integrationOrder)
+  mesh_maxNips = FEM_nQuadrature(dimPlex,integrationOrder)
-  mesh_maxNips = FE_Nips(1)
+
-  
-  write(6,*) 'mesh_maxNips',mesh_maxNips
   call mesh_FEM_build_ipCoordinates(dimPlex,FEM_quadrature_points(dimPlex,integrationOrder)%p)
   call mesh_FEM_build_ipVolumes(dimPlex)
-  
-  allocate (mesh_element (4,mesh_NcpElems)); mesh_element = 0
-  do j = 1, mesh_NcpElems
-    mesh_element( 1,j) = -1                                                                     ! DEPRECATED
-    mesh_element( 2,j) = mesh_elemType                                                               ! elem type
-    mesh_element( 3,j) = 1                                                                      ! homogenization
-    call DMGetLabelValue(geomMesh,'material',j-1,mesh_element(4,j),ierr)
-    CHKERRQ(ierr)
-  end do 
 
-  if (debug_e < 1 .or. debug_e > mesh_NcpElems) &
+  allocate(microstructureAt(mesh_NcpElems))
-    call IO_error(602,ext_msg='element')                                                        ! selected element does not exist
+  allocate(homogenizationAt(mesh_NcpElems),source=1)
-  if (debug_i < 1 .or. debug_i > FE_Nips(FE_geomtype(mesh_element(2,debug_e)))) &
+  do j = 1, mesh_NcpElems
-    call IO_error(602,ext_msg='IP')                                                             ! selected element does not have requested IP
+    call DMGetLabelValue(geomMesh,'material',j-1,microstructureAt(j),ierr)
-  
+    CHKERRQ(ierr)
-  FEsolving_execElem = [ 1,mesh_NcpElems ]                                                      ! parallel loop bounds set to comprise all DAMASK elements
+  end do
-  FEsolving_execIP   = [1,FE_Nips(FE_geomtype(mesh_element(2,1)))]
+
-  
+  if (debug_e < 1 .or. debug_e > mesh_NcpElems) call IO_error(602,ext_msg='element')
+  if (debug_i < 1 .or. debug_i > mesh_maxNips)  call IO_error(602,ext_msg='IP')
+
+  FEsolving_execElem = [1,mesh_NcpElems]                                                            ! parallel loop bounds set to comprise all DAMASK elements
+  FEsolving_execIP   = [1,mesh_maxNips]
+
   allocate(mesh_node0(3,mesh_Nnodes),source=0.0_pReal)
 
-  call discretization_init(mesh_element(3,:),mesh_element(4,:),&
+  call discretization_init(microstructureAt,homogenizationAt,&
                            reshape(mesh_ipCoordinates,[3,mesh_maxNips*mesh_NcpElems]), &
                            mesh_node0)
- 
+
 end subroutine discretization_mesh_init
 
 
@@ -191,23 +176,23 @@ end subroutine discretization_mesh_init
 !> @brief Calculates IP volume. Allocates global array 'mesh_ipVolume'
 !--------------------------------------------------------------------------------------------------
 subroutine mesh_FEM_build_ipVolumes(dimPlex)
- 
+
   PetscInt,intent(in):: dimPlex
   PetscReal          :: vol
   PetscReal, pointer,dimension(:) :: pCent, pNorm
   PetscInt           :: cellStart, cellEnd, cell
   PetscErrorCode     :: ierr
- 
+
   allocate(mesh_ipVolume(mesh_maxNips,mesh_NcpElems),source=0.0_pReal)
- 
+
   call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
   allocate(pCent(dimPlex))
   allocate(pNorm(dimPlex))
   do cell = cellStart, cellEnd-1
-    call  DMPlexComputeCellGeometryFVM(geomMesh,cell,vol,pCent,pNorm,ierr) 
+    call  DMPlexComputeCellGeometryFVM(geomMesh,cell,vol,pCent,pNorm,ierr)
     CHKERRQ(ierr)
     mesh_ipVolume(:,cell+1) = vol/real(mesh_maxNips,pReal)
-  enddo  
+  enddo
 
 end subroutine mesh_FEM_build_ipVolumes
 
@@ -219,20 +204,20 @@ subroutine mesh_FEM_build_ipCoordinates(dimPlex,qPoints)
 
   PetscInt,      intent(in) :: dimPlex
   PetscReal,     intent(in) :: qPoints(mesh_maxNips*dimPlex)
-  
+
   PetscReal,        pointer,dimension(:) :: pV0, pCellJ, pInvcellJ
   PetscReal                 :: detJ
   PetscInt                  :: cellStart, cellEnd, cell, qPt, dirI, dirJ, qOffset
   PetscErrorCode            :: ierr
- 
+
- 
+
   allocate(mesh_ipCoordinates(3,mesh_maxNips,mesh_NcpElems),source=0.0_pReal)
- 
+
   allocate(pV0(dimPlex))
   allocatE(pCellJ(dimPlex**2))
   allocatE(pinvCellJ(dimPlex**2))
   call DMPlexGetHeightStratum(geomMesh,0,cellStart,cellEnd,ierr); CHKERRQ(ierr)
-  do cell = cellStart, cellEnd-1                                                                     !< loop over all elements 
+  do cell = cellStart, cellEnd-1                                                                     !< loop over all elements
     call DMPlexComputeCellGeometryAffineFEM(geomMesh,cell,pV0,pCellJ,pInvcellJ,detJ,ierr)
     CHKERRQ(ierr)
     qOffset = 0
@@ -246,7 +231,7 @@ subroutine mesh_FEM_build_ipCoordinates(dimPlex,qPoints)
       enddo
       qOffset = qOffset + dimPlex
     enddo
-  enddo  
+  enddo
 
 end subroutine mesh_FEM_build_ipCoordinates
 

src/quaternions.f90

@@ -101,6 +101,8 @@ module quaternions
     assignment(=), &
     conjg, aimag, &
     log, exp, &
+    abs, dot_product, &
+    inverse, &
     real
 
 contains