vectorized and cleaned

This commit is contained in:
Martin Diehl 2020-06-30 21:43:57 +02:00
parent 9a83b11a99
commit 49d448dced
8 changed files with 194 additions and 316 deletions

View File

@ -29,7 +29,7 @@ class Symmetry:
raise KeyError(f'Crystal system "{system}" is unknown') raise KeyError(f'Crystal system "{system}" is unknown')
self.system = system.lower() if isinstance(system,str) else system self.system = system.lower() if isinstance(system,str) else system
self.lattice = self.system # for compatibility self.lattice = self.system # ToDo: for compatibility
def __copy__(self): def __copy__(self):
@ -82,85 +82,10 @@ class Symmetry:
otherOrder = self.crystal_systems.index(other.system) otherOrder = self.crystal_systems.index(other.system)
return (myOrder > otherOrder) - (myOrder < otherOrder) return (myOrder > otherOrder) - (myOrder < otherOrder)
def symmetryOperations(self,members=[]):
"""List (or single element) of symmetry operations as rotations."""
if self.system == 'cubic':
symQuats = [
[ 1.0, 0.0, 0.0, 0.0 ],
[ 0.0, 1.0, 0.0, 0.0 ],
[ 0.0, 0.0, 1.0, 0.0 ],
[ 0.0, 0.0, 0.0, 1.0 ],
[ 0.0, 0.0, 0.5*np.sqrt(2), 0.5*np.sqrt(2) ],
[ 0.0, 0.0, 0.5*np.sqrt(2),-0.5*np.sqrt(2) ],
[ 0.0, 0.5*np.sqrt(2), 0.0, 0.5*np.sqrt(2) ],
[ 0.0, 0.5*np.sqrt(2), 0.0, -0.5*np.sqrt(2) ],
[ 0.0, 0.5*np.sqrt(2),-0.5*np.sqrt(2), 0.0 ],
[ 0.0, -0.5*np.sqrt(2),-0.5*np.sqrt(2), 0.0 ],
[ 0.5, 0.5, 0.5, 0.5 ],
[-0.5, 0.5, 0.5, 0.5 ],
[-0.5, 0.5, 0.5, -0.5 ],
[-0.5, 0.5, -0.5, 0.5 ],
[-0.5, -0.5, 0.5, 0.5 ],
[-0.5, -0.5, 0.5, -0.5 ],
[-0.5, -0.5, -0.5, 0.5 ],
[-0.5, 0.5, -0.5, -0.5 ],
[-0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ],
[ 0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ],
[-0.5*np.sqrt(2), 0.0, 0.5*np.sqrt(2), 0.0 ],
[-0.5*np.sqrt(2), 0.0, -0.5*np.sqrt(2), 0.0 ],
[-0.5*np.sqrt(2), 0.5*np.sqrt(2), 0.0, 0.0 ],
[-0.5*np.sqrt(2),-0.5*np.sqrt(2), 0.0, 0.0 ],
]
elif self.system == 'hexagonal':
symQuats = [
[ 1.0, 0.0, 0.0, 0.0 ],
[-0.5*np.sqrt(3), 0.0, 0.0, -0.5 ],
[ 0.5, 0.0, 0.0, 0.5*np.sqrt(3) ],
[ 0.0, 0.0, 0.0, 1.0 ],
[-0.5, 0.0, 0.0, 0.5*np.sqrt(3) ],
[-0.5*np.sqrt(3), 0.0, 0.0, 0.5 ],
[ 0.0, 1.0, 0.0, 0.0 ],
[ 0.0, -0.5*np.sqrt(3), 0.5, 0.0 ],
[ 0.0, 0.5, -0.5*np.sqrt(3), 0.0 ],
[ 0.0, 0.0, 1.0, 0.0 ],
[ 0.0, -0.5, -0.5*np.sqrt(3), 0.0 ],
[ 0.0, 0.5*np.sqrt(3), 0.5, 0.0 ],
]
elif self.system == 'tetragonal':
symQuats = [
[ 1.0, 0.0, 0.0, 0.0 ],
[ 0.0, 1.0, 0.0, 0.0 ],
[ 0.0, 0.0, 1.0, 0.0 ],
[ 0.0, 0.0, 0.0, 1.0 ],
[ 0.0, 0.5*np.sqrt(2), 0.5*np.sqrt(2), 0.0 ],
[ 0.0, -0.5*np.sqrt(2), 0.5*np.sqrt(2), 0.0 ],
[ 0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ],
[-0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ],
]
elif self.system == 'orthorhombic':
symQuats = [
[ 1.0,0.0,0.0,0.0 ],
[ 0.0,1.0,0.0,0.0 ],
[ 0.0,0.0,1.0,0.0 ],
[ 0.0,0.0,0.0,1.0 ],
]
else:
symQuats = [
[ 1.0,0.0,0.0,0.0 ],
]
symOps = list(map(Rotation,
np.array(symQuats)[np.atleast_1d(members) if members != [] else range(len(symQuats))]))
try:
iter(members) # asking for (even empty) list of members?
except TypeError:
return symOps[0] # no, return rotation object
else:
return symOps # yes, return list of rotations
@property @property
def symmetry_operations(self): def symmetry_operations(self):
"""Symmetry operations as Rotations.""" """Symmetry operations as Quaternions."""
if self.system == 'cubic': if self.system == 'cubic':
symQuats = [ symQuats = [
[ 1.0, 0.0, 0.0, 0.0 ], [ 1.0, 0.0, 0.0, 0.0 ],
@ -228,42 +153,40 @@ class Symmetry:
return np.array(symQuats) return np.array(symQuats)
def inFZ(self,rodrigues): def in_FZ(self,rho):
""" """
Check whether given Rodrigues-Frank vector falls into fundamental zone of own symmetry. Check whether given Rodrigues-Frank vector falls into fundamental zone.
Fundamental zone in Rodrigues space is point symmetric around origin. Fundamental zone in Rodrigues space is point symmetric around origin.
""" """
if (len(rodrigues) != 3): if(rho.shape[-1] != 3):
raise ValueError('Input is not a Rodrigues-Frank vector.\n') raise ValueError('Input is not a Rodrigues-Frank vector.')
if np.any(rodrigues == np.inf): return False # ToDo: MD: not sure if needed rho_abs = np.abs(rho)
Rabs = abs(rodrigues)
with np.errstate(invalid='ignore'):
# using '*'/prod for 'and'
if self.system == 'cubic': if self.system == 'cubic':
return np.sqrt(2.0)-1.0 >= Rabs[0] \ return np.where(np.prod(np.sqrt(2)-1. >= rho_abs,axis=-1) * \
and np.sqrt(2.0)-1.0 >= Rabs[1] \ (1. >= np.sum(rho_abs,axis=-1)),True,False)
and np.sqrt(2.0)-1.0 >= Rabs[2] \
and 1.0 >= Rabs[0] + Rabs[1] + Rabs[2]
elif self.system == 'hexagonal': elif self.system == 'hexagonal':
return 1.0 >= Rabs[0] and 1.0 >= Rabs[1] and 1.0 >= Rabs[2] \ return np.where(np.prod(1. >= rho_abs,axis=-1) * \
and 2.0 >= np.sqrt(3)*Rabs[0] + Rabs[1] \ (2. >= np.sqrt(3)*rho_abs[...,0] + rho_abs[...,1]) * \
and 2.0 >= np.sqrt(3)*Rabs[1] + Rabs[0] \ (2. >= np.sqrt(3)*rho_abs[...,1] + rho_abs[...,0]) * \
and 2.0 >= np.sqrt(3) + Rabs[2] (2. >= np.sqrt(3) + rho_abs[...,2]),True,False)
elif self.system == 'tetragonal': elif self.system == 'tetragonal':
return 1.0 >= Rabs[0] and 1.0 >= Rabs[1] \ return np.where(np.prod(1. >= rho_abs[...,:2],axis=-1) * \
and np.sqrt(2.0) >= Rabs[0] + Rabs[1] \ (np.sqrt(2) >= rho_abs[...,0] + rho_abs[...,1]) * \
and np.sqrt(2.0) >= Rabs[2] + 1.0 (np.sqrt(2) >= rho_abs[...,2] + 1.),True,False)
elif self.system == 'orthorhombic': elif self.system == 'orthorhombic':
return 1.0 >= Rabs[0] and 1.0 >= Rabs[1] and 1.0 >= Rabs[2] return np.where(np.prod(1. >= rho_abs,axis=-1),True,False)
else: else:
return True return np.where(np.all(np.isfinite(rho_abs),axis=-1),True,False)
def inDisorientationSST(self,rodrigues): def in_disorientation_SST(self,rho):
""" """
Check whether given Rodrigues-Frank vector (of misorientation) falls into standard stereographic triangle of own symmetry. Check whether given Rodrigues-Frank vector (of misorientation) falls into standard stereographic triangle.
References References
---------- ----------
@ -271,115 +194,32 @@ class Symmetry:
https://doi.org/10.1107/S0108767391006864 https://doi.org/10.1107/S0108767391006864
""" """
if (len(rodrigues) != 3): if(rho.shape[-1] != 3):
raise ValueError('Input is not a Rodrigues-Frank vector.\n') raise ValueError('Input is not a Rodrigues-Frank vector.')
R = rodrigues
epsilon = 0.0 with np.errstate(invalid='ignore'):
# using '*' for 'and'
if self.system == 'cubic': if self.system == 'cubic':
return R[0] >= R[1]+epsilon and R[1] >= R[2]+epsilon and R[2] >= epsilon return np.where((rho[...,0] >= rho[...,1]) * \
(rho[...,1] >= rho[...,2]) * \
(rho[...,2] >= 0),True,False)
elif self.system == 'hexagonal': elif self.system == 'hexagonal':
return R[0] >= np.sqrt(3)*(R[1]-epsilon) and R[1] >= epsilon and R[2] >= epsilon return np.where((rho[...,0] >= rho[...,1]*np.sqrt(3)) * \
(rho[...,1] >= 0) * \
(rho[...,2] >= 0),True,False)
elif self.system == 'tetragonal': elif self.system == 'tetragonal':
return R[0] >= R[1]-epsilon and R[1] >= epsilon and R[2] >= epsilon return np.where((rho[...,0] >= rho[...,1]) * \
(rho[...,1] >= 0) * \
(rho[...,2] >= 0),True,False)
elif self.system == 'orthorhombic': elif self.system == 'orthorhombic':
return R[0] >= epsilon and R[1] >= epsilon and R[2] >= epsilon return np.where((rho[...,0] >= 0) * \
(rho[...,1] >= 0) * \
(rho[...,2] >= 0),True,False)
else: else:
return True return np.ones_like(rho[...,0],dtype=bool)
def inSST(self, def in_SST(self,vector,proper=False,color=False):
vector,
proper = False,
color = False):
"""
Check whether given vector falls into standard stereographic triangle of own symmetry.
proper considers only vectors with z >= 0, hence uses two neighboring SSTs.
Return inverse pole figure color if requested.
Bases are computed from
>>> basis = {'cubic' : np.linalg.inv(np.array([[0.,0.,1.], # direction of red
... [1.,0.,1.]/np.sqrt(2.), # direction of green
... [1.,1.,1.]/np.sqrt(3.)]).T), # direction of blue
... 'hexagonal' : np.linalg.inv(np.array([[0.,0.,1.], # direction of red
... [1.,0.,0.], # direction of green
... [np.sqrt(3.),1.,0.]/np.sqrt(4.)]).T), # direction of blue
... 'tetragonal' : np.linalg.inv(np.array([[0.,0.,1.], # direction of red
... [1.,0.,0.], # direction of green
... [1.,1.,0.]/np.sqrt(2.)]).T), # direction of blue
... 'orthorhombic': np.linalg.inv(np.array([[0.,0.,1.], # direction of red
... [1.,0.,0.], # direction of green
... [0.,1.,0.]]).T), # direction of blue
... }
"""
if self.system == 'cubic':
basis = {'improper':np.array([ [-1. , 0. , 1. ],
[ np.sqrt(2.) , -np.sqrt(2.) , 0. ],
[ 0. , np.sqrt(3.) , 0. ] ]),
'proper':np.array([ [ 0. , -1. , 1. ],
[-np.sqrt(2.) , np.sqrt(2.) , 0. ],
[ np.sqrt(3.) , 0. , 0. ] ]),
}
elif self.system == 'hexagonal':
basis = {'improper':np.array([ [ 0. , 0. , 1. ],
[ 1. , -np.sqrt(3.) , 0. ],
[ 0. , 2. , 0. ] ]),
'proper':np.array([ [ 0. , 0. , 1. ],
[-1. , np.sqrt(3.) , 0. ],
[ np.sqrt(3.) , -1. , 0. ] ]),
}
elif self.system == 'tetragonal':
basis = {'improper':np.array([ [ 0. , 0. , 1. ],
[ 1. , -1. , 0. ],
[ 0. , np.sqrt(2.) , 0. ] ]),
'proper':np.array([ [ 0. , 0. , 1. ],
[-1. , 1. , 0. ],
[ np.sqrt(2.) , 0. , 0. ] ]),
}
elif self.system == 'orthorhombic':
basis = {'improper':np.array([ [ 0., 0., 1.],
[ 1., 0., 0.],
[ 0., 1., 0.] ]),
'proper':np.array([ [ 0., 0., 1.],
[-1., 0., 0.],
[ 0., 1., 0.] ]),
}
else: # direct exit for unspecified symmetry
if color:
return (True,np.zeros(3,'d'))
else:
return True
v = np.array(vector,dtype=float)
if proper: # check both improper ...
theComponents = np.around(np.dot(basis['improper'],v),12)
inSST = np.all(theComponents >= 0.0)
if not inSST: # ... and proper SST
theComponents = np.around(np.dot(basis['proper'],v),12)
inSST = np.all(theComponents >= 0.0)
else:
v[2] = abs(v[2]) # z component projects identical
theComponents = np.around(np.dot(basis['improper'],v),12) # for positive and negative values
inSST = np.all(theComponents >= 0.0)
if color: # have to return color array
if inSST:
rgb = np.power(theComponents/np.linalg.norm(theComponents),0.5) # smoothen color ramps
rgb = np.minimum(np.ones(3,dtype=float),rgb) # limit to maximum intensity
rgb /= max(rgb) # normalize to (HS)V = 1
else:
rgb = np.zeros(3,dtype=float)
return (inSST,rgb)
else:
return inSST
def in_SST(self,
vector,
proper = False,
color = False):
""" """
Check whether given vector falls into standard stereographic triangle of own symmetry. Check whether given vector falls into standard stereographic triangle of own symmetry.
@ -503,7 +343,7 @@ class Lattice: # ToDo: Make a subclass of Symmetry!
# transition to subclass # transition to subclass
self.system = self.symmetry.system self.system = self.symmetry.system
self.in_SST = self.symmetry.in_SST self.in_SST = self.symmetry.in_SST
self.inFZ = self.symmetry.inFZ self.in_FZ = self.symmetry.in_FZ
def __repr__(self): def __repr__(self):

View File

@ -40,8 +40,6 @@ class Orientation: # ToDo: make subclass of lattice and Rotation
self.rotation = Rotation.from_quaternion(rotation) # assume quaternion self.rotation = Rotation.from_quaternion(rotation) # assume quaternion
def __getitem__(self,item): def __getitem__(self,item):
if isinstance(item,tuple) and len(item) >= len(self):
raise IndexError('Too many indices')
return self.__class__(self.rotation[item],self.lattice) return self.__class__(self.rotation[item],self.lattice)
@ -61,8 +59,8 @@ class Orientation: # ToDo: make subclass of lattice and Rotation
if self.lattice.symmetry != other.lattice.symmetry: if self.lattice.symmetry != other.lattice.symmetry:
raise NotImplementedError('disorientation between different symmetry classes not supported yet.') raise NotImplementedError('disorientation between different symmetry classes not supported yet.')
mySymEqs = self.equivalentOrientations() if SST else self.equivalentOrientations([0]) # take all or only first sym operation mySymEqs = self.equivalent if SST else self.equivalent[0] # take all or only first sym operation
otherSymEqs = other.equivalentOrientations() otherSymEqs = other.equivalent
for i,sA in enumerate(mySymEqs): for i,sA in enumerate(mySymEqs):
aInv = sA.rotation.inversed() aInv = sA.rotation.inversed()
@ -71,7 +69,7 @@ class Orientation: # ToDo: make subclass of lattice and Rotation
r = b*aInv r = b*aInv
for k in range(2): for k in range(2):
r.inverse() r.inverse()
breaker = self.lattice.inFZ(r.as_Rodrigues(vector=True)) \ breaker = self.in_FZ \
and (not SST or other.lattice.symmetry.inDisorientationSST(r.as_Rodrigues(vector=True))) and (not SST or other.lattice.symmetry.inDisorientationSST(r.as_Rodrigues(vector=True)))
if breaker: break if breaker: break
if breaker: break if breaker: break
@ -81,17 +79,10 @@ class Orientation: # ToDo: make subclass of lattice and Rotation
# ... own sym, other sym, # ... own sym, other sym,
# self-->other: True, self<--other: False # self-->other: True, self<--other: False
def inFZ_vec(self):
def in_FZ(self):
"""Check if orientations fall into Fundamental Zone.""" """Check if orientations fall into Fundamental Zone."""
if not self.rotation.shape: return self.lattice.in_FZ(self.rotation.as_Rodrigues(vector=True))
return self.lattice.inFZ(self.rotation.as_Rodrigues(vector=True))
else:
return [self.lattice.inFZ(\
self.rotation.as_Rodrigues(vector=True)[l]) for l in range(self.rotation.shape[0])]
def inFZ(self):
return self.lattice.inFZ(self.rotation.as_Rodrigues(vector=True))
@property @property
def equivalent(self): def equivalent(self):
@ -112,16 +103,6 @@ class Orientation: # ToDo: make subclass of lattice and Rotation
return self.__class__(s@r,self.lattice) return self.__class__(s@r,self.lattice)
def equivalentOrientations(self,members=[]):
"""List of orientations which are symmetrically equivalent."""
try:
iter(members) # asking for (even empty) list of members?
except TypeError:
return self.__class__(self.lattice.symmetry.symmetryOperations(members)*self.rotation,self.lattice) # no, return rotation object
else:
return [self.__class__(q*self.rotation,self.lattice) \
for q in self.lattice.symmetry.symmetryOperations(members)] # yes, return list of rotations
def relatedOrientations_vec(self,model): def relatedOrientations_vec(self,model):
"""List of orientations related by the given orientation relationship.""" """List of orientations related by the given orientation relationship."""
h = self.lattice.relationOperations(model) h = self.lattice.relationOperations(model)
@ -149,7 +130,7 @@ class Orientation: # ToDo: make subclass of lattice and Rotation
r= 1 if not self.rotation.shape else equi.shape[1] #number of rotations r= 1 if not self.rotation.shape else equi.shape[1] #number of rotations
num_equi=equi.shape[0] #number of equivalente orientations num_equi=equi.shape[0] #number of equivalente orientations
quat= np.reshape( equi.as_quaternion(), (r*num_equi,4) ,order='F') #equivalents are listed in intiuitive order quat= np.reshape( equi.as_quaternion(), (r*num_equi,4) ,order='F') #equivalents are listed in intiuitive order
boolean=Orientation(quat, self.lattice).inFZ_vec() #check which ones are in FZ boolean=Orientation(quat, self.lattice).in_FZ() #check which ones are in FZ
if sum(boolean) == r: if sum(boolean) == r:
return self.__class__(quat[boolean],self.lattice) return self.__class__(quat[boolean],self.lattice)
@ -162,13 +143,10 @@ class Orientation: # ToDo: make subclass of lattice and Rotation
return self.__class__(quat[index],self.lattice) return self.__class__(quat[index],self.lattice)
def reduced(self): def reduced(self):
"""Transform orientation to fall into fundamental zone according to symmetry.""" """Transform orientation to fall into fundamental zone according to symmetry."""
for me in self.equivalentOrientations(): for me in self.equivalent:
if self.lattice.inFZ(me.rotation.as_Rodrigues(vector=True)): break if self.lattice.in_FZ(me.rotation.as_Rodrigues(vector=True)): break
return self.__class__(me.rotation,self.lattice) return self.__class__(me.rotation,self.lattice)
@ -179,35 +157,23 @@ class Orientation: # ToDo: make subclass of lattice and Rotation
SST = True): SST = True):
"""Axis rotated according to orientation (using crystal symmetry to ensure location falls into SST).""" """Axis rotated according to orientation (using crystal symmetry to ensure location falls into SST)."""
if SST: # pole requested to be within SST if SST: # pole requested to be within SST
for i,o in enumerate(self.equivalentOrientations()): # test all symmetric equivalent quaternions for i,o in enumerate(self.equivalent): # test all symmetric equivalent quaternions
pole = o.rotation@axis # align crystal direction to axis pole = o.rotation@axis # align crystal direction to axis
if self.lattice.symmetry.inSST(pole,proper): break # found SST version if self.lattice.in_SST(pole,proper): break # found SST version
else: else:
pole = self.rotation@axis # align crystal direction to axis pole = self.rotation@axis # align crystal direction to axis
return (pole,i if SST else 0) return (pole,i if SST else 0)
def IPFcolor(self,axis): def IPF_color(self,axis): #ToDo axis or direction?
"""TSL color of inverse pole figure for given axis.""" """TSL color of inverse pole figure for given axis."""
color = np.zeros(3,'d')
for o in self.equivalent:
pole = o.rotation@axis # align crystal direction to axis
inSST,color = self.lattice.symmetry.inSST(pole,color=True)
if inSST: break
return color
def IPF_color(self,axis):
"""TSL color of inverse pole figure for given axis. Not for hex or triclinic lattices."""
eq = self.equivalent eq = self.equivalent
pole = eq.rotation @ np.broadcast_to(axis/np.linalg.norm(axis),eq.rotation.shape+(3,)) pole = eq.rotation @ np.broadcast_to(axis/np.linalg.norm(axis),eq.rotation.shape+(3,))
in_SST, color = self.lattice.in_SST(pole,color=True) in_SST, color = self.lattice.in_SST(pole,color=True)
# ignore duplicates (occur for highly symmetric orientations) # ignore duplicates (occur for highly symmetric orientations)
found = np.zeros_like(in_SST[1],dtype=bool) found = np.zeros_like(in_SST[0],dtype=bool)
c = np.empty(color.shape[1:]) c = np.empty(color.shape[1:])
for s in range(in_SST.shape[0]): for s in range(in_SST.shape[0]):
c = np.where(np.expand_dims(np.logical_and(in_SST[s],~found),-1),color[s],c) c = np.where(np.expand_dims(np.logical_and(in_SST[s],~found),-1),color[s],c)
@ -220,17 +186,18 @@ class Orientation: # ToDo: make subclass of lattice and Rotation
def fromAverage(orientations, def fromAverage(orientations,
weights = []): weights = []):
"""Create orientation from average of list of orientations.""" """Create orientation from average of list of orientations."""
# further read: Orientation distribution analysis in deformed grains, https://doi.org/10.1107/S0021889801003077 # further read: Orientation distribution analysis in deformed grains
# https://doi.org/10.1107/S0021889801003077
if not all(isinstance(item, Orientation) for item in orientations): if not all(isinstance(item, Orientation) for item in orientations):
raise TypeError("Only instances of Orientation can be averaged.") raise TypeError("Only instances of Orientation can be averaged.")
closest = [] closest = []
ref = orientations[0] ref = orientations[0]
for o in orientations: for o in orientations:
closest.append(o.equivalentOrientations( closest.append(o.equivalent[
ref.disorientation(o, ref.disorientation(o,
SST = False, # select (o[ther]'s) sym orientation SST = False, # select (o[ther]'s) sym orientation
symmetries = True)[2]).rotation) # with lowest misorientation symmetries = True)[2]].rotation) # with lowest misorientation
return Orientation(Rotation.fromAverage(closest,weights),ref.lattice) return Orientation(Rotation.fromAverage(closest,weights),ref.lattice)

View File

@ -733,7 +733,7 @@ class Result:
@staticmethod @staticmethod
def _add_IPFcolor(q,l): def _add_IPF_color(q,l):
m = util.scale_to_coprime(np.array(l)) m = util.scale_to_coprime(np.array(l))
o = Orientation(Rotation(rfn.structured_to_unstructured(q['data'])), o = Orientation(Rotation(rfn.structured_to_unstructured(q['data'])),
@ -749,7 +749,7 @@ class Result:
'Creator': inspect.stack()[0][3][1:] 'Creator': inspect.stack()[0][3][1:]
} }
} }
def add_IPFcolor(self,q,l): def add_IPF_color(self,q,l):
""" """
Add RGB color tuple of inverse pole figure (IPF) color. Add RGB color tuple of inverse pole figure (IPF) color.
@ -761,7 +761,7 @@ class Result:
Lab frame direction for inverse pole figure. Lab frame direction for inverse pole figure.
""" """
self._add_generic_pointwise(self._add_IPFcolor,{'q':q},{'l':l}) self._add_generic_pointwise(self._add_IPF_color,{'q':q},{'l':l})
@staticmethod @staticmethod

View File

@ -470,7 +470,7 @@ class Rotation:
elif hasattr(shape, '__iter__'): elif hasattr(shape, '__iter__'):
r = np.random.random(tuple(shape)+(3,)) r = np.random.random(tuple(shape)+(3,))
else: else:
r = np.random.random((shape,3)) r = np.random.rand(shape,3)
A = np.sqrt(r[...,2]) A = np.sqrt(r[...,2])
B = np.sqrt(1.0-r[...,2]) B = np.sqrt(1.0-r[...,2])

View File

@ -3,10 +3,118 @@ import random
import pytest import pytest
import numpy as np import numpy as np
from damask import Orientation
from damask import Rotation
from damask import Symmetry from damask import Symmetry
def in_FZ(system,rho):
"""Non-vectorized version of 'in_FZ'."""
rho_abs = abs(rho)
if system == 'cubic':
return np.sqrt(2.0)-1.0 >= rho_abs[0] \
and np.sqrt(2.0)-1.0 >= rho_abs[1] \
and np.sqrt(2.0)-1.0 >= rho_abs[2] \
and 1.0 >= rho_abs[0] + rho_abs[1] + rho_abs[2]
elif system == 'hexagonal':
return 1.0 >= rho_abs[0] and 1.0 >= rho_abs[1] and 1.0 >= rho_abs[2] \
and 2.0 >= np.sqrt(3)*rho_abs[0] + rho_abs[1] \
and 2.0 >= np.sqrt(3)*rho_abs[1] + rho_abs[0] \
and 2.0 >= np.sqrt(3) + rho_abs[2]
elif system == 'tetragonal':
return 1.0 >= rho_abs[0] and 1.0 >= rho_abs[1] \
and np.sqrt(2.0) >= rho_abs[0] + rho_abs[1] \
and np.sqrt(2.0) >= rho_abs[2] + 1.0
elif system == 'orthorhombic':
return 1.0 >= rho_abs[0] and 1.0 >= rho_abs[1] and 1.0 >= rho_abs[2]
else:
return np.all(np.isfinite(rho_abs))
def in_disorientation_SST(system,rho):
"""Non-vectorized version of 'in_Disorientation_SST'."""
epsilon = 0.0
if system == 'cubic':
return rho[0] >= rho[1]+epsilon and rho[1] >= rho[2]+epsilon and rho[2] >= epsilon
elif system == 'hexagonal':
return rho[0] >= np.sqrt(3)*(rho[1]-epsilon) and rho[1] >= epsilon and rho[2] >= epsilon
elif system == 'tetragonal':
return rho[0] >= rho[1]-epsilon and rho[1] >= epsilon and rho[2] >= epsilon
elif system == 'orthorhombic':
return rho[0] >= epsilon and rho[1] >= epsilon and rho[2] >= epsilon
else:
return True
def in_SST(system,vector,proper = False):
"""Non-vectorized version of 'in_SST'."""
if system == 'cubic':
basis = {'improper':np.array([ [-1. , 0. , 1. ],
[ np.sqrt(2.) , -np.sqrt(2.) , 0. ],
[ 0. , np.sqrt(3.) , 0. ] ]),
'proper':np.array([ [ 0. , -1. , 1. ],
[-np.sqrt(2.) , np.sqrt(2.) , 0. ],
[ np.sqrt(3.) , 0. , 0. ] ]),
}
elif system == 'hexagonal':
basis = {'improper':np.array([ [ 0. , 0. , 1. ],
[ 1. , -np.sqrt(3.) , 0. ],
[ 0. , 2. , 0. ] ]),
'proper':np.array([ [ 0. , 0. , 1. ],
[-1. , np.sqrt(3.) , 0. ],
[ np.sqrt(3.) , -1. , 0. ] ]),
}
elif system == 'tetragonal':
basis = {'improper':np.array([ [ 0. , 0. , 1. ],
[ 1. , -1. , 0. ],
[ 0. , np.sqrt(2.) , 0. ] ]),
'proper':np.array([ [ 0. , 0. , 1. ],
[-1. , 1. , 0. ],
[ np.sqrt(2.) , 0. , 0. ] ]),
}
elif system == 'orthorhombic':
basis = {'improper':np.array([ [ 0., 0., 1.],
[ 1., 0., 0.],
[ 0., 1., 0.] ]),
'proper':np.array([ [ 0., 0., 1.],
[-1., 0., 0.],
[ 0., 1., 0.] ]),
}
else:
return True
v = np.array(vector,dtype=float)
if proper:
theComponents = np.around(np.dot(basis['improper'],v),12)
inSST = np.all(theComponents >= 0.0)
if not inSST:
theComponents = np.around(np.dot(basis['proper'],v),12)
inSST = np.all(theComponents >= 0.0)
else:
v[2] = abs(v[2])
theComponents = np.around(np.dot(basis['improper'],v),12)
inSST = np.all(theComponents >= 0.0)
return inSST
@pytest.fixture
def set_of_rodrigues(set_of_quaternions):
return Rotation(set_of_quaternions).as_Rodrigues(vector=True)[:200]
class TestSymmetry: class TestSymmetry:
@pytest.mark.parametrize('system',Symmetry.crystal_systems)
def test_in_FZ_vectorize(self,set_of_rodrigues,system):
for i,in_FZ_ in enumerate(Symmetry(system).in_FZ(set_of_rodrigues)):
assert in_FZ_ == in_FZ(system,set_of_rodrigues[i])
@pytest.mark.parametrize('system',Symmetry.crystal_systems)
def test_in_disorientation_SST_vectorize(self,set_of_rodrigues,system):
for i,in_disorientation_SST_ in enumerate(Symmetry(system).in_disorientation_SST(set_of_rodrigues)):
assert in_disorientation_SST_ == in_disorientation_SST(system,set_of_rodrigues[i])
@pytest.mark.parametrize('invalid_symmetry',['fcc','bcc','hello']) @pytest.mark.parametrize('invalid_symmetry',['fcc','bcc','hello'])
def test_invalid_symmetry(self,invalid_symmetry): def test_invalid_symmetry(self,invalid_symmetry):
with pytest.raises(KeyError): with pytest.raises(KeyError):
@ -22,19 +130,19 @@ class TestSymmetry:
assert Symmetry(symmetries[0]) != Symmetry(symmetries[1]) assert Symmetry(symmetries[0]) != Symmetry(symmetries[1])
@pytest.mark.parametrize('system',Symmetry.crystal_systems) @pytest.mark.parametrize('system',Symmetry.crystal_systems)
def test_inFZ(self,system): def test_in_FZ(self,system):
assert Symmetry(system).inFZ(np.zeros(3)) assert Symmetry(system).in_FZ(np.zeros(3))
@pytest.mark.parametrize('system',Symmetry.crystal_systems) @pytest.mark.parametrize('system',Symmetry.crystal_systems)
def test_inDisorientationSST(self,system): def test_in_disorientation_SST(self,system):
assert Symmetry(system).inDisorientationSST(np.zeros(3)) assert Symmetry(system).in_disorientation_SST(np.zeros(3))
@pytest.mark.parametrize('system',Symmetry.crystal_systems) @pytest.mark.parametrize('system',Symmetry.crystal_systems)
@pytest.mark.parametrize('proper',[True,False]) @pytest.mark.parametrize('proper',[True,False])
def test_inSST(self,system,proper): def test_in_SST(self,system,proper):
assert Symmetry(system).inSST(np.zeros(3),proper) assert Symmetry(system).in_SST(np.zeros(3),proper)
@pytest.mark.parametrize('function',['inFZ','inDisorientationSST']) @pytest.mark.parametrize('function',['in_FZ','in_disorientation_SST'])
def test_invalid_argument(self,function): def test_invalid_argument(self,function):
s = Symmetry() # noqa s = Symmetry() # noqa
with pytest.raises(ValueError): with pytest.raises(ValueError):

View File

@ -40,7 +40,7 @@ class TestOrientation:
@pytest.mark.parametrize('lattice',Lattice.lattices) @pytest.mark.parametrize('lattice',Lattice.lattices)
def test_IPF_equivalent(self,set_of_quaternions,lattice): def test_IPF_equivalent(self,set_of_quaternions,lattice):
direction = np.random.random(3)*2.0-1 direction = np.random.random(3)*2.0-1
for ori in Orientation(Rotation(set_of_quaternions),lattice)[200]: for ori in Orientation(Rotation(set_of_quaternions),lattice)[:200]:
color = ori.IPF_color(direction) color = ori.IPF_color(direction)
for equivalent in ori.equivalent: for equivalent in ori.equivalent:
assert np.allclose(color,equivalent.IPF_color(direction)) assert np.allclose(color,equivalent.IPF_color(direction))
@ -48,9 +48,10 @@ class TestOrientation:
@pytest.mark.parametrize('lattice',Lattice.lattices) @pytest.mark.parametrize('lattice',Lattice.lattices)
def test_IPF_vectorize(self,set_of_quaternions,lattice): def test_IPF_vectorize(self,set_of_quaternions,lattice):
for ori in Orientation(Rotation(set_of_quaternions),lattice)[200]:
direction = np.random.random(3)*2.0-1 direction = np.random.random(3)*2.0-1
assert np.allclose(ori.IPF_color(direction),IPF_color(ori,direction)) oris = Orientation(Rotation(set_of_quaternions),lattice)[:200]
for i,color in enumerate(oris.IPF_color(direction)):
assert np.allclose(color,IPF_color(oris[i],direction))
@pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch']) @pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch'])

