Merge remote-tracking branch 'origin/development' into new-gmsh-version

LICENSE

@@ -1,4 +1,4 @@
-Copyright 2011-20 Max-Planck-Institut für Eisenforschung GmbH
+Copyright 2011-21 Max-Planck-Institut für Eisenforschung GmbH
 
 DAMASK is free software: you can redistribute it and/or modify
 it under the terms of the GNU General Public License as published by

PRIVATE

@@ -1 +1 @@
-Subproject commit f6fd3227ec0f9c03fbf991bf7f8732b22ae96530
+Subproject commit b1a31a79cc90d458494068a96cfd3e9497aa330c

VERSION

@@ -1 +1 @@
-v3.0.0-alpha2-230-g0fc670d01
+v3.0.0-alpha2-279-g8182c9c54

python/damask/_config.py

@@ -1,3 +1,4 @@
+import copy
 from io import StringIO
 import abc
 
@@ -35,6 +36,14 @@ class Config(dict):
         output.seek(0)
         return ''.join(output.readlines())
 
+
+    def __copy__(self):
+        """Create deep copy."""
+        return copy.deepcopy(self)
+
+    copy = __copy__
+
+
     @classmethod
     def load(cls,fname):
         """
@@ -52,6 +61,7 @@ class Config(dict):
             fhandle = fname
         return cls(yaml.safe_load(fhandle))
 
+
     def save(self,fname,**kwargs):
         """
         Save to yaml file.
@@ -89,12 +99,37 @@ class Config(dict):
             fhandle.write(yaml.dump(self,Dumper=NiceDumper,**kwargs))
 
 
+    def add(self,d):
+        """
+        Add dictionary.
+
+        d : dict
+            Dictionary to append.
+        """
+        duplicate = self.copy()
+        duplicate.update(d)
+        return duplicate
+
+
+    def delete(self,key):
+        """
+        Delete item.
+
+        key : str or scalar
+            Label of the key to remove.
+        """
+        duplicate = self.copy()
+        del duplicate[key]
+        return duplicate
+
+
     @property
     @abc.abstractmethod
     def is_complete(self):
         """Check for completeness."""
         pass
 
+
     @property
     @abc.abstractmethod
     def is_valid(self):

python/damask/_configmaterial.py

@@ -1,5 +1,3 @@
-import copy
-
 import numpy as np
 
 from . import Config
@@ -13,11 +11,10 @@ class ConfigMaterial(Config):
                  'homogenization': {},
                  'phase': {}}
 
-    def __init__(self,d={}):
+    def __init__(self,d=_defaults):
         """Initialize object with default dictionary keys."""
         super().__init__(d)
-        for k,v in self._defaults.items():
+
-            if k not in self: self[k] = v
 
     def save(self,fname='material.yaml',**kwargs):
         """
@@ -204,7 +201,7 @@ class ConfigMaterial(Config):
             Limit renaming to selected constituents.
 
         """
-        dup = copy.deepcopy(self)
+        dup = self.copy()
         for i,m in enumerate(dup['material']):
             if ID and i not in ID: continue
             for c in m['constituents']:
