diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 7b33574de..d82c7c8ac 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -213,6 +213,14 @@ Post_OrientationConversion: - master - release +Post_OrientationAverageMisorientation: + stage: postprocessing + script: + - OrientationAverageMisorientation/test.py + except: + - master + - release + ################################################################################################### grid_mech_compile_Intel: stage: compilePETSc diff --git a/VERSION b/VERSION index 9096149dc..c24732ae9 100644 --- a/VERSION +++ b/VERSION @@ -1 +1 @@ -v2.0.3-198-g1c762860 +v2.0.3-229-g90fa809f diff --git a/python/damask/orientation.py b/python/damask/orientation.py index 252f39420..b19037b0c 100644 --- a/python/damask/orientation.py +++ b/python/damask/orientation.py @@ -3,211 +3,8 @@ import math import numpy as np from . import Lambert - -P = -1 - -#################################################################################################### -class Quaternion: - u""" - Quaternion with basic operations - - q is the real part, p = (x, y, z) are the imaginary parts. - Defintion of multiplication depends on variable P, P ∉ {-1,1}. - """ - - def __init__(self, - q = 0.0, - p = np.zeros(3,dtype=float)): - """Initializes to identity unless specified""" - self.q = q - self.p = np.array(p) - - - def __copy__(self): - """Copy""" - return self.__class__(q=self.q, - p=self.p.copy()) - - copy = __copy__ - - - def __iter__(self): - """Components""" - return iter(self.asList()) - - def asArray(self): - """As numpy array""" - return np.array((self.q,self.p[0],self.p[1],self.p[2])) - - def asList(self): - return [self.q]+list(self.p) - - - 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 __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 - else: - return NotImplemented - - def __pos__(self): - """Unary positive operator""" - 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 - else: - return NotImplemented - - def __neg__(self): - """Unary positive operator""" - self.q *= -1.0 - self.p *= -1.0 - return self - - - def __mul__(self, other): - """Multiplication with quaternion or scalar""" - 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)): - return self.__class__(q=self.q*other, - p=self.p*other) - else: - return NotImplemented - - def __imul__(self, other): - """In-place multiplication with quaternion or scalar""" - 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 - elif isinstance(other, (int, float)): - self *= other - return self - else: - return NotImplemented - - - def __truediv__(self, other): - """Divsion with quaternion or scalar""" - if isinstance(other, Quaternion): - s = other.conjugate()/abs(other)**2. - return self.__class__(q=self.q * s, - p=self.p * s) - elif isinstance(other, (int, float)): - self.q /= other - self.p /= other - return self - else: - return NotImplemented - - def __itruediv__(self, other): - """In-place divsion with quaternion or scalar""" - if isinstance(other, Quaternion): - s = other.conjugate()/abs(other)**2. - self *= s - return self - elif isinstance(other, (int, float)): - self.q /= other - return self - else: - return NotImplemented - - - def __pow__(self, exponent): - """Power""" - if isinstance(exponent, (int, float)): - omega = np.acos(self.q) - return self.__class__(q= np.cos(exponent*omega), - p=self.p * np.sin(exponent*omega)/np.sin(omega)) - else: - return NotImplemented - - def __ipow__(self, exponent): - """In-place power""" - if isinstance(exponent, (int, float)): - omega = np.acos(self.q) - self.q = np.cos(exponent*omega) - self.p *= np.sin(exponent*omega)/np.sin(omega) - else: - return NotImplemented - - - 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 normalize(self): - d = self.magnitude() - if d > 0.0: - self.q /= d - self.p /= d - return self - - def normalized(self): - return self.copy().normalize() - - - def conjugate(self): - self.p = -self.p - return self - - def conjugated(self): - return self.copy().conjugate() - - - def homomorph(self): - if self.q < 0.0: - self.q = -self.q - self.p = -self.p - return self - - def homomorphed(self): - return self.copy().homomorph() - +from .quaternion import Quaternion +from .quaternion import P #################################################################################################### @@ -247,7 +44,13 @@ class Rotation: self.quaternion = quaternion.copy() else: self.quaternion = Quaternion(q=quaternion[0],p=quaternion[1:4]) - self.quaternion.homomorph() # ToDo: Needed? + + def __copy__(self): + """Copy""" + return self.__class__(self.quaternion) + + copy = __copy__ + def __repr__(self): """Value in selected representation""" @@ -257,6 +60,78 @@ class Rotation: 'Bunge Eulers / deg: {}'.format('\t'.join(list(map(str,self.asEulers(degrees=True)))) ), ]) + def __mul__(self, other): + """ + Multiplication + + Rotation: Details needed (active/passive), rotation of (3,3,3,3)-matrix should be considered + """ + if isinstance(other, Rotation): # rotate a rotation + return self.__class__(self.quaternion * other.quaternion).standardize() + elif isinstance(other, np.ndarray): + if other.shape == (3,): # rotate a single (3)-vector + ( x, y, z) = self.quaternion.p + (Vx,Vy,Vz) = other[0:3] + A = self.quaternion.q*self.quaternion.q - np.dot(self.quaternion.p,self.quaternion.p) + B = 2.0 * (x*Vx + y*Vy + z*Vz) + C = 2.0 * P*self.quaternion.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), + ]) + elif other.shape == (3,3,): # rotate a single (3x3)-matrix + return np.dot(self.asMatrix(),np.dot(other,self.asMatrix().T)) + elif other.shape == (3,3,3,3,): + raise NotImplementedError + else: + return NotImplemented + elif isinstance(other, tuple): # used to rotate a meshgrid-tuple + ( x, y, z) = self.quaternion.p + (Vx,Vy,Vz) = other[0:3] + A = self.quaternion.q*self.quaternion.q - np.dot(self.quaternion.p,self.quaternion.p) + B = 2.0 * (x*Vx + y*Vy + z*Vz) + C = 2.0 * P*self.quaternion.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), + ]) + else: + return NotImplemented + + + def inverse(self): + """In-place inverse rotation/backward rotation""" + self.quaternion.conjugate() + return self + + def inversed(self): + """Inverse rotation/backward rotation""" + return self.copy().inverse() + + + def standardize(self): + """In-place quaternion representation with positive q""" + if self.quaternion.q < 0.0: self.quaternion.homomorph() + return self + + def standardized(self): + """Quaternion representation with positive q""" + return self.copy().standardize() + + + def misorientation(self,other): + """Misorientation""" + return other*self.inversed() + + + def average(self,other): + """Calculate the average rotation""" + return Rotation.fromAverage([self,other]) + ################################################################################################ # convert to different orientation representations (numpy arrays) @@ -283,9 +158,10 @@ class Rotation: """Rotation matrix""" return qu2om(self.quaternion.asArray()) - def asRodrigues(self): + def asRodrigues(self,vector=False): """Rodrigues-Frank vector: ([n_1, n_2, n_3], tan(ω/2))""" - return qu2ro(self.quaternion.asArray()) + ro = qu2ro(self.quaternion.asArray()) + return ro[:3]*ro[3] if vector else ro def asHomochoric(self): """Homochoric vector: (h_1, h_2, h_3)""" @@ -294,6 +170,10 @@ class Rotation: def asCubochoric(self): return qu2cu(self.quaternion.asArray()) + def asM(self): + """Intermediate representation supporting quaternion averaging (see F. Landis Markley et al.)""" + return self.quaternion.asM() + ################################################################################################ # static constructors. The input data needs to follow the convention, options allow to @@ -402,63 +282,47 @@ class Rotation: return cls(ho2qu(ho)) - def __mul__(self, other): + @classmethod + def fromAverage(cls, + rotations, + weights = []): """ - Multiplication - - Rotation: Details needed (active/passive), rotation of (3,3,3,3)-matrix should be considered + Average rotation + + 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: input a list of rotations and optional weights """ - if isinstance(other, Rotation): # rotate a rotation - return self.__class__((self.quaternion * other.quaternion).asArray()) - elif isinstance(other, np.ndarray): - if other.shape == (3,): # rotate a single (3)-vector - ( x, y, z) = self.quaternion.p - (Vx,Vy,Vz) = other[0:3] - A = self.quaternion.q*self.quaternion.q - np.dot(self.quaternion.p,self.quaternion.p) - B = 2.0 * (x*Vx + y*Vy + z*Vz) - C = 2.0 * P*self.quaternion.q + if not all(isinstance(item, Rotation) for item in rotations): + raise TypeError("Only instances of Rotation can be averaged.") - 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), - ]) - elif other.shape == (3,3,): # rotate a single (3x3)-matrix - return np.dot(self.asMatrix(),np.dot(other,self.asMatrix().T)) - elif other.shape == (3,3,3,3): - raise NotImplementedError - else: - return NotImplemented - elif isinstance(other, tuple): # used to rotate a meshgrid-tuple - ( x, y, z) = self.quaternion.p - (Vx,Vy,Vz) = other[0:3] - A = self.quaternion.q*self.quaternion.q - np.dot(self.quaternion.p,self.quaternion.p) - B = 2.0 * (x*Vx + y*Vy + z*Vz) - C = 2.0 * P*self.quaternion.q + N = len(rotations) + if weights == [] or not weights: + weights = np.ones(N,dtype='i') - 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), - ]) - else: - return NotImplemented - - - def inverse(self): - """Inverse rotation/backward rotation""" - self.quaternion.conjugate() - return self - - def inversed(self): - """In-place inverse rotation/backward rotation""" - return self.__class__(self.quaternion.conjugated()) + for i,(r,n) in enumerate(zip(rotations,weights)): + M = r.asM() * n if i == 0 \ + else M + r.asM() * n # noqa add (multiples) of this rotation to average noqa + eig, vec = np.linalg.eig(M/N) + + return cls.fromQuaternion(np.real(vec.T[eig.argmax()]),acceptHomomorph = True) + + + @classmethod + def fromRandom(cls): + r = np.random.random(3) + A = np.sqrt(r[2]) + B = np.sqrt(1.0-r[2]) + w = np.cos(2.0*np.pi*r[0])*A + x = np.sin(2.0*np.pi*r[1])*B + y = np.cos(2.0*np.pi*r[1])*B + z = np.sin(2.0*np.pi*r[0])*A + return cls.fromQuaternion([w,x,y,z],acceptHomomorph=True) - def misorientation(self,other): - """Misorientation""" - return self.__class__(other.quaternion*self.quaternion.conjugated()) - # ****************************************************************************************** class Symmetry: @@ -503,60 +367,60 @@ class Symmetry: otherOrder = Symmetry.lattices.index(other.lattice) return (myOrder > otherOrder) - (myOrder < otherOrder) - def symmetryOperations(self): - """List of symmetry operations as quaternions.""" + def symmetryOperations(self,members=[]): + """List (or single element) of symmetry operations as rotations.""" if self.lattice == 'cubic': symQuats = [ - [ 1.0, 0.0, 0.0, 0.0 ], - [ 0.0, 1.0, 0.0, 0.0 ], - [ 0.0, 0.0, 1.0, 0.0 ], - [ 0.0, 0.0, 0.0, 1.0 ], - [ 0.0, 0.0, 0.5*math.sqrt(2), 0.5*math.sqrt(2) ], - [ 0.0, 0.0, 0.5*math.sqrt(2),-0.5*math.sqrt(2) ], - [ 0.0, 0.5*math.sqrt(2), 0.0, 0.5*math.sqrt(2) ], - [ 0.0, 0.5*math.sqrt(2), 0.0, -0.5*math.sqrt(2) ], - [ 0.0, 0.5*math.sqrt(2),-0.5*math.sqrt(2), 0.0 ], - [ 0.0, -0.5*math.sqrt(2),-0.5*math.sqrt(2), 0.0 ], - [ 0.5, 0.5, 0.5, 0.5 ], - [-0.5, 0.5, 0.5, 0.5 ], - [-0.5, 0.5, 0.5, -0.5 ], - [-0.5, 0.5, -0.5, 0.5 ], - [-0.5, -0.5, 0.5, 0.5 ], - [-0.5, -0.5, 0.5, -0.5 ], - [-0.5, -0.5, -0.5, 0.5 ], - [-0.5, 0.5, -0.5, -0.5 ], - [-0.5*math.sqrt(2), 0.0, 0.0, 0.5*math.sqrt(2) ], - [ 0.5*math.sqrt(2), 0.0, 0.0, 0.5*math.sqrt(2) ], - [-0.5*math.sqrt(2), 0.0, 0.5*math.sqrt(2), 0.0 ], - [-0.5*math.sqrt(2), 0.0, -0.5*math.sqrt(2), 0.0 ], - [-0.5*math.sqrt(2), 0.5*math.sqrt(2), 0.0, 0.0 ], - [-0.5*math.sqrt(2),-0.5*math.sqrt(2), 0.0, 0.0 ], + [ 1.0, 0.0, 0.0, 0.0 ], + [ 0.0, 1.0, 0.0, 0.0 ], + [ 0.0, 0.0, 1.0, 0.0 ], + [ 0.0, 0.0, 0.0, 1.0 ], + [ 0.0, 0.0, 0.5*np.sqrt(2), 0.5*np.sqrt(2) ], + [ 0.0, 0.0, 0.5*np.sqrt(2),-0.5*np.sqrt(2) ], + [ 0.0, 0.5*np.sqrt(2), 0.0, 0.5*np.sqrt(2) ], + [ 0.0, 0.5*np.sqrt(2), 0.0, -0.5*np.sqrt(2) ], + [ 0.0, 0.5*np.sqrt(2),-0.5*np.sqrt(2), 0.0 ], + [ 0.0, -0.5*np.sqrt(2),-0.5*np.sqrt(2), 0.0 ], + [ 0.5, 0.5, 0.5, 0.5 ], + [-0.5, 0.5, 0.5, 0.5 ], + [-0.5, 0.5, 0.5, -0.5 ], + [-0.5, 0.5, -0.5, 0.5 ], + [-0.5, -0.5, 0.5, 0.5 ], + [-0.5, -0.5, 0.5, -0.5 ], + [-0.5, -0.5, -0.5, 0.5 ], + [-0.5, 0.5, -0.5, -0.5 ], + [-0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ], + [ 0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ], + [-0.5*np.sqrt(2), 0.0, 0.5*np.sqrt(2), 0.0 ], + [-0.5*np.sqrt(2), 0.0, -0.5*np.sqrt(2), 0.0 ], + [-0.5*np.sqrt(2), 0.5*np.sqrt(2), 0.0, 0.0 ], + [-0.5*np.sqrt(2),-0.5*np.sqrt(2), 0.0, 0.0 ], ] elif self.lattice == 'hexagonal': symQuats = [ - [ 1.0,0.0,0.0,0.0 ], - [-0.5*math.sqrt(3), 0.0, 0.0,-0.5 ], - [ 0.5, 0.0, 0.0, 0.5*math.sqrt(3) ], - [ 0.0,0.0,0.0,1.0 ], - [-0.5, 0.0, 0.0, 0.5*math.sqrt(3) ], - [-0.5*math.sqrt(3), 0.0, 0.0, 0.5 ], - [ 0.0,1.0,0.0,0.0 ], - [ 0.0,-0.5*math.sqrt(3), 0.5, 0.0 ], - [ 0.0, 0.5,-0.5*math.sqrt(3), 0.0 ], - [ 0.0,0.0,1.0,0.0 ], - [ 0.0,-0.5,-0.5*math.sqrt(3), 0.0 ], - [ 0.0, 0.5*math.sqrt(3), 0.5, 0.0 ], + [ 1.0, 0.0, 0.0, 0.0 ], + [-0.5*np.sqrt(3), 0.0, 0.0, -0.5 ], + [ 0.5, 0.0, 0.0, 0.5*np.sqrt(3) ], + [ 0.0, 0.0, 0.0, 1.0 ], + [-0.5, 0.0, 0.0, 0.5*np.sqrt(3) ], + [-0.5*np.sqrt(3), 0.0, 0.0, 0.5 ], + [ 0.0, 1.0, 0.0, 0.0 ], + [ 0.0, -0.5*np.sqrt(3), 0.5, 0.0 ], + [ 0.0, 0.5, -0.5*np.sqrt(3), 0.0 ], + [ 0.0, 0.0, 1.0, 0.0 ], + [ 0.0, -0.5, -0.5*np.sqrt(3), 0.0 ], + [ 0.0, 0.5*np.sqrt(3), 0.5, 0.0 ], ] elif self.lattice == 'tetragonal': symQuats = [ - [ 1.0,0.0,0.0,0.0 ], - [ 0.0,1.0,0.0,0.0 ], - [ 0.0,0.0,1.0,0.0 ], - [ 0.0,0.0,0.0,1.0 ], - [ 0.0, 0.5*math.sqrt(2), 0.5*math.sqrt(2), 0.0 ], - [ 0.0,-0.5*math.sqrt(2), 0.5*math.sqrt(2), 0.0 ], - [ 0.5*math.sqrt(2), 0.0, 0.0, 0.5*math.sqrt(2) ], - [-0.5*math.sqrt(2), 0.0, 0.0, 0.5*math.sqrt(2) ], + [ 1.0, 0.0, 0.0, 0.0 ], + [ 0.0, 1.0, 0.0, 0.0 ], + [ 0.0, 0.0, 1.0, 0.0 ], + [ 0.0, 0.0, 0.0, 1.0 ], + [ 0.0, 0.5*np.sqrt(2), 0.5*np.sqrt(2), 0.0 ], + [ 0.0, -0.5*np.sqrt(2), 0.5*np.sqrt(2), 0.0 ], + [ 0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ], + [-0.5*np.sqrt(2), 0.0, 0.0, 0.5*np.sqrt(2) ], ] elif self.lattice == 'orthorhombic': symQuats = [ @@ -570,16 +434,28 @@ class Symmetry: [ 1.0,0.0,0.0,0.0 ], ] - return [Rotation(q) for q in symQuats] + 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 - def inFZ(self,R): + def inFZ(self,rodrigues): """ 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 (len(rodrigues) != 3): + raise ValueError('Input is not a Rodriques-Frank vector.\n') + + if np.any(rodrigues == np.inf): return False + + Rabs = abs(rodrigues) if self.lattice == 'cubic': return math.sqrt(2.0)-1.0 >= Rabs[0] \ @@ -609,14 +485,10 @@ class Symmetry: Representation of Orientation and Disorientation Data for Cubic, Hexagonal, Tetragonal and Orthorhombic Crystals Acta Cryst. (1991). A47, 780-789 """ - if isinstance(rodrigues, Quaternion): - R = rodrigues.asRodrigues() # translate accidentially passed quaternion - else: - R = rodrigues + if (len(rodrigues) != 3): + raise ValueError('Input is not a Rodriques-Frank vector.\n') + R = rodrigues - if R.shape[0]==4: # transition old (length not stored separately) to new - R = (R[0:3]*R[3]) - epsilon = 0.0 if self.lattice == 'cubic': return R[0] >= R[1]+epsilon and R[1] >= R[2]+epsilon and R[2] >= epsilon @@ -864,7 +736,7 @@ class Lattice: # Greninger--Troiano' orientation relationship for fcc <-> bcc transformation # from Y. He et al. Journal of Applied Crystallography 39 (2006) 72-81 - GTdash = {'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]], @@ -977,7 +849,7 @@ class Lattice: [[ 1, 1, 0],[ 1, 1, -1]]],dtype='float')} # Bain orientation relationship for fcc <-> bcc transformation - # from Y. He et al. Journal of Applied Crystallography 39 (2006) 72-81 + # 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]], @@ -990,7 +862,7 @@ class Lattice: def relationOperations(self,model): - models={'KS':self.KS, 'GT':self.GT, "GT'":self.GTdash, + models={'KS':self.KS, 'GT':self.GT, "GT'":self.GTprime, 'NW':self.NW, 'Pitsch': self.Pitsch, 'Bain':self.Bain} try: relationship = models[model] @@ -1050,7 +922,8 @@ class Orientation: def disorientation(self, other, - SST = True): + SST = True, + symmetries = False): """ Disorientation between myself and given other orientation. @@ -1058,33 +931,42 @@ class Orientation: (Currently requires same symmetry for both orientations. Look into A. Heinz and P. Neumann 1991 for cases with differing sym.) """ - #if self.lattice.symmetry != other.lattice.symmetry: - # raise NotImplementedError('disorientation between different symmetry classes not supported yet.') + if self.lattice.symmetry != other.lattice.symmetry: + 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([0]) # take all or only first sym operation otherSymEqs = other.equivalentOrientations() for i,sA in enumerate(mySymEqs): + aInv = sA.rotation.inversed() for j,sB in enumerate(otherSymEqs): - r = sB.rotation*mis*sA.rotation.inversed() + b = sB.rotation + r = b*aInv for k in range(2): - r.inversed() - breaker = self.lattice.symmetry.inFZ(r.asRodrigues()) \ - and (not SST or other.lattice.symmetry.inDisorientationSST(r.asRodrigues())) + r.inverse() + breaker = self.lattice.symmetry.inFZ(r.asRodrigues(vector=True)) \ + and (not SST or other.lattice.symmetry.inDisorientationSST(r.asRodrigues(vector=True))) if breaker: break if breaker: break if breaker: break - return r + 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.asRodrigues()) + return self.lattice.symmetry.inFZ(self.rotation.asRodrigues(vector=True)) - def equivalentOrientations(self): + def equivalentOrientations(self,members=[]): """List of orientations which are symmetrically equivalent""" - return [self.__class__(q*self.rotation,self.lattice) \ - for q in self.lattice.symmetry.symmetryOperations()] + try: + iter(members) # asking for (even empty) list of members? + except TypeError: + return self.__class__(self.lattice.symmetry.symmetryOperations(members)*self.rotation,self.lattice) # no, return rotation object + else: + return [self.__class__(q*self.rotation,self.lattice) \ + for q in self.lattice.symmetry.symmetryOperations(members)] # yes, return list of rotations def relatedOrientations(self,model): """List of orientations related by the given orientation relationship""" @@ -1094,7 +976,7 @@ class Orientation: 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.asRodrigues()): break + if self.lattice.symmetry.inFZ(me.rotation.asRodrigues(vector=True)): break return self.__class__(me.rotation,self.lattice) @@ -1125,37 +1007,28 @@ class Orientation: return color - # @classmethod - # def average(cls, - # orientations, - # multiplicity = []): - # """ - # Average orientation + @classmethod + def fromAverage(cls, + orientations, + weights = []): + """Create orientation from average of list of orientations""" + if not all(isinstance(item, Orientation) for item in orientations): + raise TypeError("Only instances of Orientation can be averaged.") - # 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.") + closest = [] + ref = orientations[0] + for o in orientations: + closest.append(o.equivalentOrientations( + ref.disorientation(o, + SST = False, # select (o[ther]'s) sym orientation + symmetries = True)[2]).rotation) # with lowest misorientation - # N = len(orientations) - # if multiplicity == [] or not multiplicity: - # multiplicity = np.ones(N,dtype='i') + return Orientation(Rotation.fromAverage(closest,weights),ref.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) + def average(self,other): + """Calculate the average rotation""" + return Orientation.fromAverage([self,other]) #################################################################################################### @@ -1195,7 +1068,167 @@ def isone(a): def iszero(a): return np.isclose(a,0.0,atol=1.0e-12,rtol=0.0) +#---------- quaternion ---------- + +def qu2om(qu): + """Quaternion to orientation matrix""" + qq = qu[0]**2-(qu[1]**2 + qu[2]**2 + qu[3]**2) + om = np.diag(qq + 2.0*np.array([qu[1],qu[2],qu[3]])**2) + + om[1,0] = 2.0*(qu[2]*qu[1]+qu[0]*qu[3]) + om[0,1] = 2.0*(qu[1]*qu[2]-qu[0]*qu[3]) + om[2,1] = 2.0*(qu[3]*qu[2]+qu[0]*qu[1]) + om[1,2] = 2.0*(qu[2]*qu[3]-qu[0]*qu[1]) + om[0,2] = 2.0*(qu[1]*qu[3]+qu[0]*qu[2]) + om[2,0] = 2.0*(qu[3]*qu[1]-qu[0]*qu[2]) + return om if P > 0.0 else om.T + +def qu2eu(qu): + """Quaternion to Euler angles""" + q03 = qu[0]**2+qu[3]**2 + q12 = qu[1]**2+qu[2]**2 + chi = np.sqrt(q03*q12) + + if iszero(chi): + 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: + 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 )]) + + # reduce Euler angles to definition range, i.e a lower limit of 0.0 + eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu) + return eu + + +def qu2ax(qu): + """ + Quaternion to axis angle + + Modified version of the original formulation, should be numerically more stable + """ + if iszero(qu[1]**2+qu[2]**2+qu[3]**2): # set axis to [001] if the angle is 0/360 + ax = [ 0.0, 0.0, 1.0, 0.0 ] + elif not iszero(qu[0]): + s = np.sign(qu[0])/np.sqrt(qu[1]**2+qu[2]**2+qu[3]**2) + omega = 2.0 * np.arccos(np.clip(qu[0],-1.0,1.0)) + ax = [ qu[1]*s, qu[2]*s, qu[3]*s, omega ] + else: + ax = [ qu[1], qu[2], qu[3], np.pi] + + return np.array(ax) + + +def qu2ro(qu): + """Quaternion to Rodrigues vector""" + if iszero(qu[0]): + ro = [qu[1], qu[2], qu[3], np.inf] + else: + s = np.linalg.norm([qu[1],qu[2],qu[3]]) + ro = [0.0,0.0,P,0.0] if iszero(s) else \ + [ qu[1]/s, qu[2]/s, qu[3]/s, np.tan(np.arccos(np.clip(qu[0],-1.0,1.0)))] # avoid numerical difficulties + + return np.array(ro) + + +def qu2ho(qu): + """Quaternion to homochoric""" + omega = 2.0 * np.arccos(np.clip(qu[0],-1.0,1.0)) # avoid numerical difficulties + + if iszero(omega): + ho = np.array([ 0.0, 0.0, 0.0 ]) + else: + ho = np.array([qu[1], qu[2], qu[3]]) + f = 0.75 * ( omega - np.sin(omega) ) + ho = ho/np.linalg.norm(ho) * f**(1./3.) + + return ho + + +def qu2cu(qu): + """Quaternion to cubochoric""" + return ho2cu(qu2ho(qu)) + + +#---------- orientation matrix ---------- + +def om2qu(om): + """ + Orientation matrix to quaternion + + The original formulation (direct conversion) had (numerical?) issues + """ + return eu2qu(om2eu(om)) + + +def om2eu(om): + """Orientation matrix to Euler angles""" + if abs(om[2,2]) < 1.0: + zeta = 1.0/np.sqrt(1.0-om[2,2]**2) + eu = np.array([np.arctan2(om[2,0]*zeta,-om[2,1]*zeta), + np.arccos(om[2,2]), + np.arctan2(om[0,2]*zeta, om[1,2]*zeta)]) + else: + eu = np.array([np.arctan2( om[0,1],om[0,0]), np.pi*0.5*(1-om[2,2]),0.0]) # following the paper, not the reference implementation + + # reduce Euler angles to definition range, i.e a lower limit of 0.0 + eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu) + return eu + + +def om2ax(om): + """Orientation matrix to axis angle""" + ax=np.empty(4) + + # first get the rotation angle + t = 0.5*(om.trace() -1.0) + ax[3] = np.arccos(np.clip(t,-1.0,1.0)) + + if iszero(ax[3]): + ax = [ 0.0, 0.0, 1.0, 0.0] + else: + w,vr = np.linalg.eig(om) + # next, find the eigenvalue (1,0j) + i = np.where(np.isclose(w,1.0+0.0j))[0][0] + ax[0:3] = np.real(vr[0:3,i]) + diagDelta = np.array([om[1,2]-om[2,1],om[2,0]-om[0,2],om[0,1]-om[1,0]]) + ax[0:3] = np.where(iszero(diagDelta), ax[0:3],np.abs(ax[0:3])*np.sign(-P*diagDelta)) + + return np.array(ax) + + +def om2ro(om): + """Orientation matrix to Rodriques vector""" + return eu2ro(om2eu(om)) + + +def om2ho(om): + """Orientation matrix to homochoric""" + return ax2ho(om2ax(om)) + + +def om2cu(om): + """Orientation matrix to cubochoric""" + return ho2cu(om2ho(om)) + + +#---------- Euler angles ---------- + +def eu2qu(eu): + """Euler angles to quaternion""" + ee = 0.5*eu + cPhi = np.cos(ee[1]) + sPhi = np.sin(ee[1]) + qu = np.array([ cPhi*np.cos(ee[0]+ee[2]), + -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*=-1 + return qu + + def eu2om(eu): """Euler angles to orientation matrix""" c = np.cos(eu) @@ -1241,32 +1274,28 @@ def eu2ro(eu): return ro -def eu2qu(eu): - """Euler angles to quaternion""" - ee = 0.5*eu - cPhi = np.cos(ee[1]) - sPhi = np.sin(ee[1]) - qu = np.array([ cPhi*np.cos(ee[0]+ee[2]), - -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 - return qu +def eu2ho(eu): + """Euler angles to homochoric""" + return ax2ho(eu2ax(eu)) -def om2eu(om): - """Euler angles to orientation matrix""" - if isone(om[2,2]**2): - eu = np.array([np.arctan2( om[0,1],om[0,0]), np.pi*0.5*(1-om[2,2]),0.0]) # following the paper, not the reference implementation +def eu2cu(eu): + """Euler angles to cubochoric""" + return ho2cu(eu2ho(eu)) + + +#---------- axis angle ---------- + +def ax2qu(ax): + """Axis angle to quaternion""" + if iszero(ax[3]): + qu = np.array([ 1.0, 0.0, 0.0, 0.0 ]) else: - zeta = 1.0/np.sqrt(1.0-om[2,2]**2) - eu = np.array([np.arctan2(om[2,0]*zeta,-om[2,1]*zeta), - np.arccos(om[2,2]), - np.arctan2(om[0,2]*zeta, om[1,2]*zeta)]) + c = np.cos(ax[3]*0.5) + s = np.sin(ax[3]*0.5) + qu = np.array([ c, ax[0]*s, ax[1]*s, ax[2]*s ]) - # reduce Euler angles to definition range, i.e a lower limit of 0.0 - eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu) - return eu + return qu def ax2om(ax): @@ -1283,24 +1312,23 @@ def ax2om(ax): return om if P < 0.0 else om.T + +def ax2eu(ax): + """Orientation matrix to Euler angles""" + return om2eu(ax2om(ax)) -def qu2eu(qu): - """Quaternion to Euler angles""" - q03 = qu[0]**2+qu[3]**2 - q12 = qu[1]**2+qu[2]**2 - chi = np.sqrt(q03*q12) - - if iszero(chi): - 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]) + +def ax2ro(ax): + """Axis angle to Rodrigues vector""" + if iszero(ax[3]): + ro = [ 0.0, 0.0, P, 0.0 ] else: - 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 )]) - - # reduce Euler angles to definition range, i.e a lower limit of 0.0 - eu = np.where(eu<0, (eu+2.0*np.pi)%np.array([2.0*np.pi,np.pi,2.0*np.pi]),eu) - return eu + ro = [ax[0], ax[1], ax[2]] + # 180 degree case + ro += [np.inf] if np.isclose(ax[3],np.pi,atol=1.0e-15,rtol=0.0) else \ + [np.tan(ax[3]*0.5)] + + return np.array(ro) def ax2ho(ax): @@ -1310,6 +1338,77 @@ def ax2ho(ax): return ho +def ax2cu(ax): + """Axis angle to cubochoric""" + return ho2cu(ax2ho(ax)) + + +#---------- Rodrigues--Frank ---------- + +def ro2qu(ro): + """Rodrigues vector to quaternion""" + return ax2qu(ro2ax(ro)) + + +def ro2om(ro): + """Rodgrigues vector to orientation matrix""" + return ax2om(ro2ax(ro)) + + +def ro2eu(ro): + """Rodrigues vector to orientation matrix""" + return om2eu(ro2om(ro)) + + +def ro2ax(ro): + """Rodrigues vector to axis angle""" + ta = ro[3] + + if iszero(ta): + ax = [ 0.0, 0.0, 1.0, 0.0 ] + elif not np.isfinite(ta): + ax = [ ro[0], ro[1], ro[2], np.pi ] + else: + angle = 2.0*np.arctan(ta) + ta = 1.0/np.linalg.norm(ro[0:3]) + ax = [ ro[0]/ta, ro[1]/ta, ro[2]/ta, angle ] + + return np.array(ax) + + +def ro2ho(ro): + """Rodrigues vector to homochoric""" + if iszero(np.sum(ro[0:3]**2.0)): + ho = [ 0.0, 0.0, 0.0 ] + else: + f = 2.0*np.arctan(ro[3]) -np.sin(2.0*np.arctan(ro[3])) if np.isfinite(ro[3]) else np.pi + ho = ro[0:3] * (0.75*f)**(1.0/3.0) + + return np.array(ho) + + +def ro2cu(ro): + """Rodrigues vector to cubochoric""" + return ho2cu(ro2ho(ro)) + + +#---------- homochoric ---------- + +def ho2qu(ho): + """Homochoric to quaternion""" + return ax2qu(ho2ax(ho)) + + +def ho2om(ho): + """Homochoric to orientation matrix""" + return ax2om(ho2ax(ho)) + + +def ho2eu(ho): + """Homochoric to Euler angles""" + return ax2eu(ho2ax(ho)) + + def ho2ax(ho): """Homochoric to axis angle""" tfit = np.array([+1.0000000000018852, -0.5000000002194847, @@ -1336,135 +1435,9 @@ def ho2ax(ho): return ax -def om2ax(om): - """Orientation matrix to axis angle""" - ax=np.empty(4) - - # first get the rotation angle - t = 0.5*(om.trace() -1.0) - ax[3] = np.arccos(np.clip(t,-1.0,1.0)) - - if iszero(ax[3]): - ax = [ 0.0, 0.0, 1.0, 0.0] - else: - w,vr = np.linalg.eig(om) - # next, find the eigenvalue (1,0j) - i = np.where(np.isclose(w,1.0+0.0j))[0][0] - ax[0:3] = np.real(vr[0:3,i]) - diagDelta = np.array([om[1,2]-om[2,1],om[2,0]-om[0,2],om[0,1]-om[1,0]]) - ax[0:3] = np.where(iszero(diagDelta), ax[0:3],np.abs(ax[0:3])*np.sign(-P*diagDelta)) - - return np.array(ax) - - -def ro2ax(ro): - """Rodrigues vector to axis angle""" - ta = ro[3] - - if iszero(ta): - ax = [ 0.0, 0.0, 1.0, 0.0 ] - elif not np.isfinite(ta): - ax = [ ro[0], ro[1], ro[2], np.pi ] - else: - angle = 2.0*np.arctan(ta) - ta = 1.0/np.linalg.norm(ro[0:3]) - ax = [ ro[0]/ta, ro[1]/ta, ro[2]/ta, angle ] - - return np.array(ax) - - -def ax2ro(ax): - """Axis angle to Rodrigues vector""" - if iszero(ax[3]): - ro = [ 0.0, 0.0, P, 0.0 ] - else: - ro = [ax[0], ax[1], ax[2]] - # 180 degree case - ro += [np.inf] if np.isclose(ax[3],np.pi,atol=1.0e-15,rtol=0.0) else \ - [np.tan(ax[3]*0.5)] - - return np.array(ro) - - -def ax2qu(ax): - """Axis angle to quaternion""" - if iszero(ax[3]): - qu = np.array([ 1.0, 0.0, 0.0, 0.0 ]) - else: - c = np.cos(ax[3]*0.5) - s = np.sin(ax[3]*0.5) - qu = np.array([ c, ax[0]*s, ax[1]*s, ax[2]*s ]) - - return qu - - -def ro2ho(ro): - """Rodrigues vector to homochoric""" - if iszero(np.sum(ro[0:3]**2.0)): - ho = [ 0.0, 0.0, 0.0 ] - else: - f = 2.0*np.arctan(ro[3]) -np.sin(2.0*np.arctan(ro[3])) if np.isfinite(ro[3]) else np.pi - ho = ro[0:3] * (0.75*f)**(1.0/3.0) - - return np.array(ho) - - -def qu2om(qu): - """Quaternion to orientation matrix""" - qq = qu[0]**2-(qu[1]**2 + qu[2]**2 + qu[3]**2) - om = np.diag(qq + 2.0*np.array([qu[1],qu[2],qu[3]])**2) - - om[1,0] = 2.0*(qu[2]*qu[1]+qu[0]*qu[3]) - om[0,1] = 2.0*(qu[1]*qu[2]-qu[0]*qu[3]) - om[2,1] = 2.0*(qu[3]*qu[2]+qu[0]*qu[1]) - om[1,2] = 2.0*(qu[2]*qu[3]-qu[0]*qu[1]) - om[0,2] = 2.0*(qu[1]*qu[3]+qu[0]*qu[2]) - om[2,0] = 2.0*(qu[3]*qu[1]-qu[0]*qu[2]) - return om if P > 0.0 else om.T - - -def qu2ax(qu): - """ - Quaternion to axis angle - - Modified version of the original formulation, should be numerically more stable - """ - if isone(abs(qu[0])): # set axis to [001] if the angle is 0/360 - ax = [ 0.0, 0.0, 1.0, 0.0 ] - elif not iszero(qu[0]): - omega = 2.0 * np.arccos(qu[0]) - s = np.sign(qu[0])/np.sqrt(qu[1]**2+qu[2]**2+qu[3]**2) - ax = [ qu[1]*s, qu[2]*s, qu[3]*s, omega ] - else: - ax = [ qu[1], qu[2], qu[3], np.pi] - - return np.array(ax) - - -def qu2ro(qu): - """Quaternion to Rodrigues vector""" - if iszero(qu[0]): - ro = [qu[1], qu[2], qu[3], np.inf] - else: - s = np.linalg.norm([qu[1],qu[2],qu[3]]) - ro = [0.0,0.0,P,0.0] if iszero(s) else \ - [ qu[1]/s, qu[2]/s, qu[3]/s, np.tan(np.arccos(np.clip(qu[0],-1.0,1.0)))] # avoid numerical difficulties - - return np.array(ro) - - -def qu2ho(qu): - """Quaternion to homochoric""" - omega = 2.0 * np.arccos(np.clip(qu[0],-1.0,1.0)) # avoid numerical difficulties - - if iszero(omega): - ho = np.array([ 0.0, 0.0, 0.0 ]) - else: - ho = np.array([qu[1], qu[2], qu[3]]) - f = 0.75 * ( omega - np.sin(omega) ) - ho = ho/np.linalg.norm(ho) * f**(1./3.) - - return ho +def ho2ro(ho): + """Axis angle to Rodriques vector""" + return ax2ro(ho2ax(ho)) def ho2cu(ho): @@ -1472,103 +1445,11 @@ def ho2cu(ho): return Lambert.BallToCube(ho) -def cu2ho(cu): - """Cubochoric to homochoric""" - return Lambert.CubeToBall(cu) +#---------- cubochoric ---------- - -def ro2eu(ro): - """Rodrigues vector to orientation matrix""" - return om2eu(ro2om(ro)) - - -def eu2ho(eu): - """Euler angles to homochoric""" - return ax2ho(eu2ax(eu)) - - -def om2ro(om): - """Orientation matrix to Rodriques vector""" - return eu2ro(om2eu(om)) - - -def om2ho(om): - """Orientation matrix to homochoric""" - return ax2ho(om2ax(om)) - - -def ax2eu(ax): - """Orientation matrix to Euler angles""" - return om2eu(ax2om(ax)) - - -def ro2om(ro): - """Rodgrigues vector to orientation matrix""" - return ax2om(ro2ax(ro)) - - -def ro2qu(ro): - """Rodrigues vector to quaternion""" - return ax2qu(ro2ax(ro)) - - -def ho2eu(ho): - """Homochoric to Euler angles""" - return ax2eu(ho2ax(ho)) - - -def ho2om(ho): - """Homochoric to orientation matrix""" - return ax2om(ho2ax(ho)) - - -def ho2ro(ho): - """Axis angle to Rodriques vector""" - return ax2ro(ho2ax(ho)) - - -def ho2qu(ho): - """Homochoric to quaternion""" - return ax2qu(ho2ax(ho)) - - -def eu2cu(eu): - """Euler angles to cubochoric""" - return ho2cu(eu2ho(eu)) - - -def om2cu(om): - """Orientation matrix to cubochoric""" - return ho2cu(om2ho(om)) - - -def om2qu(om): - """ - Orientation matrix to quaternion - - The original formulation (direct conversion) had numerical issues - """ - return ax2qu(om2ax(om)) - - -def ax2cu(ax): - """Axis angle to cubochoric""" - return ho2cu(ax2ho(ax)) - - -def ro2cu(ro): - """Rodrigues vector to cubochoric""" - return ho2cu(ro2ho(ro)) - - -def qu2cu(qu): - """Quaternion to cubochoric""" - return ho2cu(qu2ho(qu)) - - -def cu2eu(cu): - """Cubochoric to Euler angles""" - return ho2eu(cu2ho(cu)) +def cu2qu(cu): + """Cubochoric to quaternion""" + return ho2qu(cu2ho(cu)) def cu2om(cu): @@ -1576,6 +1457,11 @@ def cu2om(cu): return ho2om(cu2ho(cu)) +def cu2eu(cu): + """Cubochoric to Euler angles""" + return ho2eu(cu2ho(cu)) + + def cu2ax(cu): """Cubochoric to axis angle""" return ho2ax(cu2ho(cu)) @@ -1586,6 +1472,6 @@ def cu2ro(cu): return ho2ro(cu2ho(cu)) -def cu2qu(cu): - """Cubochoric to quaternion""" - return ho2qu(cu2ho(cu)) +def cu2ho(cu): + """Cubochoric to homochoric""" + return Lambert.CubeToBall(cu) diff --git a/python/damask/quaternion.py b/python/damask/quaternion.py new file mode 100644 index 000000000..fd681fe9b --- /dev/null +++ b/python/damask/quaternion.py @@ -0,0 +1,214 @@ +# -*- coding: UTF-8 no BOM -*- + +import numpy as np + +P = -1 # convention (see DOI:10.1088/0965-0393/23/8/083501) + +#################################################################################################### +class Quaternion: + u""" + Quaternion with basic operations + + q is the real part, p = (x, y, z) are the imaginary parts. + Definition of multiplication depends on variable P, P ∈ {-1,1}. + """ + + def __init__(self, + q = 0.0, + p = np.zeros(3,dtype=float)): + """Initializes to identity unless specified""" + self.q = float(q) + self.p = np.array(p,dtype=float) + + + def __copy__(self): + """Copy""" + return self.__class__(q=self.q, + p=self.p.copy()) + + copy = __copy__ + + + def __iter__(self): + """Components""" + return iter(self.asList()) + + 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 __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 + else: + return NotImplemented + + def __pos__(self): + """Unary positive operator""" + 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 + else: + return NotImplemented + + def __neg__(self): + """Unary negative operator""" + return self * -1.0 + + + def __mul__(self, other): + """Multiplication with quaternion or scalar""" + 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)): + return self.__class__(q=self.q*other, + p=self.p*other) + else: + return NotImplemented + + def __imul__(self, other): + """In-place multiplication with quaternion or scalar""" + 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 + elif isinstance(other, (int, float)): + self.q *= other + self.p *= other + return self + else: + return NotImplemented + + + def __truediv__(self, other): + """Divsion with quaternion or scalar""" + if isinstance(other, Quaternion): + s = other.conjugate()/abs(other)**2. + return self.__class__(q=self.q * s, + p=self.p * s) + elif isinstance(other, (int, float)): + return self.__class__(q=self.q / other, + p=self.p / other) + else: + return NotImplemented + + def __itruediv__(self, other): + """In-place divsion with quaternion or scalar""" + if isinstance(other, Quaternion): + s = other.conjugate()/abs(other)**2. + self *= s + return self + elif isinstance(other, (int, float)): + self.q /= other + self.p /= other + return self + else: + return NotImplemented + + + def __pow__(self, exponent): + """Power""" + if isinstance(exponent, (int, float)): + omega = np.acos(self.q) + return self.__class__(q= np.cos(exponent*omega), + p=self.p * np.sin(exponent*omega)/np.sin(omega)) + else: + return NotImplemented + + def __ipow__(self, exponent): + """In-place power""" + if isinstance(exponent, (int, float)): + omega = np.acos(self.q) + self.q = np.cos(exponent*omega) + self.p *= np.sin(exponent*omega)/np.sin(omega) + return self + else: + return NotImplemented + + + def __abs__(self): + """Norm""" + return np.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 asM(self): + """Intermediate representation useful for quaternion averaging (see F. Landis Markley et al.)""" + return np.outer(self.asArray(),self.asArray()) + + def asArray(self): + """As numpy array""" + return np.array((self.q,self.p[0],self.p[1],self.p[2])) + + def asList(self): + """As list""" + return [self.q]+list(self.p) + + + def normalize(self): + """Normalizes in-place""" + d = self.magnitude() + if d > 0.0: self /= d + return self + + def normalized(self): + """Returns normalized copy""" + return self.copy().normalize() + + + def conjugate(self): + """Conjugates in-place""" + self.p = -self.p + return self + + def conjugated(self): + """Returns conjugated copy""" + return self.copy().conjugate() + + + def homomorph(self): + """Homomorphs in-place""" + self *= -1.0 + return self + + def homomorphed(self): + """Returns homomorphed copy""" + return self.copy().homomorph() diff --git a/src/plastic_disloUCLA.f90 b/src/plastic_disloUCLA.f90 index efab7bc0b..56d910011 100644 --- a/src/plastic_disloUCLA.f90 +++ b/src/plastic_disloUCLA.f90 @@ -121,15 +121,7 @@ subroutine plastic_disloUCLA_init() math_expand use IO, only: & IO_error - use material, only: & - phase_plasticity, & - phase_plasticityInstance, & - phase_Noutput, & - material_allocatePlasticState, & - PLASTICITY_DISLOUCLA_label, & - PLASTICITY_DISLOUCLA_ID, & - material_phase, & - plasticState + use material use config, only: & config_phase use lattice diff --git a/src/plastic_dislotwin.f90 b/src/plastic_dislotwin.f90 index e8aca7f89..7db1f5f7f 100644 --- a/src/plastic_dislotwin.f90 +++ b/src/plastic_dislotwin.f90 @@ -185,15 +185,7 @@ subroutine plastic_dislotwin_init PI use IO, only: & IO_error - use material, only: & - phase_plasticity, & - phase_plasticityInstance, & - phase_Noutput, & - material_allocatePlasticState, & - PLASTICITY_DISLOTWIN_label, & - PLASTICITY_DISLOTWIN_ID, & - material_phase, & - plasticState + use material use config, only: & config_phase use lattice diff --git a/src/plastic_isotropic.f90 b/src/plastic_isotropic.f90 index 4f2892f57..c572f0ded 100644 --- a/src/plastic_isotropic.f90 +++ b/src/plastic_isotropic.f90 @@ -91,19 +91,7 @@ subroutine plastic_isotropic_init debug_levelBasic use IO, only: & IO_error - use material, only: & - phase_plasticity, & - phase_plasticityInstance, & - phase_Noutput, & - material_allocatePlasticState, & - PLASTICITY_ISOTROPIC_label, & - PLASTICITY_ISOTROPIC_ID, & - material_phase, & - plasticState -#ifdef DEBUG - use material, only: & - phasememberAt -#endif + use material use config, only: & config_phase use lattice @@ -464,7 +452,7 @@ function plastic_isotropic_postResults(Mp,instance,of) result(postResults) c = c + 1 case (dot_gamma_ID) postResults(c+1) = prm%dot_gamma_0 & - * (sqrt(1.5_pReal) * norm_Mp /(prm%M * stt%xi(of)))**prm%n + * (sqrt(1.5_pReal) * norm_Mp /(prm%M * stt%xi(of)))**prm%n c = c + 1 end select @@ -479,11 +467,10 @@ end function plastic_isotropic_postResults !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- subroutine plastic_isotropic_results(instance,group) -#if defined(PETSc) || defined(DAMASK_HDF5) - use results, only: & - results_writeDataset +#if defined(PETSc) || defined(DAMASKHDF5) + use results - integer, intent(in) :: instance + integer, intent(in) :: instance character(len=*), intent(in) :: group integer :: o @@ -496,10 +483,9 @@ subroutine plastic_isotropic_results(instance,group) end select enddo outputsLoop end associate - #else - integer, intent(in) :: instance - character(len=*), intent(in) :: group + integer, intent(in) :: instance + character(len=*) :: group #endif end subroutine plastic_isotropic_results diff --git a/src/plastic_kinematichardening.f90 b/src/plastic_kinematichardening.f90 index a255572e1..861b98da3 100644 --- a/src/plastic_kinematichardening.f90 +++ b/src/plastic_kinematichardening.f90 @@ -112,19 +112,7 @@ subroutine plastic_kinehardening_init math_expand use IO, only: & IO_error - use material, only: & - phase_plasticity, & - phase_plasticityInstance, & - phase_Noutput, & - material_allocatePlasticState, & - PLASTICITY_kinehardening_label, & - PLASTICITY_kinehardening_ID, & - material_phase, & - plasticState -#ifdef DEBUG - use material, only: & - phasememberAt -#endif + use material use config, only: & config_phase use lattice diff --git a/src/plastic_none.f90 b/src/plastic_none.f90 index 6f2ef2230..4b14266f1 100644 --- a/src/plastic_none.f90 +++ b/src/plastic_none.f90 @@ -23,13 +23,7 @@ subroutine plastic_none_init debug_level, & debug_constitutive, & debug_levelBasic - use material, only: & - phase_plasticity, & - material_allocatePlasticState, & - PLASTICITY_NONE_label, & - PLASTICITY_NONE_ID, & - material_phase, & - plasticState + use material integer :: & Ninstance, & diff --git a/src/plastic_nonlocal.f90 b/src/plastic_nonlocal.f90 index caaa0e4a2..6097bbbc8 100644 --- a/src/plastic_nonlocal.f90 +++ b/src/plastic_nonlocal.f90 @@ -248,15 +248,7 @@ subroutine plastic_nonlocal_init debug_levelBasic use mesh, only: & theMesh - use material, only: & - phase_plasticity, & - phase_plasticityInstance, & - phase_Noutput, & - PLASTICITY_NONLOCAL_label, & - PLASTICITY_NONLOCAL_ID, & - plasticState, & - material_phase, & - material_allocatePlasticState + use material use config use lattice diff --git a/src/plastic_phenopowerlaw.f90 b/src/plastic_phenopowerlaw.f90 index f8ebae68d..196129f64 100644 --- a/src/plastic_phenopowerlaw.f90 +++ b/src/plastic_phenopowerlaw.f90 @@ -116,15 +116,7 @@ subroutine plastic_phenopowerlaw_init math_expand use IO, only: & IO_error - use material, only: & - phase_plasticity, & - phase_plasticityInstance, & - phase_Noutput, & - material_allocatePlasticState, & - PLASTICITY_PHENOPOWERLAW_LABEL, & - PLASTICITY_PHENOPOWERLAW_ID, & - material_phase, & - plasticState + use material use config, only: & config_phase use lattice diff --git a/src/rotations.f90 b/src/rotations.f90 index 7c29b03bc..69529ed24 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -77,7 +77,6 @@ contains !--------------------------------------------------------------------------------------------------- function asQuaternion(self) - implicit none class(rotation), intent(in) :: self real(pReal), dimension(4) :: asQuaternion @@ -87,7 +86,6 @@ end function asQuaternion !--------------------------------------------------------------------------------------------------- function asEulerAngles(self) - implicit none class(rotation), intent(in) :: self real(pReal), dimension(3) :: asEulerAngles @@ -97,7 +95,6 @@ end function asEulerAngles !--------------------------------------------------------------------------------------------------- function asAxisAnglePair(self) - implicit none class(rotation), intent(in) :: self real(pReal), dimension(4) :: asAxisAnglePair @@ -107,7 +104,6 @@ end function asAxisAnglePair !--------------------------------------------------------------------------------------------------- function asRotationMatrix(self) - implicit none class(rotation), intent(in) :: self real(pReal), dimension(3,3) :: asRotationMatrix @@ -117,7 +113,6 @@ end function asRotationMatrix !--------------------------------------------------------------------------------------------------- function asRodriguesFrankVector(self) - implicit none class(rotation), intent(in) :: self real(pReal), dimension(4) :: asRodriguesFrankVector @@ -127,7 +122,6 @@ end function asRodriguesFrankVector !--------------------------------------------------------------------------------------------------- function asHomochoric(self) - implicit none class(rotation), intent(in) :: self real(pReal), dimension(3) :: asHomochoric @@ -140,7 +134,6 @@ end function asHomochoric !--------------------------------------------------------------------------------------------------- subroutine fromRotationMatrix(self,om) - implicit none class(rotation), intent(out) :: self real(pReal), dimension(3,3), intent(in) :: om @@ -158,7 +151,6 @@ function rotVector(self,v,active) use prec, only: & dEq - implicit none real(pReal), dimension(3) :: rotVector class(rotation), intent(in) :: self real(pReal), intent(in), dimension(3) :: v @@ -198,7 +190,6 @@ end function rotVector !--------------------------------------------------------------------------------------------------- function rotTensor(self,m,active) - implicit none real(pReal), dimension(3,3) :: rotTensor class(rotation), intent(in) :: self real(pReal), intent(in), dimension(3,3) :: m @@ -226,7 +217,6 @@ end function rotTensor !--------------------------------------------------------------------------------------------------- function misorientation(self,other) - implicit none type(rotation) :: misorientation class(rotation), intent(in) :: self, other @@ -235,469 +225,12 @@ function misorientation(self,other) end function misorientation -!--------------------------------------------------------------------------------------------------- -! The following routines convert between different representations -!--------------------------------------------------------------------------------------------------- - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief Euler angles to orientation matrix -!--------------------------------------------------------------------------------------------------- -pure function eu2om(eu) result(om) - use prec, only: & - dEq0 - - implicit none - real(pReal), intent(in), dimension(3) :: eu - real(pReal), dimension(3,3) :: om - - real(pReal), dimension(3) :: c, s - - c = cos(eu) - s = sin(eu) - - om(1,1) = c(1)*c(3)-s(1)*s(3)*c(2) - om(1,2) = s(1)*c(3)+c(1)*s(3)*c(2) - om(1,3) = s(3)*s(2) - om(2,1) = -c(1)*s(3)-s(1)*c(3)*c(2) - om(2,2) = -s(1)*s(3)+c(1)*c(3)*c(2) - om(2,3) = c(3)*s(2) - om(3,1) = s(1)*s(2) - om(3,2) = -c(1)*s(2) - om(3,3) = c(2) - - where(dEq0(om)) om = 0.0_pReal - -end function eu2om - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert euler to axis angle -!--------------------------------------------------------------------------------------------------- -pure function eu2ax(eu) result(ax) - use prec, only: & - dEq0, & - dEq - use math, only: & - PI - - implicit none - real(pReal), intent(in), dimension(3) :: eu - real(pReal), dimension(4) :: ax - - real(pReal) :: t, delta, tau, alpha, sigma - - t = tan(eu(2)*0.5) - sigma = 0.5*(eu(1)+eu(3)) - delta = 0.5*(eu(1)-eu(3)) - tau = sqrt(t**2+sin(sigma)**2) - - alpha = merge(PI, 2.0*atan(tau/cos(sigma)), dEq(sigma,PI*0.5_pReal,tol=1.0e-15_pReal)) - - if (dEq0(alpha)) then ! return a default identity axis-angle pair - ax = [ 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal ] - else - ax(1:3) = -P/tau * [ t*cos(delta), t*sin(delta), sin(sigma) ] ! passive axis-angle pair so a minus sign in front - ax(4) = alpha - if (alpha < 0.0) ax = -ax ! ensure alpha is positive - end if - -end function eu2ax - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief Euler angles to Rodrigues vector -!--------------------------------------------------------------------------------------------------- -pure function eu2ro(eu) result(ro) - use prec, only: & - dEq0 - use, intrinsic :: IEEE_ARITHMETIC, only: & - IEEE_value, & - IEEE_positive_inf - use math, only: & - PI - - implicit none - real(pReal), intent(in), dimension(3) :: eu - real(pReal), dimension(4) :: ro - - ro = eu2ax(eu) - if (ro(4) >= PI) then - ro(4) = IEEE_value(ro(4),IEEE_positive_inf) - elseif(dEq0(ro(4))) then - ro = [ 0.0_pReal, 0.0_pReal, P, 0.0_pReal ] - else - ro(4) = tan(ro(4)*0.5) - end if - -end function eu2ro - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief Euler angles to unit quaternion -!--------------------------------------------------------------------------------------------------- -pure function eu2qu(eu) result(qu) - - implicit none - real(pReal), intent(in), dimension(3) :: eu - type(quaternion) :: qu - real(pReal), dimension(3) :: ee - real(pReal) :: cPhi, sPhi - - ee = 0.5_pReal*eu - - cPhi = cos(ee(2)) - sPhi = sin(ee(2)) - - qu = quaternion([ cPhi*cos(ee(1)+ee(3)), & - -P*sPhi*cos(ee(1)-ee(3)), & - -P*sPhi*sin(ee(1)-ee(3)), & - -P*cPhi*sin(ee(1)+ee(3))]) - if(qu%w < 0.0_pReal) qu = qu%homomorphed() - -end function eu2qu - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief orientation matrix to Euler angles -!--------------------------------------------------------------------------------------------------- -pure function om2eu(om) result(eu) - use math, only: & - PI - - implicit none - real(pReal), intent(in), dimension(3,3) :: om - real(pReal), dimension(3) :: eu - real(pReal) :: zeta - - if (abs(om(3,3))>1.0_pReal) then - eu = [ atan2( om(1,2),om(1,1)), 0.5*PI*(1-om(3,3)),0.0_pReal ] - else - zeta = 1.0_pReal/sqrt(1.0_pReal-om(3,3)**2.0_pReal) - eu = [atan2(om(3,1)*zeta,-om(3,2)*zeta), & - acos(om(3,3)), & - atan2(om(1,3)*zeta, om(2,3)*zeta)] - end if - where(eu<0.0_pReal) eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI]) - -end function om2eu - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert axis angle pair to orientation matrix -!--------------------------------------------------------------------------------------------------- -pure function ax2om(ax) result(om) - - implicit none - real(pReal), intent(in), dimension(4) :: ax - real(pReal), dimension(3,3) :: om - - real(pReal) :: q, c, s, omc - integer :: i - - c = cos(ax(4)) - s = sin(ax(4)) - omc = 1.0-c - - forall(i=1:3) om(i,i) = ax(i)**2*omc + c - - q = omc*ax(1)*ax(2) - om(1,2) = q + s*ax(3) - om(2,1) = q - s*ax(3) - - q = omc*ax(2)*ax(3) - om(2,3) = q + s*ax(1) - om(3,2) = q - s*ax(1) - - q = omc*ax(3)*ax(1) - om(3,1) = q + s*ax(2) - om(1,3) = q - s*ax(2) - - if (P > 0.0) om = transpose(om) - -end function ax2om - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert unit quaternion to Euler angles -!--------------------------------------------------------------------------------------------------- -pure function qu2eu(qu) result(eu) - use prec, only: & - dEq0 - use math, only: & - PI - - implicit none - type(quaternion), intent(in) :: qu - real(pReal), dimension(3) :: eu - - real(pReal) :: q12, q03, chi, chiInv - - q03 = qu%w**2+qu%z**2 - q12 = qu%x**2+qu%y**2 - chi = sqrt(q03*q12) - - degenerated: if (dEq0(chi)) then - eu = merge([atan2(-P*2.0*qu%w*qu%z,qu%w**2-qu%z**2), 0.0_pReal, 0.0_pReal], & - [atan2(2.0*qu%x*qu%y,qu%x**2-qu%y**2), PI, 0.0_pReal], & - dEq0(q12)) - else degenerated - chiInv = 1.0/chi - eu = [atan2((-P*qu%w*qu%y+qu%x*qu%z)*chi, (-P*qu%w*qu%x-qu%y*qu%z)*chi ), & - atan2( 2.0*chi, q03-q12 ), & - atan2(( P*qu%w*qu%y+qu%x*qu%z)*chi, (-P*qu%w*qu%x+qu%y*qu%z)*chi )] - endif degenerated - where(eu<0.0_pReal) eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI]) - -end function qu2eu - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert axis angle pair to homochoric -!--------------------------------------------------------------------------------------------------- -pure function ax2ho(ax) result(ho) - - implicit none - real(pReal), intent(in), dimension(4) :: ax - real(pReal), dimension(3) :: ho - - real(pReal) :: f - - f = 0.75 * ( ax(4) - sin(ax(4)) ) - f = f**(1.0/3.0) - ho = ax(1:3) * f - -end function ax2ho - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert homochoric to axis angle pair -!--------------------------------------------------------------------------------------------------- -pure function ho2ax(ho) result(ax) - use prec, only: & - dEq0 - - implicit none - real(pReal), intent(in), dimension(3) :: ho - real(pReal), dimension(4) :: ax - - integer :: i - real(pReal) :: hmag_squared, s, hm - real(pReal), parameter, dimension(16) :: & - tfit = [ 1.0000000000018852_pReal, -0.5000000002194847_pReal, & - -0.024999992127593126_pReal, -0.003928701544781374_pReal, & - -0.0008152701535450438_pReal, -0.0002009500426119712_pReal, & - -0.00002397986776071756_pReal, -0.00008202868926605841_pReal, & - +0.00012448715042090092_pReal, -0.0001749114214822577_pReal, & - +0.0001703481934140054_pReal, -0.00012062065004116828_pReal, & - +0.000059719705868660826_pReal, -0.00001980756723965647_pReal, & - +0.000003953714684212874_pReal, -0.00000036555001439719544_pReal ] - - ! normalize h and store the magnitude - hmag_squared = sum(ho**2.0_pReal) - if (dEq0(hmag_squared)) then - ax = [ 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal ] - else - hm = hmag_squared - - ! convert the magnitude to the rotation angle - s = tfit(1) + tfit(2) * hmag_squared - do i=3,16 - hm = hm*hmag_squared - s = s + tfit(i) * hm - end do - ax = [ho/sqrt(hmag_squared), 2.0_pReal*acos(s)] - end if - -end function ho2ax - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert orientation matrix to axis angle pair -!--------------------------------------------------------------------------------------------------- -function om2ax(om) result(ax) - use prec, only: & - dEq0, & - cEq, & - dNeq0 - use IO, only: & - IO_error - use math, only: & - math_clip, & - math_trace33 - - implicit none - real(pReal), intent(in) :: om(3,3) - real(pReal) :: ax(4) - - real(pReal) :: t - real(pReal), dimension(3) :: Wr, Wi - real(pReal), dimension(10) :: WORK - real(pReal), dimension(3,3) :: VR, devNull, o - integer :: INFO, LWORK, i - - external :: dgeev,sgeev - - o = om - - ! first get the rotation angle - t = 0.5_pReal * (math_trace33(om) - 1.0) - ax(4) = acos(math_clip(t,-1.0_pReal,1.0_pReal)) - - if (dEq0(ax(4))) then - ax(1:3) = [ 0.0, 0.0, 1.0 ] - else - ! set some initial LAPACK variables - INFO = 0 - ! first initialize the parameters for the LAPACK DGEEV routines - LWORK = 20 - - ! call the eigenvalue solver - call dgeev('N','V',3,o,3,Wr,Wi,devNull,3,VR,3,WORK,LWORK,INFO) - if (INFO /= 0) call IO_error(0,ext_msg='Error in om2ax DGEEV return not zero') - i = maxloc(merge(1.0_pReal,0.0_pReal,cEq(cmplx(Wr,Wi,pReal),cmplx(1.0_pReal,0.0_pReal,pReal),tol=1.0e-14_pReal)),dim=1) ! poor substitute for findloc - ax(1:3) = VR(1:3,i) - where ( dNeq0([om(2,3)-om(3,2), om(3,1)-om(1,3), om(1,2)-om(2,1)])) & - ax(1:3) = sign(ax(1:3),-P *[om(2,3)-om(3,2), om(3,1)-om(1,3), om(1,2)-om(2,1)]) - endif - -end function om2ax - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert Rodrigues vector to axis angle pair -!--------------------------------------------------------------------------------------------------- -pure function ro2ax(ro) result(ax) - use, intrinsic :: IEEE_ARITHMETIC, only: & - IEEE_is_finite - use prec, only: & - dEq0 - use math, only: & - PI - - implicit none - real(pReal), intent(in), dimension(4) :: ro - real(pReal), dimension(4) :: ax - - real(pReal) :: ta, angle - - ta = ro(4) - - if (dEq0(ta)) then - ax = [ 0.0, 0.0, 1.0, 0.0 ] - elseif (.not. IEEE_is_finite(ta)) then - ax = [ ro(1), ro(2), ro(3), PI ] - else - angle = 2.0*atan(ta) - ta = 1.0/norm2(ro(1:3)) - ax = [ ro(1)/ta, ro(2)/ta, ro(3)/ta, angle ] - end if - -end function ro2ax - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert axis angle pair to Rodrigues vector -!--------------------------------------------------------------------------------------------------- -pure function ax2ro(ax) result(ro) - use, intrinsic :: IEEE_ARITHMETIC, only: & - IEEE_value, & - IEEE_positive_inf - use prec, only: & - dEq0 - use math, only: & - PI - - implicit none - real(pReal), intent(in), dimension(4) :: ax - real(pReal), dimension(4) :: ro - - real(pReal), parameter :: thr = 1.0E-7 - - if (dEq0(ax(4))) then - ro = [ 0.0_pReal, 0.0_pReal, P, 0.0_pReal ] - else - ro(1:3) = ax(1:3) - ! we need to deal with the 180 degree case - ro(4) = merge(IEEE_value(ro(4),IEEE_positive_inf),tan(ax(4)*0.5 ),abs(ax(4)-PI) < thr) - end if - -end function ax2ro - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert axis angle pair to quaternion -!--------------------------------------------------------------------------------------------------- -pure function ax2qu(ax) result(qu) - use prec, only: & - dEq0 - - implicit none - real(pReal), intent(in), dimension(4) :: ax - type(quaternion) :: qu - - real(pReal) :: c, s - - - if (dEq0(ax(4))) then - qu = quaternion([ 1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal ]) - else - c = cos(ax(4)*0.5) - s = sin(ax(4)*0.5) - qu = quaternion([ c, ax(1)*s, ax(2)*s, ax(3)*s ]) - end if - -end function ax2qu - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert Rodrigues vector to homochoric -!--------------------------------------------------------------------------------------------------- -pure function ro2ho(ro) result(ho) - use, intrinsic :: IEEE_ARITHMETIC, only: & - IEEE_is_finite - use prec, only: & - dEq0 - use math, only: & - PI - - implicit none - real(pReal), intent(in), dimension(4) :: ro - real(pReal), dimension(3) :: ho - - real(pReal) :: f - - if (dEq0(norm2(ro(1:3)))) then - ho = [ 0.0, 0.0, 0.0 ] - else - f = merge(2.0*atan(ro(4)) - sin(2.0*atan(ro(4))),PI, IEEE_is_finite(ro(4))) - ho = ro(1:3) * (0.75_pReal*f)**(1.0/3.0) - end if - -end function ro2ho - - !--------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University !> @brief convert unit quaternion to rotation matrix !--------------------------------------------------------------------------------------------------- pure function qu2om(qu) result(om) - implicit none type(quaternion), intent(in) :: qu real(pReal), dimension(3,3) :: om @@ -724,39 +257,36 @@ end function qu2om !--------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University -!> @brief convert rotation matrix to a unit quaternion +!> @brief convert unit quaternion to Euler angles !--------------------------------------------------------------------------------------------------- -function om2qu(om) result(qu) +pure function qu2eu(qu) result(eu) use prec, only: & - dEq + dEq0 + use math, only: & + PI - implicit none - real(pReal), intent(in), dimension(3,3) :: om - type(quaternion) :: qu + type(quaternion), intent(in) :: qu + real(pReal), dimension(3) :: eu - real(pReal), dimension(4) :: qu_A - real(pReal), dimension(4) :: s - - s = [+om(1,1) +om(2,2) +om(3,3) +1.0_pReal, & - +om(1,1) -om(2,2) -om(3,3) +1.0_pReal, & - -om(1,1) +om(2,2) -om(3,3) +1.0_pReal, & - -om(1,1) -om(2,2) +om(3,3) +1.0_pReal] - - qu_A = sqrt(max(s,0.0_pReal))*0.5_pReal*[1.0_pReal,P,P,P] - qu_A = qu_A/norm2(qu_A) - - if(any(dEq(abs(qu_A),1.0_pReal,1.0e-15_pReal))) & - where (.not.(dEq(abs(qu_A),1.0_pReal,1.0e-15_pReal))) qu_A = 0.0_pReal - - if (om(3,2) < om(2,3)) qu_A(2) = -qu_A(2) - if (om(1,3) < om(3,1)) qu_A(3) = -qu_A(3) - if (om(2,1) < om(1,2)) qu_A(4) = -qu_A(4) + real(pReal) :: q12, q03, chi, chiInv - qu = quaternion(qu_A) - !qu_A = om2ax(om) - !if(any(qu_A(1:3) * [qu%x,qu%y,qu%z] < 0.0)) print*, 'sign error' + q03 = qu%w**2+qu%z**2 + q12 = qu%x**2+qu%y**2 + chi = sqrt(q03*q12) + + degenerated: if (dEq0(chi)) then + eu = merge([atan2(-P*2.0*qu%w*qu%z,qu%w**2-qu%z**2), 0.0_pReal, 0.0_pReal], & + [atan2(2.0*qu%x*qu%y,qu%x**2-qu%y**2), PI, 0.0_pReal], & + dEq0(q12)) + else degenerated + chiInv = 1.0/chi + eu = [atan2((-P*qu%w*qu%y+qu%x*qu%z)*chi, (-P*qu%w*qu%x-qu%y*qu%z)*chi ), & + atan2( 2.0*chi, q03-q12 ), & + atan2(( P*qu%w*qu%y+qu%x*qu%z)*chi, (-P*qu%w*qu%x+qu%y*qu%z)*chi )] + endif degenerated + where(eu<0.0_pReal) eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI]) -end function om2qu +end function qu2eu !--------------------------------------------------------------------------------------------------- @@ -771,15 +301,13 @@ pure function qu2ax(qu) result(ax) PI, & math_clip - implicit none type(quaternion), intent(in) :: qu real(pReal), dimension(4) :: ax real(pReal) :: omega, s - ! if the angle equals zero, then we return the rotation axis as [001] - if (dEq0(sqrt(qu%x**2+qu%y**2+qu%z**2))) then - ax = [ 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal ] + if (dEq0(qu%x**2+qu%y**2+qu%z**2)) then + ax = [ 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal ] ! axis = [001] elseif (dNeq0(qu%w)) then s = sign(1.0_pReal,qu%w)/sqrt(qu%x**2+qu%y**2+qu%z**2) omega = 2.0_pReal * acos(math_clip(qu%w,-1.0_pReal,1.0_pReal)) @@ -835,7 +363,6 @@ pure function qu2ho(qu) result(ho) use math, only: & math_clip - implicit none type(quaternion), intent(in) :: qu real(pReal), dimension(3) :: ho @@ -856,66 +383,109 @@ end function qu2ho !--------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University -!> @brief convert homochoric to cubochoric +!> @brief convert unit quaternion to cubochoric !--------------------------------------------------------------------------------------------------- -function ho2cu(ho) result(cu) - use Lambert, only: & - LambertBallToCube +function qu2cu(qu) result(cu) + + type(quaternion), intent(in) :: qu + real(pReal), dimension(3) :: cu - implicit none - real(pReal), intent(in), dimension(3) :: ho - real(pReal), dimension(3) :: cu + cu = ho2cu(qu2ho(qu)) - cu = LambertBallToCube(ho) +end function qu2cu -end function ho2cu + +!--------------------------------------------------------------------------------------------------- +!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH +!> @brief convert rotation matrix to cubochoric +!> @details the original formulation (direct conversion) had (numerical?) issues +!--------------------------------------------------------------------------------------------------- +function om2qu(om) result(qu) + + real(pReal), intent(in), dimension(3,3) :: om + type(quaternion) :: qu + + qu = eu2qu(om2eu(om)) + +end function om2qu !--------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University -!> @brief convert cubochoric to homochoric +!> @brief orientation matrix to Euler angles !--------------------------------------------------------------------------------------------------- -function cu2ho(cu) result(ho) - use Lambert, only: & - LambertCubeToBall +pure function om2eu(om) result(eu) + use math, only: & + PI - implicit none - real(pReal), intent(in), dimension(3) :: cu - real(pReal), dimension(3) :: ho + real(pReal), intent(in), dimension(3,3) :: om + real(pReal), dimension(3) :: eu + real(pReal) :: zeta + + if (abs(om(3,3)) < 1.0_pReal) then + zeta = 1.0_pReal/sqrt(1.0_pReal-om(3,3)**2.0_pReal) + eu = [atan2(om(3,1)*zeta,-om(3,2)*zeta), & + acos(om(3,3)), & + atan2(om(1,3)*zeta, om(2,3)*zeta)] + else + eu = [ atan2( om(1,2),om(1,1)), 0.5*PI*(1-om(3,3)),0.0_pReal ] + end if - ho = LambertCubeToBall(cu) - -end function cu2ho + where(eu<0.0_pReal) eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI]) + +end function om2eu !--------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University -!> @brief convert Rodrigues vector to Euler angles +!> @brief convert orientation matrix to axis angle pair !--------------------------------------------------------------------------------------------------- -pure function ro2eu(ro) result(eu) +function om2ax(om) result(ax) + use prec, only: & + dEq0, & + cEq, & + dNeq0 + use IO, only: & + IO_error + use math, only: & + math_clip, & + math_trace33 - implicit none - real(pReal), intent(in), dimension(4) :: ro - real(pReal), dimension(3) :: eu - - eu = om2eu(ro2om(ro)) + real(pReal), intent(in) :: om(3,3) + real(pReal) :: ax(4) + + real(pReal) :: t + real(pReal), dimension(3) :: Wr, Wi + real(pReal), dimension(10) :: WORK + real(pReal), dimension(3,3) :: VR, devNull, o + integer :: INFO, LWORK, i + + external :: dgeev,sgeev + + o = om + + ! first get the rotation angle + t = 0.5_pReal * (math_trace33(om) - 1.0) + ax(4) = acos(math_clip(t,-1.0_pReal,1.0_pReal)) + + if (dEq0(ax(4))) then + ax(1:3) = [ 0.0, 0.0, 1.0 ] + else + ! set some initial LAPACK variables + INFO = 0 + ! first initialize the parameters for the LAPACK DGEEV routines + LWORK = 20 + + ! call the eigenvalue solver + call dgeev('N','V',3,o,3,Wr,Wi,devNull,3,VR,3,WORK,LWORK,INFO) + if (INFO /= 0) call IO_error(0,ext_msg='Error in om2ax DGEEV return not zero') + i = maxloc(merge(1.0_pReal,0.0_pReal,cEq(cmplx(Wr,Wi,pReal),cmplx(1.0_pReal,0.0_pReal,pReal),tol=1.0e-14_pReal)),dim=1) ! poor substitute for findloc + ax(1:3) = VR(1:3,i) + where ( dNeq0([om(2,3)-om(3,2), om(3,1)-om(1,3), om(1,2)-om(2,1)])) & + ax(1:3) = sign(ax(1:3),-P *[om(2,3)-om(3,2), om(3,1)-om(1,3), om(1,2)-om(2,1)]) + endif -end function ro2eu - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert Euler angles to homochoric -!--------------------------------------------------------------------------------------------------- -pure function eu2ho(eu) result(ho) - - implicit none - real(pReal), intent(in), dimension(3) :: eu - real(pReal), dimension(3) :: ho - - ho = ax2ho(eu2ax(eu)) - -end function eu2ho +end function om2ax !--------------------------------------------------------------------------------------------------- @@ -924,7 +494,6 @@ end function eu2ho !--------------------------------------------------------------------------------------------------- pure function om2ro(om) result(ro) - implicit none real(pReal), intent(in), dimension(3,3) :: om real(pReal), dimension(4) :: ro @@ -939,133 +508,12 @@ end function om2ro !--------------------------------------------------------------------------------------------------- function om2ho(om) result(ho) - implicit none real(pReal), intent(in), dimension(3,3) :: om real(pReal), dimension(3) :: ho ho = ax2ho(om2ax(om)) -end function om2ho - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert axis angle pair to Euler angles -!--------------------------------------------------------------------------------------------------- -pure function ax2eu(ax) result(eu) - - implicit none - real(pReal), intent(in), dimension(4) :: ax - real(pReal), dimension(3) :: eu - - eu = om2eu(ax2om(ax)) - -end function ax2eu - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert Rodrigues vector to rotation matrix -!--------------------------------------------------------------------------------------------------- -pure function ro2om(ro) result(om) - - implicit none - real(pReal), intent(in), dimension(4) :: ro - real(pReal), dimension(3,3) :: om - - om = ax2om(ro2ax(ro)) - -end function ro2om - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert Rodrigues vector to unit quaternion -!--------------------------------------------------------------------------------------------------- -pure function ro2qu(ro) result(qu) - - implicit none - real(pReal), intent(in), dimension(4) :: ro - type(quaternion) :: qu - - qu = ax2qu(ro2ax(ro)) - -end function ro2qu - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert homochoric to Euler angles -!--------------------------------------------------------------------------------------------------- -pure function ho2eu(ho) result(eu) - - implicit none - real(pReal), intent(in), dimension(3) :: ho - real(pReal), dimension(3) :: eu - - eu = ax2eu(ho2ax(ho)) - -end function ho2eu - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert homochoric to rotation matrix -!--------------------------------------------------------------------------------------------------- -pure function ho2om(ho) result(om) - - implicit none - real(pReal), intent(in), dimension(3) :: ho - real(pReal), dimension(3,3) :: om - - om = ax2om(ho2ax(ho)) - -end function ho2om - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert homochoric to Rodrigues vector -!--------------------------------------------------------------------------------------------------- -pure function ho2ro(ho) result(ro) - - implicit none - real(pReal), intent(in), dimension(3) :: ho - real(pReal), dimension(4) :: ro - - ro = ax2ro(ho2ax(ho)) - -end function ho2ro - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert homochoric to unit quaternion -!--------------------------------------------------------------------------------------------------- -pure function ho2qu(ho) result(qu) - - implicit none - real(pReal), intent(in), dimension(3) :: ho - type(quaternion) :: qu - - qu = ax2qu(ho2ax(ho)) - -end function ho2qu - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert Euler angles to cubochoric -!--------------------------------------------------------------------------------------------------- -function eu2cu(eu) result(cu) - - implicit none - real(pReal), intent(in), dimension(3) :: eu - real(pReal), dimension(3) :: cu - - cu = ho2cu(eu2ho(eu)) - -end function eu2cu +end function om2ho !--------------------------------------------------------------------------------------------------- @@ -1074,7 +522,6 @@ end function eu2cu !--------------------------------------------------------------------------------------------------- function om2cu(om) result(cu) - implicit none real(pReal), intent(in), dimension(3,3) :: om real(pReal), dimension(3) :: cu @@ -1083,13 +530,279 @@ function om2cu(om) result(cu) end function om2cu +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief Euler angles to unit quaternion +!--------------------------------------------------------------------------------------------------- +pure function eu2qu(eu) result(qu) + + real(pReal), intent(in), dimension(3) :: eu + type(quaternion) :: qu + real(pReal), dimension(3) :: ee + real(pReal) :: cPhi, sPhi + + ee = 0.5_pReal*eu + + cPhi = cos(ee(2)) + sPhi = sin(ee(2)) + + qu = quaternion([ cPhi*cos(ee(1)+ee(3)), & + -P*sPhi*cos(ee(1)-ee(3)), & + -P*sPhi*sin(ee(1)-ee(3)), & + -P*cPhi*sin(ee(1)+ee(3))]) + if(qu%w < 0.0_pReal) qu = qu%homomorphed() + +end function eu2qu + + +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief Euler angles to orientation matrix +!--------------------------------------------------------------------------------------------------- +pure function eu2om(eu) result(om) + use prec, only: & + dEq0 + + real(pReal), intent(in), dimension(3) :: eu + real(pReal), dimension(3,3) :: om + + real(pReal), dimension(3) :: c, s + + c = cos(eu) + s = sin(eu) + + om(1,1) = c(1)*c(3)-s(1)*s(3)*c(2) + om(1,2) = s(1)*c(3)+c(1)*s(3)*c(2) + om(1,3) = s(3)*s(2) + om(2,1) = -c(1)*s(3)-s(1)*c(3)*c(2) + om(2,2) = -s(1)*s(3)+c(1)*c(3)*c(2) + om(2,3) = c(3)*s(2) + om(3,1) = s(1)*s(2) + om(3,2) = -c(1)*s(2) + om(3,3) = c(2) + + where(dEq0(om)) om = 0.0_pReal + +end function eu2om + + +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief convert euler to axis angle +!--------------------------------------------------------------------------------------------------- +pure function eu2ax(eu) result(ax) + use prec, only: & + dEq0, & + dEq + use math, only: & + PI + + real(pReal), intent(in), dimension(3) :: eu + real(pReal), dimension(4) :: ax + + real(pReal) :: t, delta, tau, alpha, sigma + + t = tan(eu(2)*0.5) + sigma = 0.5*(eu(1)+eu(3)) + delta = 0.5*(eu(1)-eu(3)) + tau = sqrt(t**2+sin(sigma)**2) + + alpha = merge(PI, 2.0*atan(tau/cos(sigma)), dEq(sigma,PI*0.5_pReal,tol=1.0e-15_pReal)) + + if (dEq0(alpha)) then ! return a default identity axis-angle pair + ax = [ 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal ] + else + ax(1:3) = -P/tau * [ t*cos(delta), t*sin(delta), sin(sigma) ] ! passive axis-angle pair so a minus sign in front + ax(4) = alpha + if (alpha < 0.0) ax = -ax ! ensure alpha is positive + end if + +end function eu2ax + + +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief Euler angles to Rodrigues vector +!--------------------------------------------------------------------------------------------------- +pure function eu2ro(eu) result(ro) + use prec, only: & + dEq0 + use, intrinsic :: IEEE_ARITHMETIC, only: & + IEEE_value, & + IEEE_positive_inf + use math, only: & + PI + + real(pReal), intent(in), dimension(3) :: eu + real(pReal), dimension(4) :: ro + + ro = eu2ax(eu) + if (ro(4) >= PI) then + ro(4) = IEEE_value(ro(4),IEEE_positive_inf) + elseif(dEq0(ro(4))) then + ro = [ 0.0_pReal, 0.0_pReal, P, 0.0_pReal ] + else + ro(4) = tan(ro(4)*0.5) + end if + +end function eu2ro + + +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief convert Euler angles to homochoric +!--------------------------------------------------------------------------------------------------- +pure function eu2ho(eu) result(ho) + + real(pReal), intent(in), dimension(3) :: eu + real(pReal), dimension(3) :: ho + + ho = ax2ho(eu2ax(eu)) + +end function eu2ho + + +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief convert Euler angles to cubochoric +!--------------------------------------------------------------------------------------------------- +function eu2cu(eu) result(cu) + + real(pReal), intent(in), dimension(3) :: eu + real(pReal), dimension(3) :: cu + + cu = ho2cu(eu2ho(eu)) + +end function eu2cu + + +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief convert axis angle pair to quaternion +!--------------------------------------------------------------------------------------------------- +pure function ax2qu(ax) result(qu) + use prec, only: & + dEq0 + + real(pReal), intent(in), dimension(4) :: ax + type(quaternion) :: qu + + real(pReal) :: c, s + + + if (dEq0(ax(4))) then + qu = quaternion([ 1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal ]) + else + c = cos(ax(4)*0.5) + s = sin(ax(4)*0.5) + qu = quaternion([ c, ax(1)*s, ax(2)*s, ax(3)*s ]) + end if + +end function ax2qu + + +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief convert axis angle pair to orientation matrix +!--------------------------------------------------------------------------------------------------- +pure function ax2om(ax) result(om) + + real(pReal), intent(in), dimension(4) :: ax + real(pReal), dimension(3,3) :: om + + real(pReal) :: q, c, s, omc + integer :: i + + c = cos(ax(4)) + s = sin(ax(4)) + omc = 1.0-c + + forall(i=1:3) om(i,i) = ax(i)**2*omc + c + + q = omc*ax(1)*ax(2) + om(1,2) = q + s*ax(3) + om(2,1) = q - s*ax(3) + + q = omc*ax(2)*ax(3) + om(2,3) = q + s*ax(1) + om(3,2) = q - s*ax(1) + + q = omc*ax(3)*ax(1) + om(3,1) = q + s*ax(2) + om(1,3) = q - s*ax(2) + + if (P > 0.0) om = transpose(om) + +end function ax2om + + +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief convert axis angle pair to Euler angles +!--------------------------------------------------------------------------------------------------- +pure function ax2eu(ax) result(eu) + + real(pReal), intent(in), dimension(4) :: ax + real(pReal), dimension(3) :: eu + + eu = om2eu(ax2om(ax)) + +end function ax2eu + + +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief convert axis angle pair to Rodrigues vector +!--------------------------------------------------------------------------------------------------- +pure function ax2ro(ax) result(ro) + use, intrinsic :: IEEE_ARITHMETIC, only: & + IEEE_value, & + IEEE_positive_inf + use prec, only: & + dEq0 + use math, only: & + PI + + real(pReal), intent(in), dimension(4) :: ax + real(pReal), dimension(4) :: ro + + real(pReal), parameter :: thr = 1.0E-7 + + if (dEq0(ax(4))) then + ro = [ 0.0_pReal, 0.0_pReal, P, 0.0_pReal ] + else + ro(1:3) = ax(1:3) + ! we need to deal with the 180 degree case + ro(4) = merge(IEEE_value(ro(4),IEEE_positive_inf),tan(ax(4)*0.5 ),abs(ax(4)-PI) < thr) + end if + +end function ax2ro + + +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief convert axis angle pair to homochoric +!--------------------------------------------------------------------------------------------------- +pure function ax2ho(ax) result(ho) + + real(pReal), intent(in), dimension(4) :: ax + real(pReal), dimension(3) :: ho + + real(pReal) :: f + + f = 0.75 * ( ax(4) - sin(ax(4)) ) + f = f**(1.0/3.0) + ho = ax(1:3) * f + +end function ax2ho + + !--------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University !> @brief convert axis angle pair to cubochoric !--------------------------------------------------------------------------------------------------- function ax2cu(ax) result(cu) - implicit none real(pReal), intent(in), dimension(4) :: ax real(pReal), dimension(3) :: cu @@ -1098,13 +811,113 @@ function ax2cu(ax) result(cu) end function ax2cu +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief convert Rodrigues vector to unit quaternion +!--------------------------------------------------------------------------------------------------- +pure function ro2qu(ro) result(qu) + + real(pReal), intent(in), dimension(4) :: ro + type(quaternion) :: qu + + qu = ax2qu(ro2ax(ro)) + +end function ro2qu + + +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief convert Rodrigues vector to rotation matrix +!--------------------------------------------------------------------------------------------------- +pure function ro2om(ro) result(om) + + real(pReal), intent(in), dimension(4) :: ro + real(pReal), dimension(3,3) :: om + + om = ax2om(ro2ax(ro)) + +end function ro2om + + +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief convert Rodrigues vector to Euler angles +!--------------------------------------------------------------------------------------------------- +pure function ro2eu(ro) result(eu) + + real(pReal), intent(in), dimension(4) :: ro + real(pReal), dimension(3) :: eu + + eu = om2eu(ro2om(ro)) + +end function ro2eu + + +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief convert Rodrigues vector to axis angle pair +!--------------------------------------------------------------------------------------------------- +pure function ro2ax(ro) result(ax) + use, intrinsic :: IEEE_ARITHMETIC, only: & + IEEE_is_finite + use prec, only: & + dEq0 + use math, only: & + PI + + real(pReal), intent(in), dimension(4) :: ro + real(pReal), dimension(4) :: ax + + real(pReal) :: ta, angle + + ta = ro(4) + + if (dEq0(ta)) then + ax = [ 0.0, 0.0, 1.0, 0.0 ] + elseif (.not. IEEE_is_finite(ta)) then + ax = [ ro(1), ro(2), ro(3), PI ] + else + angle = 2.0*atan(ta) + ta = 1.0/norm2(ro(1:3)) + ax = [ ro(1)/ta, ro(2)/ta, ro(3)/ta, angle ] + end if + +end function ro2ax + + +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief convert Rodrigues vector to homochoric +!--------------------------------------------------------------------------------------------------- +pure function ro2ho(ro) result(ho) + use, intrinsic :: IEEE_ARITHMETIC, only: & + IEEE_is_finite + use prec, only: & + dEq0 + use math, only: & + PI + + real(pReal), intent(in), dimension(4) :: ro + real(pReal), dimension(3) :: ho + + real(pReal) :: f + + if (dEq0(norm2(ro(1:3)))) then + ho = [ 0.0, 0.0, 0.0 ] + else + f = merge(2.0*atan(ro(4)) - sin(2.0*atan(ro(4))),PI, IEEE_is_finite(ro(4))) + ho = ro(1:3) * (0.75_pReal*f)**(1.0/3.0) + end if + +end function ro2ho + + !--------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University !> @brief convert Rodrigues vector to cubochoric !--------------------------------------------------------------------------------------------------- function ro2cu(ro) result(cu) - implicit none real(pReal), intent(in), dimension(4) :: ro real(pReal), dimension(3) :: cu @@ -1115,32 +928,130 @@ end function ro2cu !--------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University -!> @brief convert unit quaternion to cubochoric +!> @brief convert homochoric to unit quaternion !--------------------------------------------------------------------------------------------------- -function qu2cu(qu) result(cu) - - implicit none - type(quaternion), intent(in) :: qu - real(pReal), dimension(3) :: cu +pure function ho2qu(ho) result(qu) - cu = ho2cu(qu2ho(qu)) + real(pReal), intent(in), dimension(3) :: ho + type(quaternion) :: qu -end function qu2cu + qu = ax2qu(ho2ax(ho)) + +end function ho2qu !--------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University -!> @brief convert cubochoric to Euler angles +!> @brief convert homochoric to rotation matrix !--------------------------------------------------------------------------------------------------- -function cu2eu(cu) result(eu) +pure function ho2om(ho) result(om) - implicit none - real(pReal), intent(in), dimension(3) :: cu + real(pReal), intent(in), dimension(3) :: ho + real(pReal), dimension(3,3) :: om + + om = ax2om(ho2ax(ho)) + +end function ho2om + + +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief convert homochoric to Euler angles +!--------------------------------------------------------------------------------------------------- +pure function ho2eu(ho) result(eu) + + real(pReal), intent(in), dimension(3) :: ho real(pReal), dimension(3) :: eu - eu = ho2eu(cu2ho(cu)) + eu = ax2eu(ho2ax(ho)) -end function cu2eu +end function ho2eu + + +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief convert homochoric to axis angle pair +!--------------------------------------------------------------------------------------------------- +pure function ho2ax(ho) result(ax) + use prec, only: & + dEq0 + + real(pReal), intent(in), dimension(3) :: ho + real(pReal), dimension(4) :: ax + + integer :: i + real(pReal) :: hmag_squared, s, hm + real(pReal), parameter, dimension(16) :: & + tfit = [ 1.0000000000018852_pReal, -0.5000000002194847_pReal, & + -0.024999992127593126_pReal, -0.003928701544781374_pReal, & + -0.0008152701535450438_pReal, -0.0002009500426119712_pReal, & + -0.00002397986776071756_pReal, -0.00008202868926605841_pReal, & + +0.00012448715042090092_pReal, -0.0001749114214822577_pReal, & + +0.0001703481934140054_pReal, -0.00012062065004116828_pReal, & + +0.000059719705868660826_pReal, -0.00001980756723965647_pReal, & + +0.000003953714684212874_pReal, -0.00000036555001439719544_pReal ] + + ! normalize h and store the magnitude + hmag_squared = sum(ho**2.0_pReal) + if (dEq0(hmag_squared)) then + ax = [ 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal ] + else + hm = hmag_squared + + ! convert the magnitude to the rotation angle + s = tfit(1) + tfit(2) * hmag_squared + do i=3,16 + hm = hm*hmag_squared + s = s + tfit(i) * hm + end do + ax = [ho/sqrt(hmag_squared), 2.0_pReal*acos(s)] + end if + +end function ho2ax + + +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief convert homochoric to Rodrigues vector +!--------------------------------------------------------------------------------------------------- +pure function ho2ro(ho) result(ro) + + real(pReal), intent(in), dimension(3) :: ho + real(pReal), dimension(4) :: ro + + ro = ax2ro(ho2ax(ho)) + +end function ho2ro + + +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief convert homochoric to cubochoric +!--------------------------------------------------------------------------------------------------- +function ho2cu(ho) result(cu) + use Lambert, only: & + LambertBallToCube + + real(pReal), intent(in), dimension(3) :: ho + real(pReal), dimension(3) :: cu + + cu = LambertBallToCube(ho) + +end function ho2cu + + +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief convert cubochoric to unit quaternion +!--------------------------------------------------------------------------------------------------- +function cu2qu(cu) result(qu) + + real(pReal), intent(in), dimension(3) :: cu + type(quaternion) :: qu + + qu = ho2qu(cu2ho(cu)) + +end function cu2qu !--------------------------------------------------------------------------------------------------- @@ -1149,7 +1060,6 @@ end function cu2eu !--------------------------------------------------------------------------------------------------- function cu2om(cu) result(om) - implicit none real(pReal), intent(in), dimension(3) :: cu real(pReal), dimension(3,3) :: om @@ -1158,13 +1068,26 @@ function cu2om(cu) result(om) end function cu2om +!--------------------------------------------------------------------------------------------------- +!> @author Marc De Graef, Carnegie Mellon University +!> @brief convert cubochoric to Euler angles +!--------------------------------------------------------------------------------------------------- +function cu2eu(cu) result(eu) + + real(pReal), intent(in), dimension(3) :: cu + real(pReal), dimension(3) :: eu + + eu = ho2eu(cu2ho(cu)) + +end function cu2eu + + !--------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University !> @brief convert cubochoric to axis angle pair !--------------------------------------------------------------------------------------------------- function cu2ax(cu) result(ax) - implicit none real(pReal), intent(in), dimension(3) :: cu real(pReal), dimension(4) :: ax @@ -1179,7 +1102,6 @@ end function cu2ax !--------------------------------------------------------------------------------------------------- function cu2ro(cu) result(ro) - implicit none real(pReal), intent(in), dimension(3) :: cu real(pReal), dimension(4) :: ro @@ -1190,16 +1112,18 @@ end function cu2ro !--------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University -!> @brief convert cubochoric to unit quaternion +!> @brief convert cubochoric to homochoric !--------------------------------------------------------------------------------------------------- -function cu2qu(cu) result(qu) +function cu2ho(cu) result(ho) + use Lambert, only: & + LambertCubeToBall - implicit none real(pReal), intent(in), dimension(3) :: cu - type(quaternion) :: qu + real(pReal), dimension(3) :: ho - qu = ho2qu(cu2ho(cu)) + ho = LambertCubeToBall(cu) + +end function cu2ho -end function cu2qu end module rotations