View File

@ -153,8 +153,8 @@ class TestResult:
assert np.allclose(in_memory,in_file) assert np.allclose(in_memory,in_file)
@pytest.mark.parametrize('d',[[1,0,0],[0,1,0],[0,0,1]]) @pytest.mark.parametrize('d',[[1,0,0],[0,1,0],[0,0,1]])
def test_add_IPFcolor(self,default,d): def test_add_IPF_color(self,default,d):
default.add_IPFcolor('orientation',d) default.add_IPF_color('orientation',d)
loc = {'orientation': default.get_dataset_location('orientation'), loc = {'orientation': default.get_dataset_location('orientation'),
'color': default.get_dataset_location('IPFcolor_[{} {} {}]'.format(*d))} 'color': default.get_dataset_location('IPFcolor_[{} {} {}]'.format(*d))}
qu = default.read_dataset(loc['orientation']).view(np.double).reshape(-1,4) qu = default.read_dataset(loc['orientation']).view(np.double).reshape(-1,4)
@ -162,7 +162,7 @@ class TestResult:
in_memory = np.empty((qu.shape[0],3),np.uint8) in_memory = np.empty((qu.shape[0],3),np.uint8)
for i,q in enumerate(qu): for i,q in enumerate(qu):
o = damask.Orientation(q,crystal_structure).reduced() o = damask.Orientation(q,crystal_structure).reduced()
in_memory[i] = np.uint8(o.IPFcolor(np.array(d))*255) in_memory[i] = np.uint8(o.IPF_color(np.array(d))*255)
in_file = default.read_dataset(loc['color']) in_file = default.read_dataset(loc['color'])
assert np.allclose(in_memory,in_file) assert np.allclose(in_memory,in_file)

