diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 368888436..de2fa3906 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -207,7 +207,9 @@ Post_ParaviewRelated: Post_OrientationConversion: stage: postprocessing - script: OrientationConversion/test.py + script: + - OrientationConversion/test.py + - OrientationConversion/test2.py except: - master - release diff --git a/PRIVATE b/PRIVATE index f0090997d..8deb37dd4 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit f0090997df817f0a0b5a480a60e81929875b1010 +Subproject commit 8deb37dd4526fb5e1425fe1d2360508d01b6ac3e diff --git a/processing/post/addIPFcolor.py b/processing/post/addIPFcolor.py index 9c191b3ad..c5e4d8704 100755 --- a/processing/post/addIPFcolor.py +++ b/processing/post/addIPFcolor.py @@ -41,6 +41,10 @@ parser.set_defaults(pole = (0.0,0.0,1.0), (options, filenames) = parser.parse_args() +# damask.Orientation requires Bravais lattice, but we are only interested in symmetry +symmetry2lattice={'cubic':'bcc','hexagonal':'hex','tetragonal':'bct'} +lattice = symmetry2lattice[options.symmetry] + pole = np.array(options.pole) pole /= np.linalg.norm(pole) @@ -78,8 +82,8 @@ for name in filenames: outputAlive = True while outputAlive and table.data_read(): # read next data line of ASCII table - o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4]))), - symmetry = options.symmetry).reduced() + o = damask.Orientation(np.array(list(map(float,table.data[column:column+4]))), + lattice = lattice).reduced() table.data_append(o.IPFcolor(pole)) outputAlive = table.data_write() # output processed line diff --git a/processing/post/addPole.py b/processing/post/addPole.py index 628d64d5e..5116589b4 100755 --- a/processing/post/addPole.py +++ b/processing/post/addPole.py @@ -75,9 +75,9 @@ for name in filenames: # ------------------------------------------ process data ------------------------------------------ outputAlive = True while outputAlive and table.data_read(): # read next data line of ASCII table - o = damask.Orientation(quaternion = np.array(list(map(float,table.data[column:column+4])))) + o = damask.Rotation(np.array(list(map(float,table.data[column:column+4])))) - rotatedPole = o.quaternion*pole # rotate pole according to crystal orientation + rotatedPole = o*pole # rotate pole according to crystal orientation (x,y) = rotatedPole[0:2]/(1.+abs(pole[2])) # stereographic projection table.data_append([np.sqrt(x*x+y*y),np.arctan2(y,x)] if options.polar else [x,y]) # cartesian coordinates diff --git a/processing/post/rotateData.py b/processing/post/rotateData.py index 65f5aaaa2..ae42cb54a 100755 --- a/processing/post/rotateData.py +++ b/processing/post/rotateData.py @@ -31,7 +31,7 @@ parser.add_option('--degrees', action = 'store_true', help = 'angles are given in degrees') -parser.set_defaults(rotation = (0.,1.,0.,0.), # no rotation about 1,0,0 +parser.set_defaults(rotation = (0.,1.,1.,1.), # no rotation about 1,1,1 degrees = False, ) diff --git a/processing/pre/geom_fromTable.py b/processing/pre/geom_fromTable.py index c0c4cf4d1..ad598d5b1 100755 --- a/processing/pre/geom_fromTable.py +++ b/processing/pre/geom_fromTable.py @@ -52,11 +52,8 @@ parser.add_option('--crystallite', parser.set_defaults(symmetry = [damask.Symmetry.lattices[-1]], - tolerance = 0.0, - degrees = False, homogenization = 1, crystallite = 1, - verbose = False, pos = 'pos', ) diff --git a/python/damask/__init__.py b/python/damask/__init__.py index a9209a1c6..d7ed4a9f9 100644 --- a/python/damask/__init__.py +++ b/python/damask/__init__.py @@ -13,7 +13,7 @@ from .asciitable import ASCIItable # noqa from .config import Material # noqa from .colormaps import Colormap, Color # noqa -from .orientation import Quaternion, Symmetry, Rotation, Orientation # noqa +from .orientation import Symmetry, Lattice, Rotation, Orientation # noqa #from .block import Block # only one class from .result import Result # noqa diff --git a/python/damask/orientation.py b/python/damask/orientation.py index 558e5f15d..2f9731966 100644 --- a/python/damask/orientation.py +++ b/python/damask/orientation.py @@ -1,13 +1,13 @@ # -*- coding: UTF-8 no BOM -*- -import math,os +import math import numpy as np from . import Lambert P = -1 #################################################################################################### -class Quaternion2: +class Quaternion: u""" Quaternion with basic operations @@ -50,7 +50,7 @@ class Quaternion2: def __add__(self, other): """Addition""" - if isinstance(other, Quaternion2): + if isinstance(other, Quaternion): return self.__class__(q=self.q + other.q, p=self.p + other.p) else: @@ -58,7 +58,7 @@ class Quaternion2: def __iadd__(self, other): """In-place addition""" - if isinstance(other, Quaternion2): + if isinstance(other, Quaternion): self.q += other.q self.p += other.p return self @@ -72,7 +72,7 @@ class Quaternion2: def __sub__(self, other): """Subtraction""" - if isinstance(other, Quaternion2): + if isinstance(other, Quaternion): return self.__class__(q=self.q - other.q, p=self.p - other.p) else: @@ -80,7 +80,7 @@ class Quaternion2: def __isub__(self, other): """In-place subtraction""" - if isinstance(other, Quaternion2): + if isinstance(other, Quaternion): self.q -= other.q self.p -= other.p return self @@ -96,7 +96,7 @@ class Quaternion2: def __mul__(self, other): """Multiplication with quaternion or scalar""" - if isinstance(other, Quaternion2): + if isinstance(other, Quaternion): return self.__class__(q=self.q*other.q - np.dot(self.p,other.p), p=self.q*other.p + other.q*self.p + P * np.cross(self.p,other.p)) elif isinstance(other, (int, float)): @@ -107,7 +107,7 @@ class Quaternion2: def __imul__(self, other): """In-place multiplication with quaternion or scalar""" - if isinstance(other, Quaternion2): + if isinstance(other, Quaternion): self.q = self.q*other.q - np.dot(self.p,other.p) self.p = self.q*other.p + other.q*self.p + P * np.cross(self.p,other.p) return self @@ -120,7 +120,7 @@ class Quaternion2: def __truediv__(self, other): """Divsion with quaternion or scalar""" - if isinstance(other, Quaternion2): + if isinstance(other, Quaternion): s = other.conjugate()/abs(other)**2. return self.__class__(q=self.q * s, p=self.p * s) @@ -133,7 +133,7 @@ class Quaternion2: def __itruediv__(self, other): """In-place divsion with quaternion or scalar""" - if isinstance(other, Quaternion2): + if isinstance(other, Quaternion): s = other.conjugate()/abs(other)**2. self *= s return self @@ -215,7 +215,8 @@ class Rotation: u""" Orientation stored as unit quaternion. - All methods and naming conventions based on Rowenhorst_etal2015 + Following: D Rowenhorst et al. Consistent representations of and conversions between 3D rotations + 10.1088/0965-0393/23/8/083501 Convention 1: coordinate frames are right-handed Convention 2: a rotation angle ω is taken to be positive for a counterclockwise rotation when viewing from the end point of the rotation axis towards the origin @@ -242,10 +243,10 @@ class Rotation: If a quaternion is given, it needs to comply with the convection. Use .fromQuaternion to check the input. """ - if isinstance(quaternion,Quaternion2): + if isinstance(quaternion,Quaternion): self.quaternion = quaternion.copy() else: - self.quaternion = Quaternion2(q=quaternion[0],p=quaternion[1:4]) + self.quaternion = Quaternion(q=quaternion[0],p=quaternion[1:4]) self.quaternion.homomorph() # ToDo: Needed? def __repr__(self): @@ -453,374 +454,7 @@ class Rotation: def misorientation(self,other): """Misorientation""" return self.__class__(other.quaternion*self.quaternion.conjugated()) - - -# ****************************************************************************************** -class Quaternion: - u""" - Orientation represented as unit quaternion. - - All methods and naming conventions based on Rowenhorst_etal2015 - Convention 1: coordinate frames are right-handed - Convention 2: a rotation angle ω is taken to be positive for a counterclockwise rotation - when viewing from the end point of the rotation axis towards the origin - Convention 3: rotations will be interpreted in the passive sense - Convention 4: Euler angle triplets are implemented using the Bunge convention, - with the angular ranges as [0, 2π],[0, π],[0, 2π] - Convention 5: the rotation angle ω is limited to the interval [0, π] - Convention 6: P = -1 (as default) - - w is the real part, (x, y, z) are the imaginary parts. - - Vector "a" (defined in coordinate system "A") is passively rotated - resulting in new coordinates "b" when expressed in system "B". - b = Q * a - b = np.dot(Q.asMatrix(),a) - """ - - def __init__(self, - quat = None, - q = 1.0, - p = np.zeros(3,dtype=float)): - """Initializes to identity unless specified""" - self.q = quat[0] if quat is not None else q - self.p = np.array(quat[1:4]) if quat is not None else p - self.homomorph() - - def __iter__(self): - """Components""" - return iter(self.asList()) - - def __copy__(self): - """Copy""" - return self.__class__(q=self.q,p=self.p.copy()) - - copy = __copy__ - - def __repr__(self): - """Readable string""" - return 'Quaternion(real={q:+.6f}, imag=<{p[0]:+.6f}, {p[1]:+.6f}, {p[2]:+.6f}>)'.format(q=self.q,p=self.p) - - def __pow__(self, exponent): - """Power""" - omega = math.acos(self.q) - return self.__class__(q= math.cos(exponent*omega), - p=self.p * math.sin(exponent*omega)/math.sin(omega)) - - def __ipow__(self, exponent): - """In-place power""" - omega = math.acos(self.q) - self.q = math.cos(exponent*omega) - self.p *= math.sin(exponent*omega)/math.sin(omega) - return self - - def __mul__(self, other): - """Multiplication""" - # Rowenhorst_etal2015 MSMSE: value of P is selected as -1 - P = -1.0 - try: # quaternion - return self.__class__(q=self.q*other.q - np.dot(self.p,other.p), - p=self.q*other.p + other.q*self.p + P * np.cross(self.p,other.p)) - except: pass - try: # vector (perform passive rotation) - ( x, y, z) = self.p - (Vx,Vy,Vz) = other[0:3] - A = self.q*self.q - np.dot(self.p,self.p) - B = 2.0 * (x*Vx + y*Vy + z*Vz) - C = 2.0 * P*self.q - - return np.array([ - A*Vx + B*x + C*(y*Vz - z*Vy), - A*Vy + B*y + C*(z*Vx - x*Vz), - A*Vz + B*z + C*(x*Vy - y*Vx), - ]) - except: pass - try: # scalar - return self.__class__(q=self.q*other, - p=self.p*other) - except: - return self.copy() - - def __imul__(self, other): - """In-place multiplication""" - # Rowenhorst_etal2015 MSMSE: value of P is selected as -1 - P = -1.0 - try: # Quaternion - self.q = self.q*other.q - np.dot(self.p,other.p) - self.p = self.q*other.p + other.q*self.p + P * np.cross(self.p,other.p) - except: pass - return self - - def __div__(self, other): - """Division""" - if isinstance(other, (int,float)): - return self.__class__(q=self.q / other, - p=self.p / other) - else: - return NotImplemented - - def __idiv__(self, other): - """In-place division""" - if isinstance(other, (int,float)): - self.q /= other - self.p /= other - return self - - def __add__(self, other): - """Addition""" - if isinstance(other, Quaternion): - return self.__class__(q=self.q + other.q, - p=self.p + other.p) - else: - return NotImplemented - - def __iadd__(self, other): - """In-place addition""" - if isinstance(other, Quaternion): - self.q += other.q - self.p += other.p - return self - - def __sub__(self, other): - """Subtraction""" - if isinstance(other, Quaternion): - return self.__class__(q=self.q - other.q, - p=self.p - other.p) - else: - return NotImplemented - - def __isub__(self, other): - """In-place subtraction""" - if isinstance(other, Quaternion): - self.q -= other.q - self.p -= other.p - return self - - def __neg__(self): - """Additive inverse""" - self.q = -self.q - self.p = -self.p - return self - - def __abs__(self): - """Norm""" - return math.sqrt(self.q ** 2 + np.dot(self.p,self.p)) - - magnitude = __abs__ - - def __eq__(self,other): - """Equal (sufficiently close) to each other""" - return np.isclose(( self-other).magnitude(),0.0) \ - or np.isclose((-self-other).magnitude(),0.0) - - def __ne__(self,other): - """Not equal (sufficiently close) to each other""" - return not self.__eq__(other) - - def __cmp__(self,other): - """Linear ordering""" - return (1 if np.linalg.norm(self.asRodrigues()) > np.linalg.norm(other.asRodrigues()) else 0) \ - - (1 if np.linalg.norm(self.asRodrigues()) < np.linalg.norm(other.asRodrigues()) else 0) - - def magnitude_squared(self): - return self.q ** 2 + np.dot(self.p,self.p) - - def normalize(self): - d = self.magnitude() - if d > 0.0: - self.q /= d - self.p /= d - return self - - def conjugate(self): - self.p = -self.p - return self - - def homomorph(self): - if self.q < 0.0: - self.q = -self.q - self.p = -self.p - return self - - def normalized(self): - return self.copy().normalize() - - def conjugated(self): - return self.copy().conjugate() - - def homomorphed(self): - return self.copy().homomorph() - - def asList(self): - return [self.q]+list(self.p) - - def asM(self): # to find Averaging Quaternions (see F. Landis Markley et al.) - return np.outer(self.asList(),self.asList()) - def asMatrix(self): - # Rowenhorst_etal2015 MSMSE: value of P is selected as -1 - P = -1.0 - qbarhalf = 0.5*(self.q**2 - np.dot(self.p,self.p)) - return 2.0*np.array( - [[ qbarhalf + self.p[0]**2 , - self.p[0]*self.p[1] -P* self.q*self.p[2], - self.p[0]*self.p[2] +P* self.q*self.p[1] ], - [ self.p[0]*self.p[1] +P* self.q*self.p[2], - qbarhalf + self.p[1]**2 , - self.p[1]*self.p[2] -P* self.q*self.p[0] ], - [ self.p[0]*self.p[2] -P* self.q*self.p[1], - self.p[1]*self.p[2] +P* self.q*self.p[0], - qbarhalf + self.p[2]**2 ], - ]) - - def asAngleAxis(self, - degrees = False, - flat = False): - - angle = 2.0*math.acos(self.q) - - if np.isclose(angle,0.0): - angle = 0.0 - axis = np.array([0.0,0.0,1.0]) - elif np.isclose(self.q,0.0): - angle = math.pi - axis = self.p - else: - axis = np.sign(self.q)*self.p/np.linalg.norm(self.p) - - angle = np.degrees(angle) if degrees else angle - - return np.hstack((angle,axis)) if flat else (angle,axis) - - def asRodrigues(self): - return np.inf*np.ones(3) if np.isclose(self.q,0.0) else self.p/self.q - - -# # Static constructors - @classmethod - def fromIdentity(cls): - return cls() - - - @classmethod - def fromRandom(cls,randomSeed = None): - import binascii - if randomSeed is None: - randomSeed = int(binascii.hexlify(os.urandom(4)),16) - np.random.seed(randomSeed) - r = np.random.random(3) - A = math.sqrt(max(0.0,r[2])) - B = math.sqrt(max(0.0,1.0-r[2])) - w = math.cos(2.0*math.pi*r[0])*A - x = math.sin(2.0*math.pi*r[1])*B - y = math.cos(2.0*math.pi*r[1])*B - z = math.sin(2.0*math.pi*r[0])*A - return cls(quat=[w,x,y,z]) - - - @classmethod - def fromRodrigues(cls, rodrigues): - if not isinstance(rodrigues, np.ndarray): rodrigues = np.array(rodrigues) - norm = np.linalg.norm(rodrigues) - halfangle = math.atan(norm) - s = math.sin(halfangle) - c = math.cos(halfangle) - return cls(q=c,p=s*rodrigues/norm) - - - @classmethod - def fromAngleAxis(cls, - angle, - axis, - degrees = False): - if not isinstance(axis, np.ndarray): axis = np.array(axis,dtype=float) - axis = axis.astype(float)/np.linalg.norm(axis) - angle = np.radians(angle) if degrees else angle - s = math.sin(0.5 * angle) - c = math.cos(0.5 * angle) - return cls(q=c,p=axis*s) - - - @classmethod - def fromEulers(cls, - eulers, - degrees = False): - if not isinstance(eulers, np.ndarray): eulers = np.array(eulers,dtype=float) - eulers = np.radians(eulers) if degrees else eulers - - sigma = 0.5*(eulers[0]+eulers[2]) - delta = 0.5*(eulers[0]-eulers[2]) - c = np.cos(0.5*eulers[1]) - s = np.sin(0.5*eulers[1]) - - # Rowenhorst_etal2015 MSMSE: value of P is selected as -1 - P = -1.0 - w = c * np.cos(sigma) - x = -P * s * np.cos(delta) - y = -P * s * np.sin(delta) - z = -P * c * np.sin(sigma) - return cls(quat=[w,x,y,z]) - - -# Modified Method to calculate Quaternion from Orientation Matrix, -# Source: http://www.euclideanspace.com/maths/geometry/rotations/conversions/matrixToQuaternion/ - - @classmethod - def fromMatrix(cls, m): - if m.shape != (3,3) and np.prod(m.shape) == 9: - m = m.reshape(3,3) - - # Rowenhorst_etal2015 MSMSE: value of P is selected as -1 - P = -1.0 - w = 0.5*math.sqrt(max(0.0,1.0+m[0,0]+m[1,1]+m[2,2])) - x = P*0.5*math.sqrt(max(0.0,1.0+m[0,0]-m[1,1]-m[2,2])) - y = P*0.5*math.sqrt(max(0.0,1.0-m[0,0]+m[1,1]-m[2,2])) - z = P*0.5*math.sqrt(max(0.0,1.0-m[0,0]-m[1,1]+m[2,2])) - - x *= -1 if m[2,1] < m[1,2] else 1 - y *= -1 if m[0,2] < m[2,0] else 1 - z *= -1 if m[1,0] < m[0,1] else 1 - - return cls(quat=np.array([w,x,y,z])/math.sqrt(w**2 + x**2 + y**2 + z**2)) - - - @classmethod - def new_interpolate(cls, q1, q2, t): - """ - Interpolation - - See http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/20070017872_2007014421.pdf - for (another?) way to interpolate quaternions. - """ - assert isinstance(q1, Quaternion) and isinstance(q2, Quaternion) - Q = cls() - - costheta = q1.q*q2.q + np.dot(q1.p,q2.p) - if costheta < 0.: - costheta = -costheta - q1 = q1.conjugated() - elif costheta > 1.: - costheta = 1. - - theta = math.acos(costheta) - if abs(theta) < 0.01: - Q.q = q2.q - Q.p = q2.p - return Q - - sintheta = math.sqrt(1.0 - costheta * costheta) - if abs(sintheta) < 0.01: - Q.q = (q1.q + q2.q) * 0.5 - Q.p = (q1.p + q2.p) * 0.5 - return Q - - ratio1 = math.sin((1.0 - t) * theta) / sintheta - ratio2 = math.sin( t * theta) / sintheta - - Q.q = q1.q * ratio1 + q2.q * ratio2 - Q.p = q1.p * ratio1 + q2.p * ratio2 - return Q - # ****************************************************************************************** class Symmetry: @@ -932,26 +566,16 @@ class Symmetry: [ 1.0,0.0,0.0,0.0 ], ] - return list(map(Quaternion, - np.array(symQuats)[np.atleast_1d(np.array(who)) if who != [] else range(len(symQuats))])) + return np.array(symQuats) - - def equivalentQuaternions(self, - quaternion, - who = []): - """List of symmetrically equivalent quaternions based on own symmetry.""" - return [q*quaternion for q in self.symmetryQuats(who)] - def inFZ(self,R): - """Check whether given Rodrigues vector falls into fundamental zone of own symmetry.""" - if isinstance(R, Quaternion): R = R.asRodrigues() # translate accidentally passed quaternion -# fundamental zone in Rodrigues space is point symmetric around origin - - if R.shape[0]==4: # transition old (length not stored separately) to new - Rabs = abs(R[0:3]*R[3]) - else: - Rabs = abs(R) + """ + Check whether given Rodrigues vector falls into fundamental zone of own symmetry. + + Fundamental zone in Rodrigues space is point symmetric around origin. + """ + Rabs = abs(R[0:3]*R[3]) if self.lattice == 'cubic': return math.sqrt(2.0)-1.0 >= Rabs[0] \ @@ -1184,181 +808,181 @@ class Lattice: # from Y. He et al./Journal of Applied Crystallography (2006). 39, 72-81 GT = {'mapping':{'fcc':0,'bcc':1}, 'planes': np.array([ - [[ 1, 1, 1],[ 1, 0, 1]], - [[ 1, 1, 1],[ 1, 1, 0]], - [[ 1, 1, 1],[ 0, 1, 1]], - [[ -1, -1, 1],[ -1, 0, 1]], - [[ -1, -1, 1],[ -1, -1, 0]], - [[ -1, -1, 1],[ 0, -1, 1]], - [[ -1, 1, 1],[ -1, 0, 1]], - [[ -1, 1, 1],[ -1, 1, 0]], - [[ -1, 1, 1],[ 0, 1, 1]], - [[ 1, -1, 1],[ 1, 0, 1]], - [[ 1, -1, 1],[ 1, -1, 0]], - [[ 1, -1, 1],[ 0, -1, 1]], - [[ 1, 1, 1],[ 1, 1, 0]], - [[ 1, 1, 1],[ 0, 1, 1]], - [[ 1, 1, 1],[ 1, 0, 1]], - [[ -1, -1, 1],[ -1, -1, 0]], - [[ -1, -1, 1],[ 0, -1, 1]], - [[ -1, -1, 1],[ -1, 0, 1]], - [[ -1, 1, 1],[ -1, 1, 0]], - [[ -1, 1, 1],[ 0, 1, 1]], - [[ -1, 1, 1],[ -1, 0, 1]], - [[ 1, -1, 1],[ 1, -1, 0]], - [[ 1, -1, 1],[ 0, -1, 1]], - [[ 1, -1, 1],[ 1, 0, 1]]],dtype='float'), + [[ 1, 1, 1],[ 1, 0, 1]], + [[ 1, 1, 1],[ 1, 1, 0]], + [[ 1, 1, 1],[ 0, 1, 1]], + [[ -1, -1, 1],[ -1, 0, 1]], + [[ -1, -1, 1],[ -1, -1, 0]], + [[ -1, -1, 1],[ 0, -1, 1]], + [[ -1, 1, 1],[ -1, 0, 1]], + [[ -1, 1, 1],[ -1, 1, 0]], + [[ -1, 1, 1],[ 0, 1, 1]], + [[ 1, -1, 1],[ 1, 0, 1]], + [[ 1, -1, 1],[ 1, -1, 0]], + [[ 1, -1, 1],[ 0, -1, 1]], + [[ 1, 1, 1],[ 1, 1, 0]], + [[ 1, 1, 1],[ 0, 1, 1]], + [[ 1, 1, 1],[ 1, 0, 1]], + [[ -1, -1, 1],[ -1, -1, 0]], + [[ -1, -1, 1],[ 0, -1, 1]], + [[ -1, -1, 1],[ -1, 0, 1]], + [[ -1, 1, 1],[ -1, 1, 0]], + [[ -1, 1, 1],[ 0, 1, 1]], + [[ -1, 1, 1],[ -1, 0, 1]], + [[ 1, -1, 1],[ 1, -1, 0]], + [[ 1, -1, 1],[ 0, -1, 1]], + [[ 1, -1, 1],[ 1, 0, 1]]],dtype='float'), 'directions': np.array([ - [[ -5,-12, 17],[-17, -7, 17]], - [[ 17, -5,-12],[ 17,-17, -7]], - [[-12, 17, -5],[ -7, 17,-17]], - [[ 5, 12, 17],[ 17, 7, 17]], - [[-17, 5,-12],[-17, 17, -7]], - [[ 12,-17, -5],[ 7,-17,-17]], - [[ -5, 12,-17],[-17, 7,-17]], - [[ 17, 5, 12],[ 17, 17, 7]], - [[-12,-17, 5],[ -7,-17, 17]], - [[ 5,-12,-17],[ 17, -7,-17]], - [[-17, -5, 12],[-17,-17, 7]], - [[ 12, 17, 5],[ 7, 17, 17]], - [[ -5, 17,-12],[-17, 17, -7]], - [[-12, -5, 17],[ -7,-17, 17]], - [[ 17,-12, -5],[ 17, -7,-17]], - [[ 5,-17,-12],[ 17,-17, -7]], - [[ 12, 5, 17],[ 7, 17, 17]], - [[-17, 12, -5],[-17, 7,-17]], - [[ -5,-17, 12],[-17,-17, 7]], - [[-12, 5,-17],[ -7, 17,-17]], - [[ 17, 12, 5],[ 17, 7, 17]], - [[ 5, 17, 12],[ 17, 17, 7]], - [[ 12, -5,-17],[ 7,-17,-17]], - [[-17,-12, 5],[-17, 7, 17]]],dtype='float')} + [[ -5,-12, 17],[-17, -7, 17]], + [[ 17, -5,-12],[ 17,-17, -7]], + [[-12, 17, -5],[ -7, 17,-17]], + [[ 5, 12, 17],[ 17, 7, 17]], + [[-17, 5,-12],[-17, 17, -7]], + [[ 12,-17, -5],[ 7,-17,-17]], + [[ -5, 12,-17],[-17, 7,-17]], + [[ 17, 5, 12],[ 17, 17, 7]], + [[-12,-17, 5],[ -7,-17, 17]], + [[ 5,-12,-17],[ 17, -7,-17]], + [[-17, -5, 12],[-17,-17, 7]], + [[ 12, 17, 5],[ 7, 17, 17]], + [[ -5, 17,-12],[-17, 17, -7]], + [[-12, -5, 17],[ -7,-17, 17]], + [[ 17,-12, -5],[ 17, -7,-17]], + [[ 5,-17,-12],[ 17,-17, -7]], + [[ 12, 5, 17],[ 7, 17, 17]], + [[-17, 12, -5],[-17, 7,-17]], + [[ -5,-17, 12],[-17,-17, 7]], + [[-12, 5,-17],[ -7, 17,-17]], + [[ 17, 12, 5],[ 17, 7, 17]], + [[ 5, 17, 12],[ 17, 17, 7]], + [[ 12, -5,-17],[ 7,-17,-17]], + [[-17,-12, 5],[-17, 7, 17]]],dtype='float')} # Greninger--Troiano' orientation relationship for fcc <-> bcc transformation # from Y. He et al./Journal of Applied Crystallography (2006). 39, 72-81 GTdash = {'mapping':{'fcc':0,'bcc':1}, 'planes': np.array([ - [[ 7, 17, 17],[ 12, 5, 17]], - [[ 17, 7, 17],[ 17, 12, 5]], - [[ 17, 17, 7],[ 5, 17, 12]], - [[ -7,-17, 17],[-12, -5, 17]], - [[-17, -7, 17],[-17,-12, 5]], - [[-17,-17, 7],[ -5,-17, 12]], - [[ 7,-17,-17],[ 12, -5,-17]], - [[ 17, -7,-17],[ 17,-12, -5]], - [[ 17,-17, -7],[ 5,-17,-12]], - [[ -7, 17,-17],[-12, 5,-17]], - [[-17, 7,-17],[-17, 12, -5]], - [[-17, 17, -7],[ -5, 17,-12]], - [[ 7, 17, 17],[ 12, 17, 5]], - [[ 17, 7, 17],[ 5, 12, 17]], - [[ 17, 17, 7],[ 17, 5, 12]], - [[ -7,-17, 17],[-12,-17, 5]], - [[-17, -7, 17],[ -5,-12, 17]], - [[-17,-17, 7],[-17, -5, 12]], - [[ 7,-17,-17],[ 12,-17, -5]], - [[ 17, -7,-17],[ 5, -12,-17]], - [[ 17,-17, 7],[ 17, -5,-12]], - [[ -7, 17,-17],[-12, 17, -5]], - [[-17, 7,-17],[ -5, 12,-17]], - [[-17, 17, -7],[-17, 5,-12]]],dtype='float'), + [[ 7, 17, 17],[ 12, 5, 17]], + [[ 17, 7, 17],[ 17, 12, 5]], + [[ 17, 17, 7],[ 5, 17, 12]], + [[ -7,-17, 17],[-12, -5, 17]], + [[-17, -7, 17],[-17,-12, 5]], + [[-17,-17, 7],[ -5,-17, 12]], + [[ 7,-17,-17],[ 12, -5,-17]], + [[ 17, -7,-17],[ 17,-12, -5]], + [[ 17,-17, -7],[ 5,-17,-12]], + [[ -7, 17,-17],[-12, 5,-17]], + [[-17, 7,-17],[-17, 12, -5]], + [[-17, 17, -7],[ -5, 17,-12]], + [[ 7, 17, 17],[ 12, 17, 5]], + [[ 17, 7, 17],[ 5, 12, 17]], + [[ 17, 17, 7],[ 17, 5, 12]], + [[ -7,-17, 17],[-12,-17, 5]], + [[-17, -7, 17],[ -5,-12, 17]], + [[-17,-17, 7],[-17, -5, 12]], + [[ 7,-17,-17],[ 12,-17, -5]], + [[ 17, -7,-17],[ 5, -12,-17]], + [[ 17,-17, 7],[ 17, -5,-12]], + [[ -7, 17,-17],[-12, 17, -5]], + [[-17, 7,-17],[ -5, 12,-17]], + [[-17, 17, -7],[-17, 5,-12]]],dtype='float'), 'directions': np.array([ - [[ 0, 1, -1],[ 1, 1, -1]], - [[ -1, 0, 1],[ -1, 1, 1]], - [[ 1, -1, 0],[ 1, -1, 1]], - [[ 0, -1, -1],[ -1, -1, -1]], - [[ 1, 0, 1],[ 1, -1, 1]], - [[ 1, -1, 0],[ 1, -1, -1]], - [[ 0, 1, -1],[ -1, 1, -1]], - [[ 1, 0, 1],[ 1, 1, 1]], - [[ -1, -1, 0],[ -1, -1, 1]], - [[ 0, -1, -1],[ 1, -1, -1]], - [[ -1, 0, 1],[ -1, -1, 1]], - [[ -1, -1, 0],[ -1, -1, -1]], - [[ 0, -1, 1],[ 1, -1, 1]], - [[ 1, 0, -1],[ 1, 1, -1]], - [[ -1, 1, 0],[ -1, 1, 1]], - [[ 0, 1, 1],[ -1, 1, 1]], - [[ -1, 0, -1],[ -1, -1, -1]], - [[ -1, 1, 0],[ -1, 1, -1]], - [[ 0, -1, 1],[ -1, -1, 1]], - [[ -1, 0, -1],[ -1, 1, -1]], - [[ 1, 1, 0],[ 1, 1, 1]], - [[ 0, 1, 1],[ 1, 1, 1]], - [[ 1, 0, -1],[ 1, -1, -1]], - [[ 1, 1, 0],[ 1, 1, -1]]],dtype='float')} + [[ 0, 1, -1],[ 1, 1, -1]], + [[ -1, 0, 1],[ -1, 1, 1]], + [[ 1, -1, 0],[ 1, -1, 1]], + [[ 0, -1, -1],[ -1, -1, -1]], + [[ 1, 0, 1],[ 1, -1, 1]], + [[ 1, -1, 0],[ 1, -1, -1]], + [[ 0, 1, -1],[ -1, 1, -1]], + [[ 1, 0, 1],[ 1, 1, 1]], + [[ -1, -1, 0],[ -1, -1, 1]], + [[ 0, -1, -1],[ 1, -1, -1]], + [[ -1, 0, 1],[ -1, -1, 1]], + [[ -1, -1, 0],[ -1, -1, -1]], + [[ 0, -1, 1],[ 1, -1, 1]], + [[ 1, 0, -1],[ 1, 1, -1]], + [[ -1, 1, 0],[ -1, 1, 1]], + [[ 0, 1, 1],[ -1, 1, 1]], + [[ -1, 0, -1],[ -1, -1, -1]], + [[ -1, 1, 0],[ -1, 1, -1]], + [[ 0, -1, 1],[ -1, -1, 1]], + [[ -1, 0, -1],[ -1, 1, -1]], + [[ 1, 1, 0],[ 1, 1, 1]], + [[ 0, 1, 1],[ 1, 1, 1]], + [[ 1, 0, -1],[ 1, -1, -1]], + [[ 1, 1, 0],[ 1, 1, -1]]],dtype='float')} # Nishiyama--Wassermann orientation relationship for fcc <-> bcc transformation # from H. Kitahara et al./Materials Characterization 54 (2005) 378-386 NW = {'mapping':{'fcc':0,'bcc':1}, 'planes': np.array([ - [[ 1, 1, 1],[ 0, 1, 1]], - [[ 1, 1, 1],[ 0, 1, 1]], - [[ 1, 1, 1],[ 0, 1, 1]], - [[ -1, 1, 1],[ 0, 1, 1]], - [[ -1, 1, 1],[ 0, 1, 1]], - [[ -1, 1, 1],[ 0, 1, 1]], - [[ 1, -1, 1],[ 0, 1, 1]], - [[ 1, -1, 1],[ 0, 1, 1]], - [[ 1, -1, 1],[ 0, 1, 1]], - [[ -1, -1, 1],[ 0, 1, 1]], - [[ -1, -1, 1],[ 0, 1, 1]], - [[ -1, -1, 1],[ 0, 1, 1]]],dtype='float'), + [[ 1, 1, 1],[ 0, 1, 1]], + [[ 1, 1, 1],[ 0, 1, 1]], + [[ 1, 1, 1],[ 0, 1, 1]], + [[ -1, 1, 1],[ 0, 1, 1]], + [[ -1, 1, 1],[ 0, 1, 1]], + [[ -1, 1, 1],[ 0, 1, 1]], + [[ 1, -1, 1],[ 0, 1, 1]], + [[ 1, -1, 1],[ 0, 1, 1]], + [[ 1, -1, 1],[ 0, 1, 1]], + [[ -1, -1, 1],[ 0, 1, 1]], + [[ -1, -1, 1],[ 0, 1, 1]], + [[ -1, -1, 1],[ 0, 1, 1]]],dtype='float'), 'directions': np.array([ - [[ 2, -1, -1],[ 0, -1, 1]], - [[ -1, 2, -1],[ 0, -1, 1]], - [[ -1, -1, 2],[ 0, -1, 1]], - [[ -2, -1, -1],[ 0, -1, 1]], - [[ 1, 2, -1],[ 0, -1, 1]], - [[ 1, -1, 2],[ 0, -1, 1]], - [[ 2, 1, -1],[ 0, -1, 1]], - [[ -1, -2, -1],[ 0, -1, 1]], - [[ -1, 1, 2],[ 0, -1, 1]], - [[ -1, 2, 1],[ 0, -1, 1]], - [[ -1, 2, 1],[ 0, -1, 1]], - [[ -1, -1, -2],[ 0, -1, 1]]],dtype='float')} + [[ 2, -1, -1],[ 0, -1, 1]], + [[ -1, 2, -1],[ 0, -1, 1]], + [[ -1, -1, 2],[ 0, -1, 1]], + [[ -2, -1, -1],[ 0, -1, 1]], + [[ 1, 2, -1],[ 0, -1, 1]], + [[ 1, -1, 2],[ 0, -1, 1]], + [[ 2, 1, -1],[ 0, -1, 1]], + [[ -1, -2, -1],[ 0, -1, 1]], + [[ -1, 1, 2],[ 0, -1, 1]], + [[ -1, 2, 1],[ 0, -1, 1]], + [[ -1, 2, 1],[ 0, -1, 1]], + [[ -1, -1, -2],[ 0, -1, 1]]],dtype='float')} # Pitsch orientation relationship for fcc <-> bcc transformation # from Y. He et al./Acta Materialia 53 (2005) 1179-1190 Pitsch = {'mapping':{'fcc':0,'bcc':1}, 'planes': np.array([ - [[ 0, 1, 0],[ -1, 0, 1]], - [[ 0, 0, 1],[ 1, -1, 0]], - [[ 1, 0, 0],[ 0, 1, -1]], - [[ 1, 0, 0],[ 0, -1, -1]], - [[ 0, 1, 0],[ -1, 0, -1]], - [[ 0, 0, 1],[ -1, -1, 0]], - [[ 0, 1, 0],[ -1, 0, -1]], - [[ 0, 0, 1],[ -1, -1, 0]], - [[ 1, 0, 0],[ 0, -1, -1]], - [[ 1, 0, 0],[ 0, -1, 1]], - [[ 0, 1, 0],[ 1, 0, -1]], - [[ 0, 0, 1],[ -1, 1, 0]]],dtype='float'), + [[ 0, 1, 0],[ -1, 0, 1]], + [[ 0, 0, 1],[ 1, -1, 0]], + [[ 1, 0, 0],[ 0, 1, -1]], + [[ 1, 0, 0],[ 0, -1, -1]], + [[ 0, 1, 0],[ -1, 0, -1]], + [[ 0, 0, 1],[ -1, -1, 0]], + [[ 0, 1, 0],[ -1, 0, -1]], + [[ 0, 0, 1],[ -1, -1, 0]], + [[ 1, 0, 0],[ 0, -1, -1]], + [[ 1, 0, 0],[ 0, -1, 1]], + [[ 0, 1, 0],[ 1, 0, -1]], + [[ 0, 0, 1],[ -1, 1, 0]]],dtype='float'), 'directions': np.array([ - [[ 1, 0, 1],[ 1, -1, 1]], - [[ 1, 1, 0],[ 1, 1, -1]], - [[ 0, 1, 1],[ -1, 1, 1]], - [[ 0, 1, -1],[ -1, 1, -1]], - [[ -1, 0, 1],[ -1, -1, 1]], - [[ 1, -1, 0],[ 1, -1, -1]], - [[ 1, 0, -1],[ 1, -1, -1]], - [[ -1, 1, 0],[ -1, 1, -1]], - [[ 0, -1, 1],[ -1, -1, 1]], - [[ 0, 1, 1],[ -1, 1, 1]], - [[ 1, 0, 1],[ 1, -1, 1]], - [[ 1, 1, 0],[ 1, 1, -1]]],dtype='float')} + [[ 1, 0, 1],[ 1, -1, 1]], + [[ 1, 1, 0],[ 1, 1, -1]], + [[ 0, 1, 1],[ -1, 1, 1]], + [[ 0, 1, -1],[ -1, 1, -1]], + [[ -1, 0, 1],[ -1, -1, 1]], + [[ 1, -1, 0],[ 1, -1, -1]], + [[ 1, 0, -1],[ 1, -1, -1]], + [[ -1, 1, 0],[ -1, 1, -1]], + [[ 0, -1, 1],[ -1, -1, 1]], + [[ 0, 1, 1],[ -1, 1, 1]], + [[ 1, 0, 1],[ 1, -1, 1]], + [[ 1, 1, 0],[ 1, 1, -1]]],dtype='float')} # Bain orientation relationship for fcc <-> bcc transformation # from Y. He et al./Journal of Applied Crystallography (2006). 39, 72-81 Bain = {'mapping':{'fcc':0,'bcc':1}, 'planes': np.array([ - [[ 1, 0, 0],[ 1, 0, 0]], - [[ 0, 1, 0],[ 0, 1, 0]], - [[ 0, 0, 1],[ 0, 0, 1]]],dtype='float'), + [[ 1, 0, 0],[ 1, 0, 0]], + [[ 0, 1, 0],[ 0, 1, 0]], + [[ 0, 0, 1],[ 0, 0, 1]]],dtype='float'), 'directions': np.array([ - [[ 0, 1, 0],[ 0, 1, 1]], - [[ 0, 0, 1],[ 1, 0, 1]], - [[ 1, 0, 0],[ 1, 1, 0]]],dtype='float')} + [[ 0, 1, 0],[ 0, 1, 1]], + [[ 0, 0, 1],[ 1, 0, 1]], + [[ 1, 0, 0],[ 1, 1, 0]]],dtype='float')} def relationOperations(self,model): @@ -1367,27 +991,30 @@ class Lattice: relationship = models[model] - r = {'lattice':Lattice((set(relationship['mapping'])-{self.lattice}).pop()), + r = {'lattice':Lattice((set(relationship['mapping'])-{self.lattice}).pop()), # target lattice 'rotations':[] } myPlane_id = relationship['mapping'][self.lattice] otherPlane_id = (myPlane_id+1)%2 myDir_id = myPlane_id +2 otherDir_id = otherPlane_id +2 + for miller in np.hstack((relationship['planes'],relationship['directions'])): - myPlane = miller[myPlane_id]/ np.linalg.norm(miller[myPlane_id]) - myDir = miller[myDir_id]/ np.linalg.norm(miller[myDir_id]) - otherPlane = miller[otherPlane_id]/ np.linalg.norm(miller[otherPlane_id]) - otherDir = miller[otherDir_id]/ np.linalg.norm(miller[otherDir_id]) - - myMatrix = np.array([myDir,np.cross(myPlane,myDir),myPlane]).T - otherMatrix = np.array([otherDir,np.cross(otherPlane,otherDir),otherPlane]).T + myPlane = miller[myPlane_id]/ np.linalg.norm(miller[myPlane_id]) + myDir = miller[myDir_id]/ np.linalg.norm(miller[myDir_id]) + myMatrix = np.array([myDir,np.cross(myPlane,myDir),myPlane]).T + + otherPlane = miller[otherPlane_id]/ np.linalg.norm(miller[otherPlane_id]) + otherDir = miller[otherDir_id]/ np.linalg.norm(miller[otherDir_id]) + otherMatrix = np.array([otherDir,np.cross(otherPlane,otherDir),otherPlane]).T + r['rotations'].append(Rotation.fromMatrix(np.dot(otherMatrix,myMatrix.T))) + return r -class Orientation2: +class Orientation: """ Crystallographic orientation @@ -1426,32 +1053,28 @@ class Orientation2: # raise NotImplementedError('disorientation between different symmetry classes not supported yet.') mis = other.rotation*self.rotation.inversed() - mySymEqs = self.equivalentOrientations() if SST else self.equivalentOrientations()[:1] # take all or only first sym operation + mySymEqs = self.equivalentOrientations() if SST else self.equivalentOrientations()[:1] # take all or only first sym operation otherSymEqs = other.equivalentOrientations() for i,sA in enumerate(mySymEqs): for j,sB in enumerate(otherSymEqs): - theQ = sB.rotation*mis*sA.rotation.inversed() + r = sB.rotation*mis*sA.rotation.inversed() for k in range(2): - theQ.inversed() - breaker = self.lattice.symmetry.inFZ(theQ.asRodrigues()) \ - and (not SST or other.lattice.symmetry.inDisorientationSST(theQ.asRodrigues())) + r.inversed() + breaker = self.lattice.symmetry.inFZ(r.asRodrigues()) \ + and (not SST or other.lattice.symmetry.inDisorientationSST(r.asRodrigues())) if breaker: break if breaker: break if breaker: break -# disorientation, own sym, other sym, self-->other: True, self<--other: False - return theQ + return r def inFZ(self): return self.lattice.symmetry.inFZ(self.rotation.asRodrigues()) def equivalentOrientations(self): """List of orientations which are symmetrically equivalent""" - q = self.lattice.symmetry.symmetryQuats() - q2 = [Quaternion2(q=a.asList()[0],p=a.asList()[1:4]) for a in q] # convert Quaternion to Quaternion2 - x = [self.__class__(q3*self.rotation.quaternion,self.lattice) for q3 in q2] - return x + return [self.__class__(q*self.rotation.quaternion,self.lattice) for q in self.lattice.symmetry.symmetryQuats()] def relatedOrientations(self,model): """List of orientations related by the given orientation relationship""" @@ -1464,177 +1087,69 @@ class Orientation2: if self.lattice.symmetry.inFZ(me.rotation.asRodrigues()): break return self.__class__(me.rotation,self.lattice) - -# ****************************************************************************************** -class Orientation: - - __slots__ = ['quaternion','symmetry'] - - def __init__(self, - quaternion = Quaternion.fromIdentity(), - Rodrigues = None, - Eulers = None, - random = False, # integer to have a fixed seed or True for real random - symmetry = None, - degrees = False, - ): - if random: # produce random orientation - if isinstance(random, bool ): - self.quaternion = Quaternion.fromRandom() - else: - self.quaternion = Quaternion.fromRandom(randomSeed=random) - elif isinstance(Eulers, np.ndarray) and Eulers.shape == (3,): # based on given Euler angles - self.quaternion = Quaternion.fromEulers(Eulers,degrees=degrees) - elif isinstance(Rodrigues, np.ndarray) and Rodrigues.shape == (3,): # based on given Rodrigues vector - self.quaternion = Quaternion.fromRodrigues(Rodrigues) - elif isinstance(quaternion, Quaternion): # based on given quaternion - self.quaternion = quaternion.homomorphed() - elif (isinstance(quaternion, np.ndarray) and quaternion.shape == (4,)) or \ - (isinstance(quaternion, list) and len(quaternion) == 4 ): # based on given quaternion-like array - self.quaternion = Quaternion(quat=quaternion).homomorphed() - - self.symmetry = Symmetry(symmetry) - - def __copy__(self): - """Copy""" - return self.__class__(quaternion=self.quaternion,symmetry=self.symmetry.lattice) - - copy = __copy__ - - - def __repr__(self): - """Value as all implemented representations""" - return '\n'.join([ - 'Symmetry: {}'.format(self.symmetry), - 'Quaternion: {}'.format(self.quaternion), - 'Matrix:\n{}'.format( '\n'.join(['\t'.join(list(map(str,self.asMatrix()[i,:]))) for i in range(3)]) ), - ]) - - def asRodrigues(self): - return self.quaternion.asRodrigues() - rodrigues = property(asRodrigues) - - def asAngleAxis(self, - degrees = False, - flat = False): - return self.quaternion.asAngleAxis(degrees,flat) - - def asMatrix(self): - return self.quaternion.asMatrix() - matrix = property(asMatrix) - - def inFZ(self): - return self.symmetry.inFZ(self.quaternion.asRodrigues()) - - def equivalentQuaternions(self, - who = []): - return self.symmetry.equivalentQuaternions(self.quaternion,who) - - def equivalentOrientations(self, - who = []): - return [Orientation(quaternion = q, symmetry = self.symmetry.lattice) for q in self.equivalentQuaternions(who)] - - def reduced(self): - """Transform orientation to fall into fundamental zone according to symmetry""" - for me in self.symmetry.equivalentQuaternions(self.quaternion): - if self.symmetry.inFZ(me.asRodrigues()): break - - return Orientation(quaternion=me,symmetry=self.symmetry.lattice) - - - def disorientation(self, - other, - SST = True): - """ - Disorientation between myself and given other orientation. - - Rotation axis falls into SST if SST == True. - (Currently requires same symmetry for both orientations. - Look into A. Heinz and P. Neumann 1991 for cases with differing sym.) - """ - if self.symmetry != other.symmetry: - raise NotImplementedError('disorientation between different symmetry classes not supported yet.') - - misQ = other.quaternion*self.quaternion.conjugated() - mySymQs = self.symmetry.symmetryQuats() if SST else self.symmetry.symmetryQuats()[:1] # take all or only first sym operation - otherSymQs = other.symmetry.symmetryQuats() - for i,sA in enumerate(mySymQs): - for j,sB in enumerate(otherSymQs): - theQ = sB*misQ*sA.conjugated() - for k in range(2): - theQ.conjugate() - breaker = self.symmetry.inFZ(theQ) \ - and (not SST or other.symmetry.inDisorientationSST(theQ)) - if breaker: break - if breaker: break - if breaker: break - -# disorientation, own sym, other sym, self-->other: True, self<--other: False - return (Orientation(quaternion = theQ,symmetry = self.symmetry.lattice), - i,j, k == 1) - - def inversePole(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,q in enumerate(self.symmetry.equivalentQuaternions(self.quaternion)): # test all symmetric equivalent quaternions - pole = q*axis # align crystal direction to axis - if self.symmetry.inSST(pole,proper): break # found SST version + 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 else: - pole = self.quaternion*axis # align crystal direction to axis + pole = self.rotation*axis # align crystal direction to axis return (pole,i if SST else 0) - + + def IPFcolor(self,axis): """TSL color of inverse pole figure for given axis""" color = np.zeros(3,'d') - for q in self.symmetry.equivalentQuaternions(self.quaternion): - pole = q*axis # align crystal direction to axis - inSST,color = self.symmetry.inSST(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 - return color + return color - @classmethod - def average(cls, - orientations, - multiplicity = []): - """ - Average orientation - ref: F. Landis Markley, Yang Cheng, John Lucas Crassidis, and Yaakov Oshman. - Averaging Quaternions, - Journal of Guidance, Control, and Dynamics, Vol. 30, No. 4 (2007), pp. 1193-1197. - doi: 10.2514/1.28949 - usage: - a = Orientation(Eulers=np.radians([10, 10, 0]), symmetry='hexagonal') - b = Orientation(Eulers=np.radians([20, 0, 0]), symmetry='hexagonal') - avg = Orientation.average([a,b]) - """ - if not all(isinstance(item, Orientation) for item in orientations): - raise TypeError("Only instances of Orientation can be averaged.") + # @classmethod + # def average(cls, + # orientations, + # multiplicity = []): + # """ + # Average orientation - N = len(orientations) - if multiplicity == [] or not multiplicity: - multiplicity = np.ones(N,dtype='i') + # ref: F. Landis Markley, Yang Cheng, John Lucas Crassidis, and Yaakov Oshman. + # Averaging Quaternions, + # Journal of Guidance, Control, and Dynamics, Vol. 30, No. 4 (2007), pp. 1193-1197. + # doi: 10.2514/1.28949 + # usage: + # a = Orientation(Eulers=np.radians([10, 10, 0]), symmetry='hexagonal') + # b = Orientation(Eulers=np.radians([20, 0, 0]), symmetry='hexagonal') + # avg = Orientation.average([a,b]) + # """ + # if not all(isinstance(item, Orientation) for item in orientations): + # raise TypeError("Only instances of Orientation can be averaged.") - reference = orientations[0] # take first as reference - for i,(o,n) in enumerate(zip(orientations,multiplicity)): - closest = o.equivalentOrientations(reference.disorientation(o,SST = False)[2])[0] # select sym orientation with lowest misorientation - M = closest.quaternion.asM() * n if i == 0 else M + closest.quaternion.asM() * n # noqa add (multiples) of this orientation to average noqa - eig, vec = np.linalg.eig(M/N) + # N = len(orientations) + # if multiplicity == [] or not multiplicity: + # multiplicity = np.ones(N,dtype='i') - return Orientation(quaternion = Quaternion(quat = np.real(vec.T[eig.argmax()])), - symmetry = reference.symmetry.lattice) + # reference = orientations[0] # take first as reference + # for i,(o,n) in enumerate(zip(orientations,multiplicity)): + # closest = o.equivalentOrientations(reference.disorientation(o,SST = False)[2])[0] # select sym orientation with lowest misorientation + # M = closest.quaternion.asM() * n if i == 0 else M + closest.quaternion.asM() * n # noqa add (multiples) of this orientation to average noqa + # eig, vec = np.linalg.eig(M/N) + + # return Orientation(quaternion = Quaternion(quat = np.real(vec.T[eig.argmax()])), + # symmetry = reference.symmetry.lattice) #################################################################################################### -# Code below available according to the followin conditions on https://github.com/MarDiehl/3Drotations +# Code below available according to the following conditions on https://github.com/MarDiehl/3Drotations #################################################################################################### # Copyright (c) 2017-2019, Martin Diehl/Max-Planck-Institut für Eisenforschung GmbH # Copyright (c) 2013-2014, Marc De Graef/Carnegie Mellon University @@ -1725,7 +1240,7 @@ def eu2qu(eu): -P*sPhi*np.cos(ee[0]-ee[2]), -P*sPhi*np.sin(ee[0]-ee[2]), -P*cPhi*np.sin(ee[0]+ee[2]) ]) - #if qu[0] < 0.0: qu.homomorph() !ToDo: Check with original + #if qu[0] < 0.0: qu.homomorph() !ToDo: Check with original return qu @@ -1769,7 +1284,6 @@ def qu2eu(qu): eu = np.array([np.arctan2(-P*2.0*qu[0]*qu[3],qu[0]**2-qu[3]**2), 0.0, 0.0]) if iszero(q12) else \ np.array([np.arctan2(2.0*qu[1]*qu[2],qu[1]**2-qu[2]**2), np.pi, 0.0]) else: - #chiInv = 1.0/chi ToDo: needed for what? eu = np.array([np.arctan2((-P*qu[0]*qu[2]+qu[1]*qu[3])*chi, (-P*qu[0]*qu[1]-qu[2]*qu[3])*chi ), np.arctan2( 2.0*chi, q03-q12 ), np.arctan2(( P*qu[0]*qu[2]+qu[1]*qu[3])*chi, (-P*qu[0]*qu[1]+qu[2]*qu[3])*chi )])