@@ -228,7 +225,7 @@ class ConfigMaterial(Config):
             Limit renaming to selected homogenization IDs.
 
         """
-        dup = copy.deepcopy(self)
+        dup = self.copy()
         for i,m in enumerate(dup['material']):
             if ID and i not in ID: continue
             try:
@@ -289,7 +286,7 @@ class ConfigMaterial(Config):
         c = [{} for _ in range(length)] if constituents is None else \
             [{'constituents':u} for u in ConfigMaterial._constituents(**constituents)]
 
-        if len(c) == 1: c = [copy.deepcopy(c[0]) for _ in range(length)]
+        if len(c) == 1: c = [c[0] for _ in range(length)]
 
         if length != 1 and length != len(c):
             raise ValueError('Cannot add entries of different length')
@@ -301,7 +298,7 @@ class ConfigMaterial(Config):
             else:
                 for i in range(len(c)):
                     c[i][k] = v
-        dup = copy.deepcopy(self)
+        dup = self.copy()
         dup['material'] = dup['material'] + c if 'material' in dup else c
 
         return dup

python/damask/_grid.py

@@ -57,13 +57,10 @@ class Grid:
 
 
     def __copy__(self):
-        """Copy grid."""
+        """Create deep copy."""
         return copy.deepcopy(self)
 
-
+    copy = __copy__
-    def copy(self):
-        """Copy grid."""
-        return self.__copy__()
 
 
     def diff(self,other):
@@ -766,24 +763,19 @@ class Grid:
         if fill is None: fill = np.nanmax(self.material) + 1
         dtype = float if np.isnan(fill) or int(fill) != fill or self.material.dtype==np.float else int
 
-        Eulers = R.as_Euler_angles(degrees=True)
+        material = self.material
-        material_in = self.material.copy()
-
         # These rotations are always applied in the reference coordinate system, i.e. (z,x,z) not (z,x',z'')
         # see https://www.cs.utexas.edu/~theshark/courses/cs354/lectures/cs354-14.pdf
-        for angle,axes in zip(Eulers[::-1], [(0,1),(1,2),(0,1)]):
+        for angle,axes in zip(R.as_Euler_angles(degrees=True)[::-1], [(0,1),(1,2),(0,1)]):
-            material_out = ndimage.rotate(material_in,angle,axes,order=0,
+            material_temp = ndimage.rotate(material,angle,axes,order=0,prefilter=False,output=dtype,cval=fill)
-                                                prefilter=False,output=dtype,cval=fill)
-            if np.prod(material_in.shape) == np.prod(material_out.shape):
             # avoid scipy interpolation errors for rotations close to multiples of 90°
-                material_in = np.rot90(material_in,k=np.rint(angle/90.).astype(int),axes=axes)
+            material = material_temp if np.prod(material_temp.shape) != np.prod(material.shape) else \
-            else:
+                       np.rot90(material,k=np.rint(angle/90.).astype(int),axes=axes)
-                material_in = material_out
 
-        origin = self.origin-(np.asarray(material_in.shape)-self.cells)*.5 * self.size/self.cells
+        origin = self.origin-(np.asarray(material.shape)-self.cells)*.5 * self.size/self.cells
 
-        return Grid(material = material_in,
+        return Grid(material = material,
-                    size     = self.size/self.cells*np.asarray(material_in.shape),
+                    size     = self.size/self.cells*np.asarray(material.shape),
                     origin   = origin,
                     comments = self.comments+[util.execution_stamp('Grid','rotate')],
                    )

python/damask/_orientation.py

@@ -199,7 +199,7 @@ class Orientation(Rotation):
 
 
     def __copy__(self,**kwargs):
-        """Copy."""
+        """Create deep copy."""
         return self.__class__(rotation=kwargs['rotation'] if 'rotation' in kwargs else self.quaternion,
                               lattice =kwargs['lattice']  if 'lattice'  in kwargs else self.lattice
                                                                                     if self.lattice is not None else self.family,
@@ -225,30 +225,42 @@ class Orientation(Rotation):
             Orientation to check for equality.
 
         """
-        return  super().__eq__(other) \
+        matching_type = all([hasattr(other,attr) and getattr(self,attr) == getattr(other,attr)
-            and hasattr(other, 'family')     and self.family     == other.family \
+                             for attr in ['family','lattice','parameters']])
-            and hasattr(other, 'lattice')    and self.lattice    == other.lattice \
+        return np.logical_and(super().__eq__(other),matching_type)
-            and hasattr(other, 'parameters') and self.parameters == other.parameters
 
-
+    def __ne__(self,other):
-    def __matmul__(self,other):
         """
-        Rotation of vector, second or fourth order tensor, or rotation object.
+        Not equal to other.
 
         Parameters
         ----------
-        other : numpy.ndarray, Rotation, or Orientation
+        other : Orientation
-            Vector, second or fourth order tensor, or rotation object that is rotated.
+            Orientation to check for equality.
+
+        """
+        return np.logical_not(self==other)
+
+
+    def __mul__(self,other):
+        """
+        Compose this orientation with other.
+
+        Parameters
+        ----------
+        other : Rotation or Orientation
+            Object for composition.
 
         Returns
         -------
-        other_rot : numpy.ndarray or Rotation
+        composition : Orientation
-            Rotated vector, second or fourth order tensor, or rotation object.
+            Compound rotation self*other, i.e. first other then self rotation.
 
         """