View File

@ -16,44 +16,6 @@ rot3= Rotation.from_random()
#average #average
class TestOrientation_vec: class TestOrientation_vec:
#@pytest.mark.xfail
@pytest.mark.parametrize('lattice',Lattice.lattices)
def test_equivalent_vec(self,lattice):
ori0=Orientation(rot0,lattice)
ori1=Orientation(rot1,lattice)
ori2=Orientation(rot2,lattice)
ori3=Orientation(rot3,lattice)
quat=np.array([rot0.as_quaternion(),rot1.as_quaternion(),rot2.as_quaternion(),rot3.as_quaternion()])
ori_vec=Orientation(quat,lattice)
for s in range(len(ori_vec.lattice.symmetry.symmetryOperations())):
assert all(ori_vec.equivalent.rotation.as_Eulers()[s,0] == \
ori0.equivalent[s].rotation.as_Eulers())
assert all(ori_vec.equivalent.rotation.as_quaternion()[s,1] == \
ori1.equivalent[s].rotation.as_quaternion())
assert all(ori_vec.equivalent.rotation.as_Rodrigues()[s,2] == \
ori2.equivalent[s].rotation.as_Rodrigues())
assert all(ori_vec.equivalent.rotation.as_cubochoric()[s,3] == \
ori3.equivalent[s].rotation.as_cubochoric())
@pytest.mark.parametrize('lattice',Lattice.lattices)
def test_inFZ_vec(self,lattice):
ori0=Orientation(rot0,lattice)
ori1=Orientation(rot1,lattice)
ori2=Orientation(rot2,lattice)
ori3=Orientation(rot3,lattice)
ori4=ori0.reduced() ; rot4=ori4.rotation #ensure 1 of them is in FZ
quat=np.array([rot0.as_quaternion(),rot1.as_quaternion(),\
rot2.as_quaternion(),rot3.as_quaternion(), rot4.as_quaternion()])
ori_vec=Orientation(quat,lattice)
assert ori_vec.inFZ_vec()[0] == ori0.inFZ()
assert ori_vec.inFZ_vec()[1] == ori1.inFZ()
assert ori_vec.inFZ_vec()[2] == ori2.inFZ()
assert ori_vec.inFZ_vec()[3] == ori3.inFZ()
assert ori_vec.inFZ_vec()[4] == ori4.inFZ()
@pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch']) @pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch'])