223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

Merge branch 'Vectorize-Orientation' into 'development' 

Vectorize orientation

See merge request damask/DAMASK!186
This commit is contained in:
Francisco Jose Gallardo Basile 2020-07-03 12:22:04 +02:00
parent 3526f421cc 208d5109d4
commit e1bbaac0d7
9 changed files with 641 additions and 353 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

261
python/damask/_lattice.py
View File

@ -1,11 +1,9 @@
import numpy as np
from . import Rotation
class Symmetry:
"""
Symmetry operations for lattice systems.
Symmetry-related operations for crystal systems.
References
----------
@ -13,34 +11,34 @@ class Symmetry:
"""
lattices = [None,'orthorhombic','tetragonal','hexagonal','cubic',]
crystal_systems = [None,'orthorhombic','tetragonal','hexagonal','cubic']
def __init__(self, symmetry = None):
def __init__(self, system = None):
"""
Symmetry Definition.
Parameters
----------
symmetry : str, optional
label of the crystal system
system : {None,'orthorhombic','tetragonal','hexagonal','cubic'}, optional
Name of the crystal system. Defaults to 'None'.
"""
if symmetry is not None and symmetry.lower() not in Symmetry.lattices:
raise KeyError(f'Symmetry/crystal system "{symmetry}" is unknown')
if system is not None and system.lower() not in self.crystal_systems:
raise KeyError(f'Crystal system "{system}" is unknown')
self.lattice = symmetry.lower() if isinstance(symmetry,str) else symmetry
self.system = system.lower() if isinstance(system,str) else system
def __copy__(self):
"""Copy."""
return self.__class__(self.lattice)
return self.__class__(self.system)
copy = __copy__
def __repr__(self):
"""Readable string."""
return f'{self.lattice}'
return f'{self.system}'
def __eq__(self, other):
@ -53,7 +51,7 @@ class Symmetry:
Symmetry to check for equality.
"""
return self.lattice == other.lattice
return self.system == other.system
def __neq__(self, other):
"""
@ -77,14 +75,16 @@ class Symmetry:
Symmetry to check for for order.
"""
myOrder = Symmetry.lattices.index(self.lattice)
otherOrder = Symmetry.lattices.index(other.lattice)
myOrder = self.crystal_systems.index(self.system)
otherOrder = self.crystal_systems.index(other.system)
return (myOrder > otherOrder) - (myOrder < otherOrder)
def symmetryOperations(self,members=[]):
"""List (or single element) of symmetry operations as rotations."""
if self.lattice == 'cubic':
symQuats = [
@property
def symmetry_operations(self):
"""Symmetry operations as quaternions."""
if self.system == 'cubic':
sym_quats = [
[ 1.0, 0.0, 0.0, 0.0 ],
[ 0.0, 1.0, 0.0, 0.0 ],
[ 0.0, 0.0, 1.0, 0.0 ],
@ -110,8 +110,8 @@ class Symmetry:
[-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.lattice == 'hexagonal':
symQuats = [
elif self.system == 'hexagonal':
sym_quats = [
[ 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) ],
@ -125,8 +125,8 @@ class Symmetry:
[ 0.0, -0.5, -0.5*np.sqrt(3), 0.0 ],
[ 0.0, 0.5*np.sqrt(3), 0.5, 0.0 ],
]
elif self.lattice == 'tetragonal':
symQuats = [
elif self.system == 'tetragonal':
sym_quats = [
[ 1.0, 0.0, 0.0, 0.0 ],
[ 0.0, 1.0, 0.0, 0.0 ],
[ 0.0, 0.0, 1.0, 0.0 ],
@ -136,64 +136,54 @@ class Symmetry:
[ 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.lattice == 'orthorhombic':
symQuats = [
elif self.system == 'orthorhombic':
sym_quats = [
[ 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 = [
sym_quats = [
[ 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
return np.array(sym_quats)
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.
"""
if (len(rodrigues) != 3):
raise ValueError('Input is not a Rodrigues-Frank vector.\n')
if(rho.shape[-1] != 3):
raise ValueError('Input is not a Rodrigues-Frank vector field.')
if np.any(rodrigues == np.inf): return False
rho_abs = np.abs(rho)
Rabs = abs(rodrigues)
if self.lattice == 'cubic':
return np.sqrt(2.0)-1.0 >= Rabs[0] \
and np.sqrt(2.0)-1.0 >= Rabs[1] \
and np.sqrt(2.0)-1.0 >= Rabs[2] \
and 1.0 >= Rabs[0] + Rabs[1] + Rabs[2]
elif self.lattice == 'hexagonal':
return 1.0 >= Rabs[0] and 1.0 >= Rabs[1] and 1.0 >= Rabs[2] \
and 2.0 >= np.sqrt(3)*Rabs[0] + Rabs[1] \
and 2.0 >= np.sqrt(3)*Rabs[1] + Rabs[0] \
and 2.0 >= np.sqrt(3) + Rabs[2]
elif self.lattice == 'tetragonal':
return 1.0 >= Rabs[0] and 1.0 >= Rabs[1] \
and np.sqrt(2.0) >= Rabs[0] + Rabs[1] \
and np.sqrt(2.0) >= Rabs[2] + 1.0
elif self.lattice == 'orthorhombic':
return 1.0 >= Rabs[0] and 1.0 >= Rabs[1] and 1.0 >= Rabs[2]
else:
return True
with np.errstate(invalid='ignore'):
# using '*'/prod for 'and'
if self.system == 'cubic':
return np.where(np.prod(np.sqrt(2)-1. >= rho_abs,axis=-1) * \
(1. >= np.sum(rho_abs,axis=-1)),True,False)
elif self.system == 'hexagonal':
return np.where(np.prod(1. >= rho_abs,axis=-1) * \
(2. >= np.sqrt(3)*rho_abs[...,0] + rho_abs[...,1]) * \
(2. >= np.sqrt(3)*rho_abs[...,1] + rho_abs[...,0]) * \
(2. >= np.sqrt(3) + rho_abs[...,2]),True,False)
elif self.system == 'tetragonal':
return np.where(np.prod(1. >= rho_abs[...,:2],axis=-1) * \
(np.sqrt(2) >= rho_abs[...,0] + rho_abs[...,1]) * \
(np.sqrt(2) >= rho_abs[...,2] + 1.),True,False)
elif self.system == 'orthorhombic':
return np.where(np.prod(1. >= rho_abs,axis=-1),True,False)
else:
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
----------
@ -201,27 +191,33 @@ class Symmetry:
https://doi.org/10.1107/S0108767391006864
"""
if (len(rodrigues) != 3):
raise ValueError('Input is not a Rodrigues-Frank vector.\n')
R = rodrigues
if(rho.shape[-1] != 3):
raise ValueError('Input is not a Rodrigues-Frank vector field.')
epsilon = 0.0
if self.lattice == 'cubic':
return R[0] >= R[1]+epsilon and R[1] >= R[2]+epsilon and R[2] >= epsilon
elif self.lattice == 'hexagonal':
return R[0] >= np.sqrt(3)*(R[1]-epsilon) and R[1] >= epsilon and R[2] >= epsilon
elif self.lattice == 'tetragonal':
return R[0] >= R[1]-epsilon and R[1] >= epsilon and R[2] >= epsilon
elif self.lattice == 'orthorhombic':
return R[0] >= epsilon and R[1] >= epsilon and R[2] >= epsilon
else:
return True
with np.errstate(invalid='ignore'):
# using '*' for 'and'
if self.system == 'cubic':
return np.where((rho[...,0] >= rho[...,1]) * \
(rho[...,1] >= rho[...,2]) * \
(rho[...,2] >= 0),True,False)
elif self.system == 'hexagonal':
return np.where((rho[...,0] >= rho[...,1]*np.sqrt(3)) * \
(rho[...,1] >= 0) * \
(rho[...,2] >= 0),True,False)
elif self.system == 'tetragonal':
return np.where((rho[...,0] >= rho[...,1]) * \
(rho[...,1] >= 0) * \
(rho[...,2] >= 0),True,False)
elif self.system == 'orthorhombic':
return np.where((rho[...,0] >= 0) * \
(rho[...,1] >= 0) * \
(rho[...,2] >= 0),True,False)
else:
return np.ones_like(rho[...,0],dtype=bool)
def inSST(self,
vector,
proper = False,
color = False):
#ToDo: IPF color in separate function
def in_SST(self,vector,proper=False,color=False):
"""
Check whether given vector falls into standard stereographic triangle of own symmetry.
@ -244,7 +240,10 @@ class Symmetry:
... }
"""
if self.lattice == 'cubic':
if(vector.shape[-1] != 3):
raise ValueError('Input is not a 3D vector field.')
if self.system == 'cubic':
basis = {'improper':np.array([ [-1. , 0. , 1. ],
[ np.sqrt(2.) , -np.sqrt(2.) , 0. ],
[ 0. , np.sqrt(3.) , 0. ] ]),
@ -252,7 +251,7 @@ class Symmetry:
[-np.sqrt(2.) , np.sqrt(2.) , 0. ],
[ np.sqrt(3.) , 0. , 0. ] ]),
}
elif self.lattice == 'hexagonal':
elif self.system == 'hexagonal':
basis = {'improper':np.array([ [ 0. , 0. , 1. ],
[ 1. , -np.sqrt(3.) , 0. ],
[ 0. , 2. , 0. ] ]),
@ -260,7 +259,7 @@ class Symmetry:
[-1. , np.sqrt(3.) , 0. ],
[ np.sqrt(3.) , -1. , 0. ] ]),
}
elif self.lattice == 'tetragonal':
elif self.system == 'tetragonal':
basis = {'improper':np.array([ [ 0. , 0. , 1. ],
[ 1. , -1. , 0. ],
[ 0. , np.sqrt(2.) , 0. ] ]),
@ -268,7 +267,7 @@ class Symmetry:
[-1. , 1. , 0. ],
[ np.sqrt(2.) , 0. , 0. ] ]),
}
elif self.lattice == 'orthorhombic':
elif self.system == 'orthorhombic':
basis = {'improper':np.array([ [ 0., 0., 1.],
[ 1., 0., 0.],
[ 0., 1., 0.] ]),
@ -278,43 +277,41 @@ class Symmetry:
}
else: # direct exit for unspecified symmetry
if color:
return (True,np.zeros(3,'d'))
return (np.ones_like(vector[...,0],bool),np.zeros_like(vector))
else:
return True
return np.ones_like(vector[...,0],bool)
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)
b_i = np.broadcast_to(basis['improper'],vector.shape+(3,))
if proper:
b_p = np.broadcast_to(basis['proper'], vector.shape+(3,))
improper = np.all(np.around(np.einsum('...ji,...i',b_i,vector),12)>=0.0,axis=-1,keepdims=True)
theComponents = np.where(np.broadcast_to(improper,vector.shape),
np.around(np.einsum('...ji,...i',b_i,vector),12),
np.around(np.einsum('...ji,...i',b_p,vector),12))
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)
vector_ = np.block([vector[...,0:2],np.abs(vector[...,2:3])]) # z component projects identical
theComponents = np.around(np.einsum('...ji,...i',b_i,vector_),12)
in_SST = np.all(theComponents >= 0.0,axis=-1)
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)
with np.errstate(invalid='ignore',divide='ignore'):
rgb = (theComponents/np.linalg.norm(theComponents,axis=-1,keepdims=True))**0.5 # smoothen color ramps
rgb = np.minimum(1.,rgb) # limit to maximum intensity
rgb /= np.max(rgb,axis=-1,keepdims=True) # normalize to (HS)V = 1
rgb[np.broadcast_to(~in_SST.reshape(vector[...,0].shape+(1,)),vector.shape)] = 0.0
return (in_SST,rgb)
else:
return inSST
# code derived from https://github.com/ezag/pyeuclid
# suggested reading: http://web.mit.edu/2.998/www/QuaternionReport1.pdf
return in_SST
# ******************************************************************************************
class Lattice:
class Lattice: # ToDo: Make a subclass of Symmetry!
"""
Lattice system.
Bravais lattice.
Currently, this contains only a mapping from Bravais lattice to symmetry
This contains only a mapping from Bravais lattice to symmetry
and orientation relationships. It could include twin and slip systems.
References
@ -324,15 +321,15 @@ class Lattice:
"""
lattices = {
'triclinic':{'symmetry':None},
'bct':{'symmetry':'tetragonal'},
'hex':{'symmetry':'hexagonal'},
'fcc':{'symmetry':'cubic','c/a':1.0},
'bcc':{'symmetry':'cubic','c/a':1.0},
'triclinic':{'system':None},
'bct': {'system':'tetragonal'},
'hex': {'system':'hexagonal'},
'fcc': {'system':'cubic','c/a':1.0},
'bcc': {'system':'cubic','c/a':1.0},
}
def __init__(self, lattice):
def __init__(self,lattice,c_over_a=None):
"""
New lattice of given type.
@ -343,18 +340,23 @@ class Lattice:
"""
self.lattice = lattice
self.symmetry = Symmetry(self.lattices[lattice]['symmetry'])
self.symmetry = Symmetry(self.lattices[lattice]['system'])
# transition to subclass
self.system = self.symmetry.system
self.in_SST = self.symmetry.in_SST
self.in_FZ = self.symmetry.in_FZ
self.in_disorientation_SST = self.symmetry.in_disorientation_SST
def __repr__(self):
"""Report basic lattice information."""
return f'Bravais lattice {self.lattice} ({self.symmetry} symmetry)'
return f'Bravais lattice {self.lattice} ({self.symmetry} crystal system)'
# Kurdjomov--Sachs orientation relationship for fcc <-> bcc transformation
# from S. Morito et al., Journal of Alloys and Compounds 577:s587-s592, 2013
# also see K. Kitahara et al., Acta Materialia 54:1279-1288, 2006
KS = {'mapping':{'fcc':0,'bcc':1},
_KS = {'mapping':{'fcc':0,'bcc':1},
'planes': np.array([
[[ 1, 1, 1],[ 0, 1, 1]],
[[ 1, 1, 1],[ 0, 1, 1]],
@ -408,7 +410,7 @@ class Lattice:
# Greninger--Troiano orientation relationship for fcc <-> bcc transformation
# from Y. He et al., Journal of Applied Crystallography 39:72-81, 2006
GT = {'mapping':{'fcc':0,'bcc':1},
_GT = {'mapping':{'fcc':0,'bcc':1},
'planes': np.array([
[[ 1, 1, 1],[ 1, 0, 1]],
[[ 1, 1, 1],[ 1, 1, 0]],
@ -462,7 +464,7 @@ class Lattice:
# Greninger--Troiano' orientation relationship for fcc <-> bcc transformation
# from Y. He et al., Journal of Applied Crystallography 39:72-81, 2006
GTprime = {'mapping':{'fcc':0,'bcc':1},
_GTprime = {'mapping':{'fcc':0,'bcc':1},
'planes': np.array([
[[ 7, 17, 17],[ 12, 5, 17]],
[[ 17, 7, 17],[ 17, 12, 5]],
@ -516,7 +518,7 @@ class Lattice:
# Nishiyama--Wassermann orientation relationship for fcc <-> bcc transformation
# from H. Kitahara et al., Materials Characterization 54:378-386, 2005
NW = {'mapping':{'fcc':0,'bcc':1},
_NW = {'mapping':{'fcc':0,'bcc':1},
'planes': np.array([
[[ 1, 1, 1],[ 0, 1, 1]],
[[ 1, 1, 1],[ 0, 1, 1]],
@ -546,7 +548,7 @@ class Lattice:
# Pitsch orientation relationship for fcc <-> bcc transformation
# from Y. He et al., Acta Materialia 53:1179-1190, 2005
Pitsch = {'mapping':{'fcc':0,'bcc':1},
_Pitsch = {'mapping':{'fcc':0,'bcc':1},
'planes': np.array([
[[ 0, 1, 0],[ -1, 0, 1]],
[[ 0, 0, 1],[ 1, -1, 0]],
@ -576,7 +578,7 @@ class Lattice:
# Bain orientation relationship for fcc <-> bcc transformation
# from Y. He et al., Journal of Applied Crystallography 39:72-81, 2006
Bain = {'mapping':{'fcc':0,'bcc':1},
_Bain = {'mapping':{'fcc':0,'bcc':1},
'planes': np.array([
[[ 1, 0, 0],[ 1, 0, 0]],
[[ 0, 1, 0],[ 0, 1, 0]],
@ -586,7 +588,8 @@ class Lattice:
[[ 0, 0, 1],[ 1, 0, 1]],
[[ 1, 0, 0],[ 1, 1, 0]]],dtype='float')}
def relationOperations(self,model):
def relation_operations(self,model):
"""
Crystallographic orientation relationships for phase transformations.
@ -608,8 +611,8 @@ class Lattice:
https://doi.org/10.1016/j.actamat.2004.11.021
"""
models={'KS':self.KS, 'GT':self.GT, 'GT_prime':self.GTprime,
'NW':self.NW, 'Pitsch': self.Pitsch, 'Bain':self.Bain}
models={'KS':self._KS, 'GT':self._GT, 'GT_prime':self._GTprime,
'NW':self._NW, 'Pitsch': self._Pitsch, 'Bain':self._Bain}
try:
relationship = models[model]
except KeyError :
@ -635,6 +638,8 @@ class Lattice:
otherDir = miller[otherDir_id]/ np.linalg.norm(miller[otherDir_id])
otherMatrix = np.array([otherDir,np.cross(otherPlane,otherDir),otherPlane])
r['rotations'].append(Rotation.from_matrix(np.dot(otherMatrix.T,myMatrix)))
r['rotations'].append(np.dot(otherMatrix.T,myMatrix))
r['rotations'] = np.array(r['rotations'])
return r

143
python/damask/_orientation.py
View File

@ -3,7 +3,7 @@ import numpy as np
from . import Lattice
from . import Rotation
class Orientation:
class Orientation: # ToDo: make subclass of lattice and Rotation?
"""
Crystallographic orientation.
@ -39,9 +39,12 @@ class Orientation:
else:
self.rotation = Rotation.from_quaternion(rotation) # assume quaternion
if self.rotation.quaternion.shape != (4,):
raise NotImplementedError('Support for multiple rotations missing')
def __getitem__(self,item):
"""Iterate over leading/leftmost dimension of Orientation array."""
return self.__class__(self.rotation[item],self.lattice)
# ToDo: Discuss vectorization/calling signature
def disorientation(self,
other,
SST = True,
@ -58,8 +61,8 @@ class Orientation:
if self.lattice.symmetry != other.lattice.symmetry:
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
otherSymEqs = other.equivalentOrientations()
mySymEqs = self.equivalent if SST else self.equivalent[0] #ToDo: This is just me! # take all or only first sym operation
otherSymEqs = other.equivalent
for i,sA in enumerate(mySymEqs):
aInv = sA.rotation.inversed()
@ -68,8 +71,8 @@ class Orientation:
r = b*aInv
for k in range(2):
r.inverse()
breaker = self.lattice.symmetry.inFZ(r.as_Rodrigues(vector=True)) \
and (not SST or other.lattice.symmetry.inDisorientationSST(r.as_Rodrigues(vector=True)))
breaker = self.lattice.in_FZ(r.as_Rodrigues(vector=True)) \
and (not SST or other.lattice.in_disorientation_SST(r.as_Rodrigues(vector=True)))
if breaker: break
if breaker: break
if breaker: break
@ -77,79 +80,123 @@ class Orientation:
return (Orientation(r,self.lattice), i,j, k == 1) if symmetries else r # disorientation ...
# ... own sym, other sym,
# self-->other: True, self<--other: False
def inFZ(self):
return self.lattice.symmetry.inFZ(self.rotation.as_Rodrigues(vector=True))
@property
def in_FZ(self):
"""Check if orientations fall into Fundamental Zone."""
return self.lattice.in_FZ(self.rotation.as_Rodrigues(vector=True))
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
@property
def equivalent(self):
"""
Orientations which are symmetrically equivalent.
def relatedOrientations(self,model):
"""List of orientations related by the given orientation relationship."""
r = self.lattice.relationOperations(model)
return [self.__class__(o*self.rotation,r['lattice']) for o in r['rotations']]
One dimension (length according to number of symmetrically equivalent orientations)
is added to the left of the Rotation array.
"""
o = self.lattice.symmetry.symmetry_operations
o = o.reshape(o.shape[:1]+(1,)*len(self.rotation.shape)+(4,))
o = Rotation(np.broadcast_to(o,o.shape[:1]+self.rotation.quaternion.shape))
s = np.broadcast_to(self.rotation.quaternion,o.shape[:1]+self.rotation.quaternion.shape)
return self.__class__(o@Rotation(s),self.lattice)
def related(self,model):
"""
Orientations related by the given orientation relationship.
One dimension (length according to number of related orientations)
is added to the left of the Rotation array.
"""
o = Rotation.from_matrix(self.lattice.relation_operations(model)['rotations']).as_quaternion()
o = o.reshape(o.shape[:1]+(1,)*len(self.rotation.shape)+(4,))
o = Rotation(np.broadcast_to(o,o.shape[:1]+self.rotation.quaternion.shape))
s = np.broadcast_to(self.rotation.quaternion,o.shape[:1]+self.rotation.quaternion.shape)
return self.__class__(o@Rotation(s),self.lattice.relation_operations(model)['lattice'])
@property
def reduced(self):
"""Transform orientation to fall into fundamental zone according to symmetry."""
for me in self.equivalentOrientations():
if self.lattice.symmetry.inFZ(me.rotation.as_Rodrigues(vector=True)): break
eq = self.equivalent
in_FZ = eq.in_FZ
return self.__class__(me.rotation,self.lattice)
# remove duplicates (occur for highly symmetric orientations)
found = np.zeros_like(in_FZ[0],dtype=bool)
q = self.rotation.quaternion[0]
for s in range(in_FZ.shape[0]):
#something fishy... why does q needs to be initialized?
q = np.where(np.expand_dims(np.logical_and(in_FZ[s],~found),-1),eq.rotation.quaternion[s],q)
found = np.logical_or(in_FZ[s],found)
return self.__class__(q,self.lattice)
def inversePole(self,
axis,
proper = False,
SST = True):
def inverse_pole(self,axis,proper=False,SST=True):
"""Axis rotated according to orientation (using crystal symmetry to ensure location falls into SST)."""
if SST: # pole requested to be within SST
for i,o in enumerate(self.equivalentOrientations()): # test all symmetric equivalent quaternions
pole = o.rotation*axis # align crystal direction to axis
if self.lattice.symmetry.inSST(pole,proper): break # found SST version
if SST:
eq = self.equivalent
pole = eq.rotation @ np.broadcast_to(axis/np.linalg.norm(axis),eq.rotation.shape+(3,))
in_SST = self.lattice.in_SST(pole,proper=proper)
# remove duplicates (occur for highly symmetric orientations)
found = np.zeros_like(in_SST[0],dtype=bool)
p = pole[0]
for s in range(in_SST.shape[0]):
p = np.where(np.expand_dims(np.logical_and(in_SST[s],~found),-1),pole[s],p)
found = np.logical_or(in_SST[s],found)
return p
else:
pole = self.rotation*axis # align crystal direction to axis
return (pole,i if SST else 0)
return self.rotation @ np.broadcast_to(axis/np.linalg.norm(axis),self.rotation.shape+(3,))
def IPFcolor(self,axis):
def IPF_color(self,axis): #ToDo axis or direction?
"""TSL color of inverse pole figure for given axis."""
color = np.zeros(3,'d')
eq = self.equivalent
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)
for o in self.equivalentOrientations():
pole = o.rotation*axis # align crystal direction to axis
inSST,color = self.lattice.symmetry.inSST(pole,color=True)
if inSST: break
# remove duplicates (occur for highly symmetric orientations)
found = np.zeros_like(in_SST[0],dtype=bool)
c = color[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)
found = np.logical_or(in_SST[s],found)
return color
return c
# ToDo: Discuss vectorization/calling signature
@staticmethod
def fromAverage(orientations,
def from_average(orientations,
weights = []):
"""Create orientation from average of list of orientations."""
# further read: Orientation distribution analysis in deformed grains
# https://doi.org/10.1107/S0021889801003077
if not all(isinstance(item, Orientation) for item in orientations):
raise TypeError("Only instances of Orientation can be averaged.")
closest = []
ref = orientations[0]
for o in orientations:
closest.append(o.equivalentOrientations(
closest.append(o.equivalent[
ref.disorientation(o,
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.from_average(closest,weights),ref.lattice)
# ToDo: Discuss vectorization/calling signature
def average(self,other):
"""Calculate the average rotation."""
return Orientation.fromAverage([self,other])
return Orientation.from_average([self,other])

25
python/damask/_result.py
View File

@ -11,6 +11,7 @@ from functools import partial
import h5py
import numpy as np
from numpy.lib import recfunctions as rfn
import damask
from . import VTK
@ -731,29 +732,23 @@ class Result:
@staticmethod
def _add_IPFcolor(q,l):
d = np.array(l)
d_unit = d/np.linalg.norm(d)
m = util.scale_to_coprime(d)
colors = np.empty((len(q['data']),3),np.uint8)
def _add_IPF_color(q,l):
m = util.scale_to_coprime(np.array(l))
lattice = q['meta']['Lattice']
for i,qu in enumerate(q['data']):
o = Orientation(np.array([qu['w'],qu['x'],qu['y'],qu['z']]),lattice).reduced()
colors[i] = np.uint8(o.IPFcolor(d_unit)*255)
o = Orientation(Rotation(rfn.structured_to_unstructured(q['data'])),
lattice = q['meta']['Lattice'])
return {
'data': colors,
'data': np.uint8(o.IPF_color(l)*255),
'label': 'IPFcolor_[{} {} {}]'.format(*m),
'meta' : {
'Unit': 'RGB (8bit)',
'Lattice': lattice,
'Unit': '8-bit RGB',
'Lattice': q['meta']['Lattice'],
'Description': 'Inverse Pole Figure (IPF) colors along sample direction [{} {} {}]'.format(*m),
'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.
@ -765,7 +760,7 @@ class Result:
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

44
python/damask/_rotation.py
View File

@ -53,6 +53,8 @@ class Rotation:
Use .from_quaternion to perform a sanity check.
"""
if quaternion.shape[-1] != 4:
raise ValueError('Not a quaternion')
self.quaternion = quaternion.copy()
@ -61,6 +63,7 @@ class Rotation:
return self.quaternion.shape[:-1]
# ToDo: Check difference __copy__ vs __deepcopy__
def __copy__(self):
"""Copy."""
return self.__class__(self.quaternion)
@ -71,7 +74,7 @@ class Rotation:
def __repr__(self):
"""Orientation displayed as unit quaternion, rotation matrix, and Bunge-Euler angles."""
if self.quaternion.shape != (4,):
raise NotImplementedError('Support for multiple rotations missing')
return 'Quaternions:\n'+str(self.quaternion) # ToDo: could be nicer ...
return '\n'.join([
'Quaternion: (real={:.3f}, imag=<{:+.3f}, {:+.3f}, {:+.3f}>)'.format(*(self.quaternion)),
'Matrix:\n{}'.format(np.round(self.as_matrix(),8)),
@ -79,6 +82,19 @@ class Rotation:
])
def __getitem__(self,item):
"""Iterate over leading/leftmost dimension of Rotation array."""
if self.shape == (): return self.copy()
if isinstance(item,tuple) and len(item) >= len(self):
raise IndexError('Too many indices')
return self.__class__(self.quaternion[item])
def __len__(self):
"""Length of leading/leftmost dimension of Rotation array."""
return 0 if self.shape == () else self.shape[0]
def __matmul__(self, other):
"""
Rotation of vector, second or fourth order tensor, or rotation object.
@ -88,6 +104,11 @@ class Rotation:
other : numpy.ndarray or Rotation
Vector, second or fourth order tensor, or rotation object that is rotated.
Returns
-------
other_rot : numpy.ndarray or Rotation
Rotated vector, second or fourth order tensor, or rotation object.
"""
if isinstance(other, Rotation):
q_m = self.quaternion[...,0:1]
@ -178,7 +199,7 @@ class Rotation:
"""
if self.quaternion.shape != (4,) or other.quaternion.shape != (4,):
raise NotImplementedError('Support for multiple rotations missing')
return Rotation.fromAverage([self,other])
return Rotation.from_average([self,other])
################################################################################################
@ -273,7 +294,11 @@ class Rotation:
"""
ro = Rotation._qu2ro(self.quaternion)
return ro[...,:3]*ro[...,3] if vector else ro
if vector:
with np.errstate(invalid='ignore'):
return ro[...,:3]*ro[...,3:4]
else:
return ro
def as_homochoric(self):
"""
@ -299,6 +324,7 @@ class Rotation:
"""
return Rotation._qu2cu(self.quaternion)
@property
def M(self): # ToDo not sure about the name: as_M or M? we do not have a from_M
"""
Intermediate representation supporting quaternion averaging.
@ -555,7 +581,7 @@ class Rotation:
@staticmethod
def fromAverage(rotations,weights = None):
def from_average(rotations,weights = None):
"""
Average rotation.
@ -580,8 +606,8 @@ class Rotation:
weights = np.ones(N,dtype='i')
for i,(r,n) in enumerate(zip(rotations,weights)):
M = r.M() * n if i == 0 \
else M + r.M() * n # noqa add (multiples) of this rotation to average noqa
M = r.M * n if i == 0 \
else M + r.M * n # noqa add (multiples) of this rotation to average noqa
eig, vec = np.linalg.eig(M/N)
return Rotation.from_quaternion(np.real(vec.T[eig.argmax()]),accept_homomorph = True)
@ -593,7 +619,7 @@ class Rotation:
elif hasattr(shape, '__iter__'):
r = np.random.random(tuple(shape)+(3,))
else:
r = np.random.random((shape,3))
r = np.random.rand(shape,3)
A = np.sqrt(r[...,2])
B = np.sqrt(1.0-r[...,2])
@ -604,11 +630,9 @@ class Rotation:
return Rotation(q.reshape(r.shape[:-1]+(4,)) if shape is not None else q)._standardize()
# for compatibility (old names do not follow convention)
asM = M
fromQuaternion = from_quaternion
fromEulers = from_Eulers
fromQuaternion = from_quaternion
asAxisAngle = as_axis_angle
__mul__ = __matmul__

90
python/tests/conftest.py
View File

@ -1,4 +1,5 @@
from pathlib import Path
import numpy as np
import pytest
@ -21,3 +22,92 @@ def update(request):
def reference_dir_base():
"""Directory containing reference results."""
return Path(__file__).parent/'reference'
@pytest.fixture
def set_of_quaternions():
"""A set of n random rotations."""
def random_quaternions(N):
r = np.random.rand(N,3)
A = np.sqrt(r[:,2])
B = np.sqrt(1.0-r[:,2])
qu = np.column_stack([np.cos(2.0*np.pi*r[:,0])*A,
np.sin(2.0*np.pi*r[:,1])*B,
np.cos(2.0*np.pi*r[:,1])*B,
np.sin(2.0*np.pi*r[:,0])*A])
qu[:,0]*=np.sign(qu[:,0])
return qu
n = 1100
scatter=1.e-2
specials = np.array([
[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,-1.0, 0.0, 0.0],
[0.0, 0.0,-1.0, 0.0],
[0.0, 0.0, 0.0,-1.0],
#----------------------
[1.0, 1.0, 0.0, 0.0],
[1.0, 0.0, 1.0, 0.0],
[1.0, 0.0, 0.0, 1.0],
[0.0, 1.0, 1.0, 0.0],
[0.0, 1.0, 0.0, 1.0],
[0.0, 0.0, 1.0, 1.0],
#----------------------
[1.0,-1.0, 0.0, 0.0],
[1.0, 0.0,-1.0, 0.0],
[1.0, 0.0, 0.0,-1.0],
[0.0, 1.0,-1.0, 0.0],
[0.0, 1.0, 0.0,-1.0],
[0.0, 0.0, 1.0,-1.0],
#----------------------
[0.0, 1.0,-1.0, 0.0],
[0.0, 1.0, 0.0,-1.0],
[0.0, 0.0, 1.0,-1.0],
#----------------------
[0.0,-1.0,-1.0, 0.0],
[0.0,-1.0, 0.0,-1.0],
[0.0, 0.0,-1.0,-1.0],
#----------------------
[1.0, 1.0, 1.0, 0.0],
[1.0, 1.0, 0.0, 1.0],
[1.0, 0.0, 1.0, 1.0],
[1.0,-1.0, 1.0, 0.0],
[1.0,-1.0, 0.0, 1.0],
[1.0, 0.0,-1.0, 1.0],
[1.0, 1.0,-1.0, 0.0],
[1.0, 1.0, 0.0,-1.0],
[1.0, 0.0, 1.0,-1.0],
[1.0,-1.0,-1.0, 0.0],
[1.0,-1.0, 0.0,-1.0],
[1.0, 0.0,-1.0,-1.0],
#----------------------
[0.0, 1.0, 1.0, 1.0],
[0.0, 1.0,-1.0, 1.0],
[0.0, 1.0, 1.0,-1.0],
[0.0,-1.0, 1.0, 1.0],
[0.0,-1.0,-1.0, 1.0],
[0.0,-1.0, 1.0,-1.0],
[0.0,-1.0,-1.0,-1.0],
#----------------------
[1.0, 1.0, 1.0, 1.0],
[1.0,-1.0, 1.0, 1.0],
[1.0, 1.0,-1.0, 1.0],
[1.0, 1.0, 1.0,-1.0],
[1.0,-1.0,-1.0, 1.0],
[1.0,-1.0, 1.0,-1.0],
[1.0, 1.0,-1.0,-1.0],
[1.0,-1.0,-1.0,-1.0],
])
specials /= np.linalg.norm(specials,axis=1).reshape(-1,1)
specials_scatter = specials + np.broadcast_to(np.random.rand(4)*scatter,specials.shape)
specials_scatter /= np.linalg.norm(specials_scatter,axis=1).reshape(-1,1)
specials_scatter[specials_scatter[:,0]<0]*=-1
return np.array([s for s in specials] + \
[s for s in specials_scatter] + \
[s for s in random_quaternions(n-2*len(specials))])

140
python/tests/test_Lattice.py
View File

@ -3,38 +3,154 @@ import random
import pytest
import numpy as np
from damask import Rotation
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:
@pytest.mark.parametrize('system',Symmetry.crystal_systems)
def test_in_FZ_vectorize(self,set_of_rodrigues,system):
result = Symmetry(system).in_FZ(set_of_rodrigues.reshape(50,4,3)).reshape(200)
for i,r in enumerate(result):
assert r == 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):
result = Symmetry(system).in_disorientation_SST(set_of_rodrigues.reshape(50,4,3)).reshape(200)
for i,r in enumerate(result):
assert r == in_disorientation_SST(system,set_of_rodrigues[i])
@pytest.mark.parametrize('proper',[True,False])
@pytest.mark.parametrize('system',Symmetry.crystal_systems)
def test_in_SST_vectorize(self,system,proper):
vecs = np.random.rand(20,4,3)
result = Symmetry(system).in_SST(vecs,proper).reshape(20*4)
for i,r in enumerate(result):
assert r == in_SST(system,vecs.reshape(20*4,3)[i],proper)
@pytest.mark.parametrize('invalid_symmetry',['fcc','bcc','hello'])
def test_invalid_symmetry(self,invalid_symmetry):
with pytest.raises(KeyError):
s = Symmetry(invalid_symmetry) # noqa
def test_equal(self):
symmetry = random.choice(Symmetry.lattices)
symmetry = random.choice(Symmetry.crystal_systems)
print(symmetry)
assert Symmetry(symmetry) == Symmetry(symmetry)
def test_not_equal(self):
symmetries = random.sample(Symmetry.lattices,k=2)
symmetries = random.sample(Symmetry.crystal_systems,k=2)
assert Symmetry(symmetries[0]) != Symmetry(symmetries[1])
@pytest.mark.parametrize('lattice',Symmetry.lattices)
def test_inFZ(self,lattice):
assert Symmetry(lattice).inFZ(np.zeros(3))
@pytest.mark.parametrize('system',Symmetry.crystal_systems)
def test_in_FZ(self,system):
assert Symmetry(system).in_FZ(np.zeros(3))
@pytest.mark.parametrize('lattice',Symmetry.lattices)
def test_inDisorientationSST(self,lattice):
assert Symmetry(lattice).inDisorientationSST(np.zeros(3))
@pytest.mark.parametrize('system',Symmetry.crystal_systems)
def test_in_disorientation_SST(self,system):
assert Symmetry(system).in_disorientation_SST(np.zeros(3))
@pytest.mark.parametrize('lattice',Symmetry.lattices)
@pytest.mark.parametrize('system',Symmetry.crystal_systems)
@pytest.mark.parametrize('proper',[True,False])
def test_inSST(self,lattice,proper):
assert Symmetry(lattice).inSST(np.zeros(3),proper)
def test_in_SST(self,system,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','in_SST'])
def test_invalid_argument(self,function):
s = Symmetry() # noqa
with pytest.raises(ValueError):

94
python/tests/test_Orientation.py
View File

@ -11,6 +11,25 @@ from damask import Lattice
n = 1000
def IPF_color(orientation,direction):
"""TSL color of inverse pole figure for given axis (non-vectorized)."""
for o in orientation.equivalent:
pole = o.rotation@direction
inSST,color = orientation.lattice.in_SST(pole,color=True)
if inSST: break
return color
def inverse_pole(orientation,axis,proper=False,SST=True):
if SST:
for eq in orientation.equivalent:
pole = eq.rotation @ axis/np.linalg.norm(axis)
if orientation.lattice.in_SST(pole,proper=proper):
return pole
else:
return orientation.rotation @ axis/np.linalg.norm(axis)
@pytest.fixture
def reference_dir(reference_dir_base):
"""Directory containing reference results."""
@ -19,6 +38,31 @@ def reference_dir(reference_dir_base):
class TestOrientation:
@pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch'])
@pytest.mark.parametrize('lattice',['fcc','bcc'])
def test_relationship_vectorize(self,set_of_quaternions,lattice,model):
result = Orientation(set_of_quaternions[:200].reshape(50,4,4),lattice).related(model)
ref_qu = result.rotation.quaternion.reshape(-1,200,4)
for i in range(200):
single = Orientation(set_of_quaternions[i],lattice).related(model).rotation.quaternion
assert np.allclose(ref_qu[:,i,:],single)
@pytest.mark.parametrize('lattice',Lattice.lattices)
def test_IPF_vectorize(self,set_of_quaternions,lattice):
direction = np.random.random(3)*2.0-1
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('SST',[False,True])
@pytest.mark.parametrize('proper',[True,False])
@pytest.mark.parametrize('lattice',Lattice.lattices)
def test_inverse_pole_vectorize(self,set_of_quaternions,lattice,SST,proper):
axis = np.random.random(3)*2.0-1
oris = Orientation(Rotation(set_of_quaternions),lattice)[:200]
for i,pole in enumerate(oris.inverse_pole(axis,SST=SST)):
assert np.allclose(pole,inverse_pole(oris[i],axis,SST=SST))
@pytest.mark.parametrize('color',[{'label':'red', 'RGB':[1,0,0],'direction':[0,0,1]},
{'label':'green','RGB':[0,1,0],'direction':[0,1,1]},
{'label':'blue', 'RGB':[0,0,1],'direction':[1,1,1]}])
@ -26,35 +70,63 @@ class TestOrientation:
def test_IPF_cubic(self,color,lattice):
cube = damask.Orientation(damask.Rotation(),lattice)
for direction in set(permutations(np.array(color['direction']))):
assert np.allclose(cube.IPFcolor(np.array(direction)),np.array(color['RGB']))
assert np.allclose(cube.IPF_color(np.array(direction)),np.array(color['RGB']))
@pytest.mark.parametrize('lattice',Lattice.lattices)
def test_IPF(self,lattice):
def test_IPF_equivalent(self,set_of_quaternions,lattice):
direction = np.random.random(3)*2.0-1
for rot in [Rotation.from_random() for r in range(n//100)]:
R = damask.Orientation(rot,lattice)
color = R.IPFcolor(direction)
for equivalent in R.equivalentOrientations():
assert np.allclose(color,R.IPFcolor(direction))
for ori in Orientation(Rotation(set_of_quaternions),lattice)[:200]:
color = ori.IPF_color(direction)
for equivalent in ori.equivalent:
assert np.allclose(color,equivalent.IPF_color(direction))
@pytest.mark.parametrize('lattice',Lattice.lattices)
def test_reduced(self,set_of_quaternions,lattice):
oris = Orientation(Rotation(set_of_quaternions),lattice)
reduced = oris.reduced
assert np.all(reduced.in_FZ) and oris.rotation.shape == reduced.rotation.shape
@pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch'])
@pytest.mark.parametrize('lattice',['fcc','bcc'])
def test_relationship_forward_backward(self,model,lattice):
ori = Orientation(Rotation.from_random(),lattice)
for i,r in enumerate(ori.relatedOrientations(model)):
ori2 = r.relatedOrientations(model)[i]
for i,r in enumerate(ori.related(model)):
ori2 = r.related(model)[i]
misorientation = ori.rotation.misorientation(ori2.rotation)
assert misorientation.asAxisAngle(degrees=True)[3]<1.0e-5
assert misorientation.as_axis_angle(degrees=True)[3]<1.0e-5
@pytest.mark.parametrize('model',['Bain','KS','GT','GT_prime','NW','Pitsch'])
@pytest.mark.parametrize('lattice',['fcc','bcc'])
def test_relationship_reference(self,update,reference_dir,model,lattice):
reference = os.path.join(reference_dir,f'{lattice}_{model}.txt')
ori = Orientation(Rotation(),lattice)
eu = np.array([o.rotation.as_Eulers(degrees=True) for o in ori.relatedOrientations(model)])
eu = np.array([o.rotation.as_Eulers(degrees=True) for o in ori.related(model)])
if update:
coords = np.array([(1,i+1) for i,x in enumerate(eu)])
table = damask.Table(eu,{'Eulers':(3,)})
table.add('pos',coords)
table.to_ASCII(reference)
assert np.allclose(eu,damask.Table.from_ASCII(reference).get('Eulers'))
@pytest.mark.parametrize('lattice',Lattice.lattices)
def test_disorientation360(self,lattice):
R_1 = Orientation(Rotation(),lattice)
R_2 = Orientation(damask.Rotation.from_Eulers([360,0,0],degrees=True),lattice)
assert np.allclose(R_1.disorientation(R_2).as_matrix(),np.eye(3))
@pytest.mark.parametrize('lattice',Lattice.lattices)
@pytest.mark.parametrize('angle',[10,20,30,40])
def test_average(self,angle,lattice):
R_1 = Orientation(Rotation.from_axis_angle([0,0,1,10],degrees=True),lattice)
R_2 = Orientation(Rotation.from_axis_angle([0,0,1,angle],degrees=True),lattice)
avg_angle = R_1.average(R_2).rotation.as_axis_angle(degrees=True,pair=True)[1]
assert np.isclose(avg_angle,10+(angle-10)/2.)
@pytest.mark.parametrize('lattice',Lattice.lattices)
def test_from_average(self,lattice):
R_1 = Orientation(Rotation.from_random(),lattice)
eqs = [r for r in R_1.equivalent]
R_2 = damask.Orientation.from_average(eqs)
assert np.allclose(R_1.rotation.quaternion,R_2.rotation.quaternion)

10
python/tests/test_Result.py
View File

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

187
python/tests/test_Rotation.py
View File

@ -6,94 +6,21 @@ import numpy as np
from damask import Rotation
from damask import _rotation
n = 1100
atol=1.e-4
scatter=1.e-2
@pytest.fixture
def default():
"""A set of n rotations (corner cases and random)."""
specials = np.array([
[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,-1.0, 0.0, 0.0],
[0.0, 0.0,-1.0, 0.0],
[0.0, 0.0, 0.0,-1.0],
#----------------------
[1.0, 1.0, 0.0, 0.0],
[1.0, 0.0, 1.0, 0.0],
[1.0, 0.0, 0.0, 1.0],
[0.0, 1.0, 1.0, 0.0],
[0.0, 1.0, 0.0, 1.0],
[0.0, 0.0, 1.0, 1.0],
#----------------------
[1.0,-1.0, 0.0, 0.0],
[1.0, 0.0,-1.0, 0.0],
[1.0, 0.0, 0.0,-1.0],
[0.0, 1.0,-1.0, 0.0],
[0.0, 1.0, 0.0,-1.0],
[0.0, 0.0, 1.0,-1.0],
#----------------------
[0.0, 1.0,-1.0, 0.0],
[0.0, 1.0, 0.0,-1.0],
[0.0, 0.0, 1.0,-1.0],
#----------------------
[0.0,-1.0,-1.0, 0.0],
[0.0,-1.0, 0.0,-1.0],
[0.0, 0.0,-1.0,-1.0],
#----------------------
[1.0, 1.0, 1.0, 0.0],
[1.0, 1.0, 0.0, 1.0],
[1.0, 0.0, 1.0, 1.0],
[1.0,-1.0, 1.0, 0.0],
[1.0,-1.0, 0.0, 1.0],
[1.0, 0.0,-1.0, 1.0],
[1.0, 1.0,-1.0, 0.0],
[1.0, 1.0, 0.0,-1.0],
[1.0, 0.0, 1.0,-1.0],
[1.0,-1.0,-1.0, 0.0],
[1.0,-1.0, 0.0,-1.0],
[1.0, 0.0,-1.0,-1.0],
#----------------------
[0.0, 1.0, 1.0, 1.0],
[0.0, 1.0,-1.0, 1.0],
[0.0, 1.0, 1.0,-1.0],
[0.0,-1.0, 1.0, 1.0],
[0.0,-1.0,-1.0, 1.0],
[0.0,-1.0, 1.0,-1.0],
[0.0,-1.0,-1.0,-1.0],
#----------------------
[1.0, 1.0, 1.0, 1.0],
[1.0,-1.0, 1.0, 1.0],
[1.0, 1.0,-1.0, 1.0],
[1.0, 1.0, 1.0,-1.0],
[1.0,-1.0,-1.0, 1.0],
[1.0,-1.0, 1.0,-1.0],
[1.0, 1.0,-1.0,-1.0],
[1.0,-1.0,-1.0,-1.0],
])
specials /= np.linalg.norm(specials,axis=1).reshape(-1,1)
specials_scatter = specials + np.broadcast_to(np.random.rand(4)*scatter,specials.shape)
specials_scatter /= np.linalg.norm(specials_scatter,axis=1).reshape(-1,1)
specials_scatter[specials_scatter[:,0]<0]*=-1
return [Rotation.from_quaternion(s) for s in specials] + \
[Rotation.from_quaternion(s) for s in specials_scatter] + \
[Rotation.from_random() for _ in range(n-len(specials)-len(specials_scatter))]
@pytest.fixture
def reference_dir(reference_dir_base):
"""Directory containing reference results."""
return os.path.join(reference_dir_base,'Rotation')
@pytest.fixture
def set_of_rotations(set_of_quaternions):
return [Rotation.from_quaternion(s) for s in set_of_quaternions]
####################################################################################################
# Code below available according to the following conditions on https://github.com/MarDiehl/3Drotations
# Code below available according to the following conditions
####################################################################################################
# Copyright (c) 2017-2019, Martin Diehl/Max-Planck-Institut für Eisenforschung GmbH
# Copyright (c) 2013-2014, Marc De Graef/Carnegie Mellon University
@ -567,9 +494,9 @@ class TestRotation:
(Rotation._qu2ro,Rotation._ro2qu),
(Rotation._qu2ho,Rotation._ho2qu),
(Rotation._qu2cu,Rotation._cu2qu)])
def test_quaternion_internal(self,default,forward,backward):
def test_quaternion_internal(self,set_of_rotations,forward,backward):
"""Ensure invariance of conversion from quaternion and back."""
for rot in default:
for rot in set_of_rotations:
m = rot.as_quaternion()
o = backward(forward(m))
ok = np.allclose(m,o,atol=atol)
@ -584,9 +511,9 @@ class TestRotation:
(Rotation._om2ro,Rotation._ro2om),
(Rotation._om2ho,Rotation._ho2om),
(Rotation._om2cu,Rotation._cu2om)])
def test_matrix_internal(self,default,forward,backward):
def test_matrix_internal(self,set_of_rotations,forward,backward):
"""Ensure invariance of conversion from rotation matrix and back."""
for rot in default:
for rot in set_of_rotations:
m = rot.as_matrix()
o = backward(forward(m))
ok = np.allclose(m,o,atol=atol)
@ -599,9 +526,9 @@ class TestRotation:
(Rotation._eu2ro,Rotation._ro2eu),
(Rotation._eu2ho,Rotation._ho2eu),
(Rotation._eu2cu,Rotation._cu2eu)])
def test_Eulers_internal(self,default,forward,backward):
def test_Eulers_internal(self,set_of_rotations,forward,backward):
"""Ensure invariance of conversion from Euler angles and back."""
for rot in default:
for rot in set_of_rotations:
m = rot.as_Eulers()
o = backward(forward(m))
u = np.array([np.pi*2,np.pi,np.pi*2])
@ -619,9 +546,9 @@ class TestRotation:
(Rotation._ax2ro,Rotation._ro2ax),
(Rotation._ax2ho,Rotation._ho2ax),
(Rotation._ax2cu,Rotation._cu2ax)])
def test_axis_angle_internal(self,default,forward,backward):
def test_axis_angle_internal(self,set_of_rotations,forward,backward):
"""Ensure invariance of conversion from axis angle angles pair and back."""
for rot in default:
for rot in set_of_rotations:
m = rot.as_axis_angle()
o = backward(forward(m))
ok = np.allclose(m,o,atol=atol)
@ -636,10 +563,10 @@ class TestRotation:
(Rotation._ro2ax,Rotation._ax2ro),
(Rotation._ro2ho,Rotation._ho2ro),
(Rotation._ro2cu,Rotation._cu2ro)])
def test_Rodrigues_internal(self,default,forward,backward):
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))
for rot in default:
for rot in set_of_rotations:
m = rot.as_Rodrigues()
o = backward(forward(m))
ok = np.allclose(np.clip(m,None,cutoff),np.clip(o,None,cutoff),atol=atol)
@ -653,9 +580,9 @@ class TestRotation:
(Rotation._ho2ax,Rotation._ax2ho),
(Rotation._ho2ro,Rotation._ro2ho),
(Rotation._ho2cu,Rotation._cu2ho)])
def test_homochoric_internal(self,default,forward,backward):
def test_homochoric_internal(self,set_of_rotations,forward,backward):
"""Ensure invariance of conversion from homochoric vector and back."""
for rot in default:
for rot in set_of_rotations:
m = rot.as_homochoric()
o = backward(forward(m))
ok = np.allclose(m,o,atol=atol)
@ -668,9 +595,9 @@ class TestRotation:
(Rotation._cu2ax,Rotation._ax2cu),
(Rotation._cu2ro,Rotation._ro2cu),
(Rotation._cu2ho,Rotation._ho2cu)])
def test_cubochoric_internal(self,default,forward,backward):
def test_cubochoric_internal(self,set_of_rotations,forward,backward):
"""Ensure invariance of conversion from cubochoric vector and back."""
for rot in default:
for rot in set_of_rotations:
m = rot.as_cubochoric()
o = backward(forward(m))
ok = np.allclose(m,o,atol=atol)
@ -684,9 +611,9 @@ class TestRotation:
(Rotation._qu2ax,qu2ax),
(Rotation._qu2ro,qu2ro),
(Rotation._qu2ho,qu2ho)])
def test_quaternion_vectorization(self,default,vectorized,single):
def test_quaternion_vectorization(self,set_of_quaternions,vectorized,single):
"""Check vectorized implementation for quaternion against single point calculation."""
qu = np.array([rot.as_quaternion() for rot in default])
qu = np.array(set_of_quaternions)
vectorized(qu.reshape(qu.shape[0]//2,-1,4))
co = vectorized(qu)
for q,c in zip(qu,co):
@ -697,9 +624,9 @@ class TestRotation:
@pytest.mark.parametrize('vectorized, single',[(Rotation._om2qu,om2qu),
(Rotation._om2eu,om2eu),
(Rotation._om2ax,om2ax)])
def test_matrix_vectorization(self,default,vectorized,single):
def test_matrix_vectorization(self,set_of_rotations,vectorized,single):
"""Check vectorized implementation for rotation matrix against single point calculation."""
om = np.array([rot.as_matrix() for rot in default])
om = np.array([rot.as_matrix() for rot in set_of_rotations])
vectorized(om.reshape(om.shape[0]//2,-1,3,3))
co = vectorized(om)
for o,c in zip(om,co):
@ -710,9 +637,9 @@ class TestRotation:
(Rotation._eu2om,eu2om),
(Rotation._eu2ax,eu2ax),
(Rotation._eu2ro,eu2ro)])
def test_Eulers_vectorization(self,default,vectorized,single):
def test_Eulers_vectorization(self,set_of_rotations,vectorized,single):
"""Check vectorized implementation for Euler angles against single point calculation."""
eu = np.array([rot.as_Eulers() for rot in default])
eu = np.array([rot.as_Eulers() for rot in set_of_rotations])
vectorized(eu.reshape(eu.shape[0]//2,-1,3))
co = vectorized(eu)
for e,c in zip(eu,co):
@ -723,9 +650,9 @@ class TestRotation:
(Rotation._ax2om,ax2om),
(Rotation._ax2ro,ax2ro),
(Rotation._ax2ho,ax2ho)])
def test_axis_angle_vectorization(self,default,vectorized,single):
def test_axis_angle_vectorization(self,set_of_rotations,vectorized,single):
"""Check vectorized implementation for axis angle pair against single point calculation."""
ax = np.array([rot.as_axis_angle() for rot in default])
ax = np.array([rot.as_axis_angle() for rot in set_of_rotations])
vectorized(ax.reshape(ax.shape[0]//2,-1,4))
co = vectorized(ax)
for a,c in zip(ax,co):
@ -735,9 +662,9 @@ class TestRotation:
@pytest.mark.parametrize('vectorized, single',[(Rotation._ro2ax,ro2ax),
(Rotation._ro2ho,ro2ho)])
def test_Rodrigues_vectorization(self,default,vectorized,single):
def test_Rodrigues_vectorization(self,set_of_rotations,vectorized,single):
"""Check vectorized implementation for Rodrigues-Frank vector against single point calculation."""
ro = np.array([rot.as_Rodrigues() for rot in default])
ro = np.array([rot.as_Rodrigues() for rot in set_of_rotations])
vectorized(ro.reshape(ro.shape[0]//2,-1,4))
co = vectorized(ro)
for r,c in zip(ro,co):
@ -746,9 +673,9 @@ class TestRotation:
@pytest.mark.parametrize('vectorized, single',[(Rotation._ho2ax,ho2ax),
(Rotation._ho2cu,ho2cu)])
def test_homochoric_vectorization(self,default,vectorized,single):
def test_homochoric_vectorization(self,set_of_rotations,vectorized,single):
"""Check vectorized implementation for homochoric vector against single point calculation."""
ho = np.array([rot.as_homochoric() for rot in default])
ho = np.array([rot.as_homochoric() for rot in set_of_rotations])
vectorized(ho.reshape(ho.shape[0]//2,-1,3))
co = vectorized(ho)
for h,c in zip(ho,co):
@ -756,9 +683,9 @@ class TestRotation:
assert np.allclose(single(h),c) and np.allclose(single(h),vectorized(h))
@pytest.mark.parametrize('vectorized, single',[(Rotation._cu2ho,cu2ho)])
def test_cubochoric_vectorization(self,default,vectorized,single):
def test_cubochoric_vectorization(self,set_of_rotations,vectorized,single):
"""Check vectorized implementation for cubochoric vector against single point calculation."""
cu = np.array([rot.as_cubochoric() for rot in default])
cu = np.array([rot.as_cubochoric() for rot in set_of_rotations])
vectorized(cu.reshape(cu.shape[0]//2,-1,3))
co = vectorized(cu)
for u,c in zip(cu,co):
@ -766,8 +693,8 @@ class TestRotation:
assert np.allclose(single(u),c) and np.allclose(single(u),vectorized(u))
@pytest.mark.parametrize('degrees',[True,False])
def test_Eulers(self,default,degrees):
for rot in default:
def test_Eulers(self,set_of_rotations,degrees):
for rot in set_of_rotations:
m = rot.as_quaternion()
o = Rotation.from_Eulers(rot.as_Eulers(degrees),degrees).as_quaternion()
ok = np.allclose(m,o,atol=atol)
@ -779,9 +706,9 @@ class TestRotation:
@pytest.mark.parametrize('P',[1,-1])
@pytest.mark.parametrize('normalise',[True,False])
@pytest.mark.parametrize('degrees',[True,False])
def test_axis_angle(self,default,degrees,normalise,P):
def test_axis_angle(self,set_of_rotations,degrees,normalise,P):
c = np.array([P*-1,P*-1,P*-1,1.])
for rot in default:
for rot in set_of_rotations:
m = rot.as_Eulers()
o = Rotation.from_axis_angle(rot.as_axis_angle(degrees)*c,degrees,normalise,P).as_Eulers()
u = np.array([np.pi*2,np.pi,np.pi*2])
@ -793,8 +720,8 @@ class TestRotation:
print(m,o,rot.as_quaternion())
assert ok and (np.zeros(3)-1.e-9 <= o).all() and (o <= np.array([np.pi*2.,np.pi,np.pi*2.])+1.e-9).all()
def test_matrix(self,default):
for rot in default:
def test_matrix(self,set_of_rotations):
for rot in set_of_rotations:
m = rot.as_axis_angle()
o = Rotation.from_axis_angle(rot.as_axis_angle()).as_axis_angle()
ok = np.allclose(m,o,atol=atol)
@ -805,9 +732,9 @@ class TestRotation:
@pytest.mark.parametrize('P',[1,-1])
@pytest.mark.parametrize('normalise',[True,False])
def test_Rodrigues(self,default,normalise,P):
def test_Rodrigues(self,set_of_rotations,normalise,P):
c = np.array([P*-1,P*-1,P*-1,1.])
for rot in default:
for rot in set_of_rotations:
m = rot.as_matrix()
o = Rotation.from_Rodrigues(rot.as_Rodrigues()*c,normalise,P).as_matrix()
ok = np.allclose(m,o,atol=atol)
@ -815,9 +742,9 @@ class TestRotation:
assert ok and np.isclose(np.linalg.det(o),1.0)
@pytest.mark.parametrize('P',[1,-1])
def test_homochoric(self,default,P):
def test_homochoric(self,set_of_rotations,P):
cutoff = np.tan(np.pi*.5*(1.-1e-4))
for rot in default:
for rot in set_of_rotations:
m = rot.as_Rodrigues()
o = Rotation.from_homochoric(rot.as_homochoric()*P*-1,P).as_Rodrigues()
ok = np.allclose(np.clip(m,None,cutoff),np.clip(o,None,cutoff),atol=atol)
@ -826,8 +753,8 @@ class TestRotation:
assert ok and np.isclose(np.linalg.norm(o[:3]),1.0)
@pytest.mark.parametrize('P',[1,-1])
def test_cubochoric(self,default,P):
for rot in default:
def test_cubochoric(self,set_of_rotations,P):
for rot in set_of_rotations:
m = rot.as_homochoric()
o = Rotation.from_cubochoric(rot.as_cubochoric()*P*-1,P).as_homochoric()
ok = np.allclose(m,o,atol=atol)
@ -836,9 +763,9 @@ class TestRotation:
@pytest.mark.parametrize('P',[1,-1])
@pytest.mark.parametrize('accept_homomorph',[True,False])
def test_quaternion(self,default,P,accept_homomorph):
def test_quaternion(self,set_of_rotations,P,accept_homomorph):
c = np.array([1,P*-1,P*-1,P*-1]) * (-1 if accept_homomorph else 1)
for rot in default:
for rot in set_of_rotations:
m = rot.as_cubochoric()
o = Rotation.from_quaternion(rot.as_quaternion()*c,accept_homomorph,P).as_cubochoric()
ok = np.allclose(m,o,atol=atol)
@ -848,8 +775,8 @@ class TestRotation:
assert ok and o.max() < np.pi**(2./3.)*0.5+1.e-9
@pytest.mark.parametrize('reciprocal',[True,False])
def test_basis(self,default,reciprocal):
for rot in default:
def test_basis(self,set_of_rotations,reciprocal):
for rot in set_of_rotations:
om = rot.as_matrix() + 0.1*np.eye(3)
rot = Rotation.from_basis(om,False,reciprocal=reciprocal)
assert np.isclose(np.linalg.det(rot.as_matrix()),1.0)
@ -923,8 +850,8 @@ class TestRotation:
@pytest.mark.parametrize('data',[np.random.rand(5,3),
np.random.rand(5,3,3),
np.random.rand(5,3,3,3,3)])
def test_rotate_vectorization(self,default,data):
for rot in default:
def test_rotate_vectorization(self,set_of_rotations,data):
for rot in set_of_rotations:
v = rot.broadcast_to((5,)) @ data
for i in range(data.shape[0]):
print(i-data[i])
@ -978,3 +905,15 @@ class TestRotation:
def test_misorientation(self):
R = Rotation.from_random()
assert np.allclose(R.misorientation(R).as_matrix(),np.eye(3))
def test_misorientation360(self):
R_1 = Rotation()
R_2 = Rotation.from_Eulers([360,0,0],degrees=True)
assert np.allclose(R_1.misorientation(R_2).as_matrix(),np.eye(3))
@pytest.mark.parametrize('angle',[10,20,30,40,50,60,70,80,90,100,120])
def test_average(self,angle):
R_1 = Rotation.from_axis_angle([0,0,1,10],degrees=True)
R_2 = Rotation.from_axis_angle([0,0,1,angle],degrees=True)
avg_angle = R_1.average(R_2).as_axis_angle(degrees=True,pair=True)[1]
assert np.isclose(avg_angle,10+(angle-10)/2.)