-        return self.copy(rotation=Rotation.__matmul__(self,Rotation(other.quaternion))) \
+        if isinstance(other,Orientation) or isinstance(other,Rotation):
-               if isinstance(other,self.__class__) else \
+            return self.copy(rotation=Rotation.__mul__(self,Rotation(other.quaternion)))
-               Rotation.__matmul__(self,other)
+        else:
+            raise TypeError('Use "O@b", i.e. matmul, to apply Orientation "O" to object "b"')
 
 
     @classmethod
@@ -429,7 +441,7 @@ class Orientation(Rotation):
             raise ValueError('Missing crystal symmetry')
 
         o = self.symmetry_operations.broadcast_to(self.symmetry_operations.shape+self.shape,mode='right')
-        return self.copy(rotation=o@Rotation(self.quaternion).broadcast_to(o.shape,mode='left'))
+        return self.copy(rotation=o*Rotation(self.quaternion).broadcast_to(o.shape,mode='left'))
 
 
     @property
@@ -608,7 +620,7 @@ class Orientation(Rotation):
         o,lattice = self.relation_operations(model,return_lattice=True)
         target = Orientation(lattice=lattice)
         o = o.broadcast_to(o.shape+self.shape,mode='right')
-        return self.copy(rotation=o@Rotation(self.quaternion).broadcast_to(o.shape,mode='left'),
+        return self.copy(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'],
                          c = self.c if target.ratio['c'] is None else self.a*target.ratio['c'],

python/damask/_rotation.py

@@ -35,6 +35,11 @@ class Rotation:
     - b = Q @ a
     - b = np.dot(Q.as_matrix(),a)
 
+    Compound rotations R1 (first) and R2 (second):
+
+    - R = R2 * R1
+    - R = Rotation.from_matrix(np.dot(R2.as_matrix(),R1.as_matrix())
+
     References
     ----------
     D. Rowenhorst et al., Modelling and Simulation in Materials Science and Engineering 23:083501, 2015
@@ -65,22 +70,13 @@ class Rotation:
 
 
     def __repr__(self):
-        """Represent rotation as unit quaternion, rotation matrix, and Bunge-Euler angles."""
+        """Represent rotation as unit quaternion(s)."""
-        if self == Rotation():
+        return f'Quaternion{" " if self.quaternion.shape == (4,) else "s of shape "+str(self.quaternion.shape)+chr(10)}'\
-           return 'Rotation()'
+               + str(self.quaternion)
-        else:
-            return f'Quaternions {self.shape}:\n'+str(self.quaternion) \
-                   if self.quaternion.shape != (4,) else \
-                   '\n'.join([
-                   'Quaternion: (real={:.3f}, imag=<{:+.3f}, {:+.3f}, {:+.3f}>)'.format(*(self.quaternion)),
-                   'Matrix:\n{}'.format(np.round(self.as_matrix(),8)),
-                   'Bunge Eulers / deg: ({:3.2f}, {:3.2f}, {:3.2f})'.format(*self.as_Euler_angles(degrees=True)),
-                    ])
 
 
-    # ToDo: Check difference __copy__ vs __deepcopy__
     def __copy__(self,**kwargs):
-        """Copy."""
+        """Create deep copy."""
         return self.__class__(rotation=kwargs['rotation'] if 'rotation' in kwargs else self.quaternion)
 
     copy = __copy__
@@ -97,6 +93,26 @@ class Rotation:
         """
         Equal to other.
 
+        Equality is determined taking limited floating point precision into account.
+        See numpy.allclose for details.
+
+        Parameters
+        ----------
+        other : Rotation
+            Rotation to check for equality.
+
+        """
+        s = self.quaternion
+        o = other.quaternion
+        if self.shape == () == other.shape:
+            return np.allclose(s,o) or (np.isclose(s[0],0.0) and np.allclose(s,-1.0*o))
+        else:
+            return np.all(np.isclose(s,o),-1) + np.all(np.isclose(s,-1.0*o),-1) * np.isclose(s[...,0],0.0)
+
+    def __ne__(self,other):
+        """
+        Not equal to other.
+
         Equality is determined taking limited floating point precision into
         account. See numpy.allclose for details.
 
@@ -106,8 +122,7 @@ class Rotation:
             Rotation to check for equality.
 
         """
-        return      np.prod(self.shape,dtype=int) == np.prod(other.shape,dtype=int) \
+        return np.logical_not(self==other)
-                and np.allclose(self.quaternion,other.quaternion)
 
 
     @property
@@ -127,41 +142,46 @@ class Rotation:
         return dup
 
 
-    def __pow__(self,pwr):
+    def __pow__(self,exp):
         """
-        Raise quaternion to power.
+        Perform the rotation 'exp' times.
-
-        Equivalent to performing the rotation 'pwr' times.
 
         Parameters
         ----------
-        pwr : float
+        exp : float
-            Power to raise quaternion to.
+            Exponent.
 
         """
         phi = np.arccos(self.quaternion[...,0:1])
         p = self.quaternion[...,1:]/np.linalg.norm(self.quaternion[...,1:],axis=-1,keepdims=True)
-        return self.copy(rotation=Rotation(np.block([np.cos(pwr*phi),np.sin(pwr*phi)*p]))._standardize())
+        return self.copy(rotation=Rotation(np.block([np.cos(exp*phi),np.sin(exp*phi)*p]))._standardize())
 
-
+    def __ipow__(self,exp):
-    def __mul__(self,other):
-        """Standard multiplication is not implemented."""
-        raise NotImplementedError('Use "R@b", i.e. matmul, to apply rotation "R" to object "b"')
-
-
-    def __matmul__(self,other):
         """
-        Rotation of vector, second or fourth order tensor, or rotation object.
+        Perform the rotation 'exp' times (in-place).
 
         Parameters
         ----------
-        other : numpy.ndarray or Rotation
+        exp : float
-            Vector, second or fourth order tensor, or rotation object that is rotated.
+            Exponent.
+
+        """
+        return self**exp
+
+
+    def __mul__(self,other):
+        """
+        Compose this rotation with other.
+
+        Parameters
+        ----------
+        other : Rotation of shape(self.shape)
+            Rotation for composition.
 
         Returns
         -------
-        other_rot : numpy.ndarray or Rotation
+        composition : Rotation
-            Rotated vector, second or fourth order tensor, or rotation object.
+            Compound rotation self*other, i.e. first other then self rotation.
 
         """
         if isinstance(other,Rotation):
@@ -172,8 +192,71 @@ class Rotation:
             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 Rotation(np.block([q,p]))._standardize()
+        else:
+            raise TypeError('Use "R@b", i.e. matmul, to apply rotation "R" to object "b"')
 
-        elif isinstance(other,np.ndarray):
+    def __imul__(self,other):
+        """
+        Compose this rotation with other (in-place).
+
+        Parameters
+        ----------
+        other : Rotation of shape(self.shape)
+            Rotation for composition.
+
+        """
+        return self*other
+
+
+    def __truediv__(self,other):
+        """
+        Compose this rotation with inverse of other.
+
+        Parameters
+        ----------
+        other : damask.Rotation of shape (self.shape)
+            Rotation to inverse composition.
+
+        Returns
+        -------
+        composition : Rotation
+            Compound rotation self*(~other), i.e. first inverse of other then self rotation.
+
+        """
+        if isinstance(other,Rotation):
+            return self*~other
+        else:
+            raise TypeError('Use "R@b", i.e. matmul, to apply rotation "R" to object "b"')
+
+    def __itruediv__(self,other):
+        """
+        Compose this rotation with inverse of other (in-place).
+
+        Parameters
+        ----------
+        other : Rotation of shape (self.shape)
+            Rotation to inverse composition.
+
+        """
+        return self/other
+
+
+    def __matmul__(self,other):
+        """
+        Rotation of vector, second order tensor, or fourth order tensor.
+
+        Parameters
+        ----------
+        other : numpy.ndarray of shape (...,3), (...,3,3), or (...,3,3,3,3)
+            Vector or tensor on which to apply the rotation.
+
+        Returns
+        -------
+        rotated : numpy.ndarray of shape (...,3), (...,3,3), or (...,3,3,3,3)
+            Rotated vector or tensor, i.e. transformed to frame defined by rotation.
+
+        """
+        if isinstance(other,np.ndarray):
             if self.shape + (3,) == other.shape:
                 q_m = self.quaternion[...,0]
                 p_m = self.quaternion[...,1:]
@@ -193,9 +276,13 @@ class Rotation:
                 return np.einsum('...im,...jn,...ko,...lp,...mnop',R,R,R,R,other)
             else:
                 raise ValueError('Can only rotate vectors, 2nd order tensors, and 4th order tensors')
+        elif isinstance(other,Rotation):
+            raise TypeError('Use "R1*R2", i.e. multiplication, to compose rotations "R1" and "R2"')
         else:
             raise TypeError(f'Cannot rotate {type(other)}')
 
+    apply = __matmul__
+
 
     def _standardize(self):
         """Standardize quaternion (ensure positive real hemisphere)."""
@@ -296,7 +383,7 @@ class Rotation:
             Rotation to which the misorientation is computed.
 
         """
-        return other@~self
+        return other*~self
 
 
     ################################################################################################
@@ -819,7 +906,7 @@ class Rotation:
                              np.sqrt(1-u**2)*np.sin(Theta),
                              u, omega])
 
-        return  Rotation.from_axis_angle(p) @ center
+        return  Rotation.from_axis_angle(p) * center
 
 
     @staticmethod
@@ -870,8 +957,8 @@ class Rotation:
         f[::2,:3] *= -1                                                                             # flip half the rotation axes to negative sense
 
         return R_align.broadcast_to(N) \
-             @ Rotation.from_axis_angle(p,normalize=True) \
+             * Rotation.from_axis_angle(p,normalize=True) \
-             @ Rotation.from_axis_angle(f)
+             * Rotation.from_axis_angle(f)
 
 
 ####################################################################################################
@@ -1060,7 +1147,6 @@ class Rotation:
     @staticmethod
     def _om2ax(om):
         """Rotation matrix to axis angle pair."""
-        #return Rotation._qu2ax(Rotation._om2qu(om)) # HOTFIX
         diag_delta = -_P*np.block([om[...,1,2:3]-om[...,2,1:2],
                                    om[...,2,0:1]-om[...,0,2:3],
                                    om[...,0,1:2]-om[...,1,0:1]

python/damask/_table.py

@@ -42,12 +42,10 @@ class Table:
         return len(self.data)
 
     def __copy__(self):
-        """Copy Table."""
+        """Create deep copy."""
         return copy.deepcopy(self)
 
-    def copy(self):
+    copy = __copy__
-        """Copy Table."""
-        return self.__copy__()
 
 
     def _label_discrete(self):

python/tests/test_Config.py

@@ -22,6 +22,10 @@ class TestConfig:
         with open(tmp_path/'config.yaml') as f:
             assert Config.load(f) == config
 
+    def test_add_remove(self):
+        config = Config()
+        assert config.add({'hello':'world'}).delete('hello') == config
+
     def test_repr(self,tmp_path):
         config = Config()
         config['A'] = 1

python/tests/test_Orientation.py

@@ -25,13 +25,16 @@ class TestOrientation:
     @pytest.mark.parametrize('shape',[None,5,(4,6)])
     def test_equal(self,lattice,shape):
         R = Rotation.from_random(shape)
-        assert Orientation(R,lattice) == Orientation(R,lattice)
+        assert Orientation(R,lattice) == Orientation(R,lattice) if shape is None else \
+              (Orientation(R,lattice) == Orientation(R,lattice)).all()
+
 
     @pytest.mark.parametrize('lattice',Orientation.crystal_families)
     @pytest.mark.parametrize('shape',[None,5,(4,6)])
     def test_unequal(self,lattice,shape):
         R = Rotation.from_random(shape)
-        assert not(Orientation(R,lattice) != Orientation(R,lattice))
+        assert not ( Orientation(R,lattice) != Orientation(R,lattice) if shape is None else \
+                    (Orientation(R,lattice) != Orientation(R,lattice)).any())
 
     @pytest.mark.parametrize('a,b',[
                                     (dict(rotation=[1,0,0,0]),
@@ -403,7 +406,7 @@ class TestOrientation:
     def test_relationship_vectorize(self,set_of_quaternions,lattice,model):
         r = Orientation(rotation=set_of_quaternions[:200].reshape((50,4,4)),lattice=lattice).related(model)
         for i in range(200):
-            assert r.reshape((-1,200))[:,i] == Orientation(set_of_quaternions[i],lattice).related(model)
+            assert (r.reshape((-1,200))[:,i] == Orientation(set_of_quaternions[i],lattice).related(model)).all()
 
     @pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch'])
     @pytest.mark.parametrize('lattice',['cF','cI'])

python/tests/test_Rotation.py

@@ -526,7 +526,7 @@ class TestRotation:
             o = backward(forward(m))
             u = np.array([np.pi*2,np.pi,np.pi*2])
             ok = np.allclose(m,o,atol=atol)
-            ok = ok or np.allclose(np.where(np.isclose(m,u),m-u,m),np.where(np.isclose(o,u),o-u,o),atol=atol)
+            ok |= np.allclose(np.where(np.isclose(m,u),m-u,m),np.where(np.isclose(o,u),o-u,o),atol=atol)
             if np.isclose(m[1],0.0,atol=atol) or np.isclose(m[1],np.pi,atol=atol):
                 sum_phi = np.unwrap([m[0]+m[2],o[0]+o[2]])
                 ok |= np.isclose(sum_phi[0],sum_phi[1],atol=atol)
@@ -550,19 +550,22 @@ class TestRotation:
             assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) and o[3]<=np.pi+1.e-9, f'{m},{o},{rot.as_quaternion()}'
 
     @pytest.mark.parametrize('forward,backward',[(Rotation._ro2qu,Rotation._qu2ro),
-                                                 #(Rotation._ro2om,Rotation._om2ro),
+                                                 (Rotation._ro2om,Rotation._om2ro),
-                                                 #(Rotation._ro2eu,Rotation._eu2ro),
+                                                 (Rotation._ro2eu,Rotation._eu2ro),
                                                  (Rotation._ro2ax,Rotation._ax2ro),
                                                  (Rotation._ro2ho,Rotation._ho2ro),
                                                  (Rotation._ro2cu,Rotation._cu2ro)])
     def test_Rodrigues_internal(self,set_of_rotations,forward,backward):
         """Ensure invariance of conversion from Rodrigues-Frank vector and back."""
-        cutoff = np.tan(np.pi*.5*(1.-1e-4))
+        cutoff = np.tan(np.pi*.5*(1.-1e-5))
         for rot in set_of_rotations:
             m = rot.as_Rodrigues_vector()
             o = backward(forward(m))
             ok = np.allclose(np.clip(m,None,cutoff),np.clip(o,None,cutoff),atol=atol)
-            ok = ok or np.isclose(m[3],0.0,atol=atol)
+            ok |= np.isclose(m[3],0.0,atol=atol)
+            if m[3] > cutoff:
+                ok |= np.allclose(m[:3],-1*o[:3])
+
             assert ok and np.isclose(np.linalg.norm(o[:3]),1.0), f'{m},{o},{rot.as_quaternion()}'
 
     @pytest.mark.parametrize('forward,backward',[(Rotation._ho2qu,Rotation._qu2ho),
@@ -592,7 +595,7 @@ class TestRotation:
             o = backward(forward(m))
             ok = np.allclose(m,o,atol=atol)
             if np.count_nonzero(np.isclose(np.abs(o),np.pi**(2./3.)*.5)):
-                ok = ok or np.allclose(m*-1.,o,atol=atol)
+                ok |= np.allclose(m*-1.,o,atol=atol)
             assert ok and np.max(np.abs(o)) < np.pi**(2./3.) * 0.5 + 1.e-9, f'{m},{o},{rot.as_quaternion()}'
 
     @pytest.mark.parametrize('vectorized, single',[(Rotation._qu2om,qu2om),
@@ -719,7 +722,7 @@ class TestRotation:
             o = Rotation.from_axis_angle(rot.as_axis_angle()).as_axis_angle()
             ok = np.allclose(m,o,atol=atol)
             if np.isclose(m[3],np.pi,atol=atol):
-                ok = ok or np.allclose(m*np.array([-1.,-1.,-1.,1.]),o,atol=atol)
+                ok |= np.allclose(m*np.array([-1.,-1.,-1.,1.]),o,atol=atol)
             assert ok and np.isclose(np.linalg.norm(o[:3]),1.0) \
                       and o[3]<=np.pi+1.e-9, f'{m},{o},{rot.as_quaternion()}'
 
@@ -740,7 +743,7 @@ class TestRotation:
             m = rot.as_Rodrigues_vector()
             o = Rotation.from_homochoric(rot.as_homochoric()*P*-1,P).as_Rodrigues_vector()
             ok = np.allclose(np.clip(m,None,cutoff),np.clip(o,None,cutoff),atol=atol)
-            ok = ok or np.isclose(m[3],0.0,atol=atol)
+            ok |= np.isclose(m[3],0.0,atol=atol)
             assert ok and np.isclose(np.linalg.norm(o[:3]),1.0), f'{m},{o},{rot.as_quaternion()}'
 
     @pytest.mark.parametrize('P',[1,-1])
@@ -780,8 +783,22 @@ class TestRotation:
         else:
             assert r.shape == shape
 
-    def test_equal(self):
+    @pytest.mark.parametrize('shape',[None,5,(4,6)])
-        assert Rotation.from_random(rng_seed=1) == Rotation.from_random(rng_seed=1)
+    def test_equal(self,shape):
+        R = Rotation.from_random(shape,rng_seed=1)
+        assert R == R if shape is None else (R == R).all()
+
+    @pytest.mark.parametrize('shape',[None,5,(4,6)])
+    def test_unequal(self,shape):
+        R = Rotation.from_random(shape,rng_seed=1)
+        assert not (R != R if shape is None else (R != R).any())
+
+
+    def test_equal_ambiguous(self):
+        qu = np.random.rand(10,4)
+        qu[:,0] = 0.
+        qu/=np.linalg.norm(qu,axis=1,keepdims=True)
+        assert (Rotation(qu) == Rotation(-qu)).all()
 
     def test_inversion(self):
         r = Rotation.from_random()
@@ -798,7 +815,7 @@ class TestRotation:
         p = Rotation.from_random(shape=shape)
         s = r.append(p)
         print(f'append 2x {shape} --> {s.shape}')
-        assert s[0,...] == r[0,...] and s[-1,...] == p[-1,...]
+        assert np.logical_and(s[0,...] == r[0,...], s[-1,...] == p[-1,...]).all()
 
     @pytest.mark.parametrize('shape',[None,1,(1,),(4,2),(3,3,2)])
     def test_append_list(self,shape):
@@ -806,7 +823,7 @@ class TestRotation:
         p = Rotation.from_random(shape=shape)
         s = r.append([r,p])
         print(f'append 3x {shape} --> {s.shape}')
-        assert s[0,...] == r[0,...] and s[-1,...] == p[-1,...]
+        assert np.logical_and(s[0,...] == r[0,...], s[-1,...] == p[-1,...]).all()
 
     @pytest.mark.parametrize('quat,standardized',[
                                                   ([-1,0,0,0],[1,0,0,0]),
@@ -828,7 +845,7 @@ class TestRotation:
     @pytest.mark.parametrize('order',['C','F'])
     def test_flatten_reshape(self,shape,order):
         r = Rotation.from_random(shape=shape)
-        assert r == r.flatten(order).reshape(shape,order)
+        assert (r == r.flatten(order).reshape(shape,order)).all()
 
     @pytest.mark.parametrize('function',[Rotation.from_quaternion,
                                          Rotation.from_Euler_angles,
@@ -939,7 +956,7 @@ class TestRotation:
 
     def test_rotate_inverse(self):
         R = Rotation.from_random()
-        assert np.allclose(np.eye(3),(~R@R).as_matrix())
+        assert np.allclose(np.eye(3),(~R*R).as_matrix())
 
     @pytest.mark.parametrize('data',[np.random.rand(3),
                                      np.random.rand(3,3),
@@ -973,6 +990,42 @@ class TestRotation:
         R_2 = Rotation.from_Euler_angles([360,0,0],degrees=True)
         assert np.allclose(R_1.misorientation(R_2).as_matrix(),np.eye(3))
 
+    def test_composition(self):
+        a,b = (Rotation.from_random(),Rotation.from_random())
+        c = a * b
+        a *= b
+        assert c == a
+
+    def test_composition_invalid(self):
+        with pytest.raises(TypeError):
+            Rotation()*np.ones(3)
+
+    def test_composition_inverse(self):
+        a,b = (Rotation.from_random(),Rotation.from_random())
+        c = a / b
+        a /= b
+        assert c == a
+
+    def test_composition_inverse_invalid(self):
+        with pytest.raises(TypeError):
+            Rotation()/np.ones(3)
+
+    def test_power(self):
+        a = Rotation.from_random()
+        r = (np.random.rand()-.5)*4
+        b = a**r
+        a **= r
+        assert a == b
+
+    def test_invariant(self):
+        R = Rotation.from_random()
+        assert R/R == R*R**(-1) == Rotation()
+
+    @pytest.mark.parametrize('item',[np.ones(3),np.ones((3,3)), np.ones((3,3,3,3))])
+    def test_apply(self,item):
+        r = Rotation.from_random()
+        assert (r.apply(item) == r@item).all()
+
     @pytest.mark.parametrize('angle',[10,20,30,40,50,60,70,80,90,100,120])
     def test_average(self,angle):
         R = Rotation.from_axis_angle([[0,0,1,10],[0,0,1,angle]],degrees=True)

src/C_routines.c

@@ -43,7 +43,7 @@ void gethostname_c(char hostname[], int *stat){
 
 
 void getusername_c(char username[], int *stat){
-  struct passwd *pw = getpwuid(geteuid());
+  struct passwd *pw = getpwuid(getuid());
   if(pw && strlen(pw->pw_name) <= STRLEN){
     strncpy(username,pw->pw_name,STRLEN+1);
     *stat = 0;

src/mesh/discretization_mesh.f90

@@ -114,11 +114,10 @@ subroutine discretization_mesh_init(restart)
 
   if (worldrank == 0) then
     call DMClone(globalMesh,geomMesh,ierr)
-    CHKERRQ(ierr)
   else
     call DMPlexDistribute(globalMesh,0,sf,geomMesh,ierr)
-    CHKERRQ(ierr)
   endif
+  CHKERRQ(ierr)
 
   allocate(mesh_boundaries(mesh_Nboundaries), source = 0)
   call DMGetLabelSize(globalMesh,'Face Sets',nFaceSets,ierr)
@@ -151,7 +150,7 @@ subroutine discretization_mesh_init(restart)
     call DMGetLabelValue(geomMesh,'Cell Sets',j-1,materialAt(j),ierr)
     CHKERRQ(ierr)
   end do
-  materialAt(:) = materialAt(:) + 1
+  materialAt = materialAt + 1
 
   if (debug_element < 1 .or. debug_element > mesh_NcpElems) call IO_error(602,ext_msg='element')
   if (debug_ip < 1 .or. debug_ip > mesh_maxNips)            call IO_error(602,ext_msg='IP')