diff --git a/.gitmodules b/.gitmodules index 0587fff4c..c415745bc 100644 --- a/.gitmodules +++ b/.gitmodules @@ -1,5 +1,5 @@ [submodule "PRIVATE"] - path = PRIVATE + path = PRIVATE url = ../PRIVATE.git branch = master - shallow = true + shallow = true diff --git a/LICENSE b/LICENSE index 3ffc3b9e3..4290d15bd 100644 --- a/LICENSE +++ b/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 diff --git a/PRIVATE b/PRIVATE index f6fd3227e..b1a31a79c 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit f6fd3227ec0f9c03fbf991bf7f8732b22ae96530 +Subproject commit b1a31a79cc90d458494068a96cfd3e9497aa330c diff --git a/VERSION b/VERSION index f123674b3..6e2dff56f 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v3.0.0-alpha2-230-g0fc670d01 +v3.0.0-alpha2-279-g8182c9c54 diff --git a/python/damask/_config.py b/python/damask/_config.py index 76955588f..9aa031ff0 100644 --- a/python/damask/_config.py +++ b/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): diff --git a/python/damask/_configmaterial.py b/python/damask/_configmaterial.py index b94e9897a..6415ee4dc 100644 --- a/python/damask/_configmaterial.py +++ b/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 diff --git a/python/damask/_grid.py b/python/damask/_grid.py index 8380bbc5b..0edec05f9 100644 --- a/python/damask/_grid.py +++ b/python/damask/_grid.py @@ -57,13 +57,10 @@ class Grid: def __copy__(self): - """Copy grid.""" + """Create deep copy.""" return copy.deepcopy(self) - - def copy(self): - """Copy grid.""" - return self.__copy__() + copy = __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_in = self.material.copy() - + material = self.material # 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)]): - material_out = ndimage.rotate(material_in,angle,axes,order=0, - 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) - else: - material_in = material_out + for angle,axes in zip(R.as_Euler_angles(degrees=True)[::-1], [(0,1),(1,2),(0,1)]): + material_temp = ndimage.rotate(material,angle,axes,order=0,prefilter=False,output=dtype,cval=fill) + # avoid scipy interpolation errors for rotations close to multiples of 90° + material = material_temp if np.prod(material_temp.shape) != np.prod(material.shape) else \ + np.rot90(material,k=np.rint(angle/90.).astype(int),axes=axes) - 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, - size = self.size/self.cells*np.asarray(material_in.shape), + return Grid(material = material, + size = self.size/self.cells*np.asarray(material.shape), origin = origin, comments = self.comments+[util.execution_stamp('Grid','rotate')], ) diff --git a/python/damask/_orientation.py b/python/damask/_orientation.py index e6434813e..cf31f4089 100644 --- a/python/damask/_orientation.py +++ b/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) \ - and hasattr(other, 'family') and self.family == other.family \ - and hasattr(other, 'lattice') and self.lattice == other.lattice \ - and hasattr(other, 'parameters') and self.parameters == other.parameters + matching_type = all([hasattr(other,attr) and getattr(self,attr) == getattr(other,attr) + for attr in ['family','lattice','parameters']]) + return np.logical_and(super().__eq__(other),matching_type) - - def __matmul__(self,other): + def __ne__(self,other): """ - Rotation of vector, second or fourth order tensor, or rotation object. + Not equal to other. Parameters ---------- - other : numpy.ndarray, Rotation, or Orientation - Vector, second or fourth order tensor, or rotation object that is rotated. + other : Orientation + 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 - Rotated vector, second or fourth order tensor, or rotation object. + composition : Orientation + Compound rotation self*other, i.e. first other then self rotation. """ - return self.copy(rotation=Rotation.__matmul__(self,Rotation(other.quaternion))) \ - if isinstance(other,self.__class__) else \ - Rotation.__matmul__(self,other) + if isinstance(other,Orientation) or isinstance(other,Rotation): + return self.copy(rotation=Rotation.__mul__(self,Rotation(other.quaternion))) + 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'], diff --git a/python/damask/_rotation.py b/python/damask/_rotation.py index 78029b59a..441fb5b01 100644 --- a/python/damask/_rotation.py +++ b/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.""" - if self == Rotation(): - return 'Rotation()' - 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)), - ]) + """Represent rotation as unit quaternion(s).""" + return f'Quaternion{" " if self.quaternion.shape == (4,) else "s of shape "+str(self.quaternion.shape)+chr(10)}'\ + + str(self.quaternion) - # 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) \ - and np.allclose(self.quaternion,other.quaternion) + return np.logical_not(self==other) @property @@ -127,41 +142,46 @@ class Rotation: return dup - def __pow__(self,pwr): + def __pow__(self,exp): """ - Raise quaternion to power. - - Equivalent to performing the rotation 'pwr' times. + Perform the rotation 'exp' times. Parameters ---------- - pwr : float - Power to raise quaternion to. + exp : float + 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 __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): + def __ipow__(self,exp): """ - Rotation of vector, second or fourth order tensor, or rotation object. + Perform the rotation 'exp' times (in-place). Parameters ---------- - other : numpy.ndarray or Rotation - Vector, second or fourth order tensor, or rotation object that is rotated. + exp : float + 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 - Rotated vector, second or fourth order tensor, or rotation object. + composition : Rotation + 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(f) + * Rotation.from_axis_angle(p,normalize=True) \ + * 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] diff --git a/python/damask/_table.py b/python/damask/_table.py index e6e6c4eeb..78a8a276e 100644 --- a/python/damask/_table.py +++ b/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 Table.""" - return self.__copy__() + copy = __copy__ def _label_discrete(self): diff --git a/python/tests/test_Config.py b/python/tests/test_Config.py index 67c419b3e..0319fb6de 100644 --- a/python/tests/test_Config.py +++ b/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 diff --git a/python/tests/test_Orientation.py b/python/tests/test_Orientation.py index 5ab0361a8..436b73c04 100644 --- a/python/tests/test_Orientation.py +++ b/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']) diff --git a/python/tests/test_Rotation.py b/python/tests/test_Rotation.py index 36e3a3ac9..6bee44e7f 100644 --- a/python/tests/test_Rotation.py +++ b/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._ro2eu,Rotation._eu2ro), + (Rotation._ro2om,Rotation._om2ro), + (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): - assert Rotation.from_random(rng_seed=1) == Rotation.from_random(rng_seed=1) + @pytest.mark.parametrize('shape',[None,5,(4,6)]) + 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) diff --git a/src/C_routines.c b/src/C_routines.c index 4b07c0ee0..3d62a87c2 100644 --- a/src/C_routines.c +++ b/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; diff --git a/src/mesh/discretization_mesh.f90 b/src/mesh/discretization_mesh.f90 index 1320d2609..4d8546b6a 100644 --- a/src/mesh/discretization_mesh.f90 +++ b/src/mesh/discretization_mesh.f90 @@ -82,7 +82,7 @@ subroutine discretization_mesh_init(restart) class(tNode), pointer :: & num_mesh integer :: integrationOrder !< order of quadrature rule required - + print'(/,a)', ' <<<+- discretization_mesh init -+>>>' @@ -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')