From b8285d57493e184387e1d966f42e0bc9f1b7274b Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Thu, 11 Apr 2019 17:18:32 -0400 Subject: [PATCH 01/29] restored orientation averaging capability --- python/damask/orientation.py | 102 ++++++++++++++++++++--------------- 1 file changed, 59 insertions(+), 43 deletions(-) diff --git a/python/damask/orientation.py b/python/damask/orientation.py index 7cb05af40..463dfeb1b 100644 --- a/python/damask/orientation.py +++ b/python/damask/orientation.py @@ -35,14 +35,6 @@ class Quaternion: """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) @@ -116,7 +108,7 @@ class Quaternion: return self else: return NotImplemented - + def __truediv__(self, other): """Divsion with quaternion or scalar""" @@ -180,6 +172,17 @@ class Quaternion: 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): + return [self.q]+list(self.p) + def normalize(self): d = self.magnitude() if d > 0.0: @@ -190,7 +193,7 @@ class Quaternion: def normalized(self): return self.copy().normalize() - + def conjugate(self): self.p = -self.p return self @@ -294,6 +297,10 @@ class Rotation: def asCubochoric(self): return qu2cu(self.quaternion.asArray()) + def asM(self): + """Intermediate representation useful fro 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 @@ -503,7 +510,7 @@ class Symmetry: otherOrder = Symmetry.lattices.index(other.lattice) return (myOrder > otherOrder) - (myOrder < otherOrder) - def symmetryOperations(self): + def symmetryOperations(self,who=[]): """List of symmetry operations as quaternions.""" if self.lattice == 'cubic': symQuats = [ @@ -570,7 +577,8 @@ class Symmetry: [ 1.0,0.0,0.0,0.0 ], ] - return [Rotation(q) for q in symQuats] + return list(map(Rotation, + np.array(symQuats)[np.atleast_1d(np.array(who)) if who != [] else range(len(symQuats))])) def inFZ(self,R): @@ -579,6 +587,8 @@ class Symmetry: Fundamental zone in Rodrigues space is point symmetric around origin. """ + if np.any(R == np.inf): return False + Rabs = abs(R[0:3]*R[3]) if self.lattice == 'cubic': @@ -1050,7 +1060,8 @@ class Orientation: def disorientation(self, other, - SST = True): + SST = True, + symmetries = False): """ Disorientation between myself and given other orientation. @@ -1076,15 +1087,18 @@ class Orientation: if breaker: break if breaker: break - return r + return (r, 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()) - def equivalentOrientations(self): + def equivalentOrientations(self, who=[]): """List of orientations which are symmetrically equivalent""" return [self.__class__(q*self.rotation,self.lattice) \ - for q in self.lattice.symmetry.symmetryOperations()] + for q in self.lattice.symmetry.symmetryOperations(who)] def relatedOrientations(self,model): """List of orientations related by the given orientation relationship""" @@ -1125,37 +1139,39 @@ class Orientation: return color - # @classmethod - # def average(cls, - # orientations, - # multiplicity = []): - # """ - # Average orientation + @classmethod + def average(cls, + orientations, + multiplicity = []): + """ + Average orientation - # ref: F. Landis Markley, Yang Cheng, John Lucas Crassidis, and Yaakov Oshman. - # Averaging Quaternions, - # Journal of Guidance, Control, and Dynamics, Vol. 30, No. 4 (2007), pp. 1193-1197. - # doi: 10.2514/1.28949 - # usage: - # a = Orientation(Eulers=np.radians([10, 10, 0]), symmetry='hexagonal') - # b = Orientation(Eulers=np.radians([20, 0, 0]), symmetry='hexagonal') - # avg = Orientation.average([a,b]) - # """ - # if not all(isinstance(item, Orientation) for item in orientations): - # raise TypeError("Only instances of Orientation can be averaged.") + 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): + print('got a list of {}'.format(orientations)) + raise TypeError("Only instances of Orientation can be averaged.") - # N = len(orientations) - # if multiplicity == [] or not multiplicity: - # multiplicity = np.ones(N,dtype='i') + N = len(orientations) + if multiplicity == [] or not multiplicity: + multiplicity = np.ones(N,dtype='i') - # 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) + ref = orientations[0] # take first as reference + for i,(o,n) in enumerate(zip(orientations,multiplicity)): + closest = o.equivalentOrientations(ref.disorientation(o,SST=False,symmetries=True)[2])[0] # select sym orientation with lowest misorientation + M = closest.rotation.asM() * n if i == 0 \ + else M + closest.rotation.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) + return Orientation(Rotation.fromQuaternion(np.real(vec.T[eig.argmax()]),acceptHomomorph = True), + ref.lattice) #################################################################################################### From 2190c3ef466f0ac08100276671393293d4b9b05d Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Thu, 11 Apr 2019 18:32:07 -0400 Subject: [PATCH 02/29] fixed serious disorientation bug; sorted transformation functions --- python/damask/orientation.py | 530 ++++++++++++++++++----------------- 1 file changed, 274 insertions(+), 256 deletions(-) diff --git a/python/damask/orientation.py b/python/damask/orientation.py index 463dfeb1b..f4df909c6 100644 --- a/python/damask/orientation.py +++ b/python/damask/orientation.py @@ -286,9 +286,9 @@ 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()) + return qu2ro(self.quaternion.asArray(),vector) def asHomochoric(self): """Homochoric vector: (h_1, h_2, h_3)""" @@ -1072,15 +1072,16 @@ class Orientation: #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() + r.inverse() breaker = self.lattice.symmetry.inFZ(r.asRodrigues()) \ and (not SST or other.lattice.symmetry.inDisorientationSST(r.asRodrigues())) if breaker: break @@ -1211,94 +1212,7 @@ def isone(a): def iszero(a): return np.isclose(a,0.0,atol=1.0e-12,rtol=0.0) - -def eu2om(eu): - """Euler angles to orientation matrix""" - c = np.cos(eu) - s = np.sin(eu) - - om = np.array([[+c[0]*c[2]-s[0]*s[2]*c[1], +s[0]*c[2]+c[0]*s[2]*c[1], +s[2]*s[1]], - [-c[0]*s[2]-s[0]*c[2]*c[1], -s[0]*s[2]+c[0]*c[2]*c[1], +c[2]*s[1]], - [+s[0]*s[1], -c[0]*s[1], +c[1] ]]) - - om[np.where(iszero(om))] = 0.0 - return om - - -def eu2ax(eu): - """Euler angles to axis angle""" - t = np.tan(eu[1]*0.5) - sigma = 0.5*(eu[0]+eu[2]) - delta = 0.5*(eu[0]-eu[2]) - tau = np.linalg.norm([t,np.sin(sigma)]) - alpha = np.pi if iszero(np.cos(sigma)) else \ - 2.0*np.arctan(tau/np.cos(sigma)) - - if iszero(alpha): - ax = np.array([ 0.0, 0.0, 1.0, 0.0 ]) - else: - ax = -P/tau * np.array([ t*np.cos(delta), t*np.sin(delta), np.sin(sigma) ]) # passive axis-angle pair so a minus sign in front - ax = np.append(ax,alpha) - if alpha < 0.0: ax *= -1.0 # ensure alpha is positive - - return ax - - -def eu2ro(eu): - """Euler angles to Rodrigues vector""" - ro = eu2ax(eu) # convert to axis angle representation - if ro[3] >= np.pi: # Differs from original implementation. check convention 5 - ro[3] = np.inf - elif iszero(ro[3]): - ro = np.array([ 0.0, 0.0, P, 0.0 ]) - else: - ro[3] = np.tan(ro[3]*0.5) - - 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 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 - 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)]) - - # 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 ax2om(ax): - """Axis angle to orientation matrix""" - c = np.cos(ax[3]) - s = np.sin(ax[3]) - omc = 1.0-c - om=np.diag(ax[0:3]**2*omc + c) - - for idx in [[0,1,2],[1,2,0],[2,0,1]]: - q = omc*ax[idx[0]] * ax[idx[1]] - om[idx[0],idx[1]] = q + s*ax[idx[2]] - om[idx[1],idx[0]] = q - s*ax[idx[2]] - - return om if P < 0.0 else om.T - +#---------- quaternion ---------- def qu2eu(qu): """Quaternion to Euler angles""" @@ -1318,113 +1232,6 @@ def qu2eu(qu): 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 ax2ho(ax): - """Axis angle to homochoric""" - f = (0.75 * ( ax[3] - np.sin(ax[3]) ))**(1.0/3.0) - ho = ax[0:3] * f - return ho - - -def ho2ax(ho): - """Homochoric to axis angle""" - tfit = np.array([+1.0000000000018852, -0.5000000002194847, - -0.024999992127593126, -0.003928701544781374, - -0.0008152701535450438, -0.0002009500426119712, - -0.00002397986776071756, -0.00008202868926605841, - +0.00012448715042090092, -0.0001749114214822577, - +0.0001703481934140054, -0.00012062065004116828, - +0.000059719705868660826, -0.00001980756723965647, - +0.000003953714684212874, -0.00000036555001439719544]) - # normalize h and store the magnitude - hmag_squared = np.sum(ho**2.) - if iszero(hmag_squared): - ax = np.array([ 0.0, 0.0, 1.0, 0.0 ]) - else: - hm = hmag_squared - - # convert the magnitude to the rotation angle - s = tfit[0] + tfit[1] * hmag_squared - for i in range(2,16): - hm *= hmag_squared - s += tfit[i] * hm - ax = np.append(ho/np.sqrt(hmag_squared),2.0*np.arccos(np.clip(s,-1.0,1.0))) - 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) @@ -1457,7 +1264,7 @@ def qu2ax(qu): return np.array(ax) -def qu2ro(qu): +def qu2ro(qu,vector=False): """Quaternion to Rodrigues vector""" if iszero(qu[0]): ro = [qu[1], qu[2], qu[3], np.inf] @@ -1466,7 +1273,7 @@ def qu2ro(qu): 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) + return np.array(ro[:3])*ro[3] if vector else np.array(ro) def qu2ho(qu): @@ -1483,19 +1290,120 @@ def qu2ho(qu): return ho -def ho2cu(ho): - """Homochoric to cubochoric""" - return Lambert.BallToCube(ho) +def qu2cu(qu): + """Quaternion to cubochoric""" + return ho2cu(qu2ho(qu)) -def cu2ho(cu): - """Cubochoric to homochoric""" - return Lambert.CubeToBall(cu) +#---------- orientation matrix ---------- -def ro2eu(ro): - """Rodrigues vector to orientation matrix""" - return om2eu(ro2om(ro)) +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 + 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)]) + + # 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 om2qu(om): + """ + Orientation matrix to quaternion + + The original formulation (direct conversion) had numerical issues + """ + return ax2qu(om2ax(om)) + + +def om2ro(om,vector=False): + """Orientation matrix to Rodriques vector""" + return eu2ro(om2eu(om,vector)) + + +def om2cu(om): + """Orientation matrix to cubochoric""" + return ho2cu(om2ho(om)) + + +def om2ho(om): + """Orientation matrix to homochoric""" + return ax2ho(om2ax(om)) + + +#---------- Euler angles ---------- + + +def eu2om(eu): + """Euler angles to orientation matrix""" + c = np.cos(eu) + s = np.sin(eu) + + om = np.array([[+c[0]*c[2]-s[0]*s[2]*c[1], +s[0]*c[2]+c[0]*s[2]*c[1], +s[2]*s[1]], + [-c[0]*s[2]-s[0]*c[2]*c[1], -s[0]*s[2]+c[0]*c[2]*c[1], +c[2]*s[1]], + [+s[0]*s[1], -c[0]*s[1], +c[1] ]]) + + om[np.where(iszero(om))] = 0.0 + return om + + +def eu2ax(eu): + """Euler angles to axis angle""" + t = np.tan(eu[1]*0.5) + sigma = 0.5*(eu[0]+eu[2]) + delta = 0.5*(eu[0]-eu[2]) + tau = np.linalg.norm([t,np.sin(sigma)]) + alpha = np.pi if iszero(np.cos(sigma)) else \ + 2.0*np.arctan(tau/np.cos(sigma)) + + if iszero(alpha): + ax = np.array([ 0.0, 0.0, 1.0, 0.0 ]) + else: + ax = -P/tau * np.array([ t*np.cos(delta), t*np.sin(delta), np.sin(sigma) ]) # passive axis-angle pair so a minus sign in front + ax = np.append(ax,alpha) + if alpha < 0.0: ax *= -1.0 # ensure alpha is positive + + return ax + + +def eu2ro(eu,vector=False): + """Euler angles to Rodrigues vector""" + ro = eu2ax(eu) # convert to axis angle representation + if ro[3] >= np.pi: # Differs from original implementation. check convention 5 + ro[3] = np.inf + elif iszero(ro[3]): + ro = np.array([ 0.0, 0.0, P, 0.0 ]) + else: + ro[3] = np.tan(ro[3]*0.5) + + return ro[:3] * ro[3] if vector else ro def eu2ho(eu): @@ -1503,21 +1411,108 @@ def eu2ho(eu): return ax2ho(eu2ax(eu)) -def om2ro(om): - """Orientation matrix to Rodriques vector""" - return eu2ro(om2eu(om)) +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 om2ho(om): - """Orientation matrix to homochoric""" - return ax2ho(om2ax(om)) +def eu2cu(eu): + """Euler angles to cubochoric""" + return ho2cu(eu2ho(eu)) +#---------- axis angle ---------- + + def ax2eu(ax): """Orientation matrix to Euler angles""" return om2eu(ax2om(ax)) +def ax2om(ax): + """Axis angle to orientation matrix""" + c = np.cos(ax[3]) + s = np.sin(ax[3]) + omc = 1.0-c + om=np.diag(ax[0:3]**2*omc + c) + + for idx in [[0,1,2],[1,2,0],[2,0,1]]: + q = omc*ax[idx[0]] * ax[idx[1]] + om[idx[0],idx[1]] = q + s*ax[idx[2]] + om[idx[1],idx[0]] = q - s*ax[idx[2]] + + return om if P < 0.0 else om.T + + +def ax2ro(ax,vector=False): + """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[:3])*ro[3] if vector else 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 ax2ho(ax): + """Axis angle to homochoric""" + f = (0.75 * ( ax[3] - np.sin(ax[3]) ))**(1.0/3.0) + ho = ax[0:3] * f + return ho + + +def ax2cu(ax): + """Axis angle to cubochoric""" + return ho2cu(ax2ho(ax)) + + +#---------- Rodrigues--Frank ---------- + + +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 ro2eu(ro): + """Rodrigues vector to orientation matrix""" + return om2eu(ro2om(ro)) + + def ro2om(ro): """Rodgrigues vector to orientation matrix""" return ax2om(ro2ax(ro)) @@ -1528,6 +1523,56 @@ def ro2qu(ro): return ax2qu(ro2ax(ro)) +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 ho2ax(ho): + """Homochoric to axis angle""" + tfit = np.array([+1.0000000000018852, -0.5000000002194847, + -0.024999992127593126, -0.003928701544781374, + -0.0008152701535450438, -0.0002009500426119712, + -0.00002397986776071756, -0.00008202868926605841, + +0.00012448715042090092, -0.0001749114214822577, + +0.0001703481934140054, -0.00012062065004116828, + +0.000059719705868660826, -0.00001980756723965647, + +0.000003953714684212874, -0.00000036555001439719544]) + # normalize h and store the magnitude + hmag_squared = np.sum(ho**2.) + if iszero(hmag_squared): + ax = np.array([ 0.0, 0.0, 1.0, 0.0 ]) + else: + hm = hmag_squared + + # convert the magnitude to the rotation angle + s = tfit[0] + tfit[1] * hmag_squared + for i in range(2,16): + hm *= hmag_squared + s += tfit[i] * hm + ax = np.append(ho/np.sqrt(hmag_squared),2.0*np.arccos(np.clip(s,-1.0,1.0))) + return ax + + +def ho2cu(ho): + """Homochoric to cubochoric""" + return Lambert.BallToCube(ho) + + def ho2eu(ho): """Homochoric to Euler angles""" return ax2eu(ho2ax(ho)) @@ -1538,9 +1583,9 @@ def ho2om(ho): return ax2om(ho2ax(ho)) -def ho2ro(ho): +def ho2ro(ho,vector=False): """Axis angle to Rodriques vector""" - return ax2ro(ho2ax(ho)) + return ax2ro(ho2ax(ho,vector)) def ho2qu(ho): @@ -1548,40 +1593,6 @@ def ho2qu(ho): 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)) @@ -1597,11 +1608,18 @@ def cu2ax(cu): return ho2ax(cu2ho(cu)) -def cu2ro(cu): +def cu2ro(cu,vector=False): """Cubochoric to Rodrigues vector""" - return ho2ro(cu2ho(cu)) + return ho2ro(cu2ho(cu,vector)) def cu2qu(cu): """Cubochoric to quaternion""" return ho2qu(cu2ho(cu)) + + +def cu2ho(cu): + """Cubochoric to homochoric""" + return Lambert.CubeToBall(cu) + + From 3d6eb76da3992d64a37edb170e77beb7efb6c696 Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Thu, 11 Apr 2019 19:07:41 -0400 Subject: [PATCH 03/29] equivalentOrientations accepts scalar argument to directly return single rotation object (not single element list) --- python/damask/orientation.py | 23 +++++++++++++++++------ 1 file changed, 17 insertions(+), 6 deletions(-) diff --git a/python/damask/orientation.py b/python/damask/orientation.py index f4df909c6..227f02cda 100644 --- a/python/damask/orientation.py +++ b/python/damask/orientation.py @@ -510,7 +510,7 @@ class Symmetry: otherOrder = Symmetry.lattices.index(other.lattice) return (myOrder > otherOrder) - (myOrder < otherOrder) - def symmetryOperations(self,who=[]): + def symmetryOperations(self,members=[]): """List of symmetry operations as quaternions.""" if self.lattice == 'cubic': symQuats = [ @@ -577,8 +577,14 @@ class Symmetry: [ 1.0,0.0,0.0,0.0 ], ] - return list(map(Rotation, - np.array(symQuats)[np.atleast_1d(np.array(who)) if who != [] else range(len(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): @@ -1096,10 +1102,15 @@ class Orientation: def inFZ(self): return self.lattice.symmetry.inFZ(self.rotation.asRodrigues()) - def equivalentOrientations(self, who=[]): + 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(who)] + 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""" From 3a7408a2131a825d2c4b4b2103031896c4a46e1c Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Fri, 12 Apr 2019 09:02:00 -0400 Subject: [PATCH 04/29] averaging possible for rotations and orientations; Rodrigues 3Dvector output only at top-level; code reordering --- python/damask/orientation.py | 186 ++++++++++++++++++++--------------- 1 file changed, 104 insertions(+), 82 deletions(-) diff --git a/python/damask/orientation.py b/python/damask/orientation.py index 227f02cda..fde168d27 100644 --- a/python/damask/orientation.py +++ b/python/damask/orientation.py @@ -260,6 +260,63 @@ 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).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 + + 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.__class__(self.quaternion.conjugated()) + + + def misorientation(self,other): + """Misorientation""" + return self.__class__(other.quaternion*self.quaternion.conjugated()) + ################################################################################################ # convert to different orientation representations (numpy arrays) @@ -288,7 +345,8 @@ class Rotation: def asRodrigues(self,vector=False): """Rodrigues-Frank vector: ([n_1, n_2, n_3], tan(ω/2))""" - return qu2ro(self.quaternion.asArray(),vector) + 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)""" @@ -298,7 +356,7 @@ class Rotation: return qu2cu(self.quaternion.asArray()) def asM(self): - """Intermediate representation useful fro quaternion averaging (see F. Landis Markley et al.)""" + """Intermediate representation supporting quaternion averaging (see F. Landis Markley et al.)""" return self.quaternion.asM() @@ -409,63 +467,32 @@ class Rotation: return cls(ho2qu(ho)) - def __mul__(self, other): + @classmethod + def average(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 """ - 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 Rotation.fromQuaternion(np.real(vec.T[eig.argmax()]),acceptHomomorph = True) - def misorientation(self,other): - """Misorientation""" - return self.__class__(other.quaternion*self.quaternion.conjugated()) - # ****************************************************************************************** class Symmetry: @@ -511,7 +538,7 @@ class Symmetry: return (myOrder > otherOrder) - (myOrder < otherOrder) def symmetryOperations(self,members=[]): - """List of symmetry operations as quaternions.""" + """List (or single element) of symmetry operations as rotations.""" if self.lattice == 'cubic': symQuats = [ [ 1.0, 0.0, 0.0, 0.0 ], @@ -880,7 +907,7 @@ class Lattice: # Greninger--Troiano' orientation relationship for fcc <-> bcc transformation # from Y. He et al./Journal of Applied Crystallography (2006). 39, 72-81 - GTdash = {'mapping':{'fcc':0,'bcc':1}, + GTprime = {'mapping':{'fcc':0,'bcc':1}, 'planes': np.array([ [[ 7, 17, 17],[ 12, 5, 17]], [[ 17, 7, 17],[ 17, 12, 5]], @@ -1006,7 +1033,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] @@ -1154,7 +1181,7 @@ class Orientation: @classmethod def average(cls, orientations, - multiplicity = []): + weights = []): """ Average orientation @@ -1168,22 +1195,17 @@ class Orientation: avg = Orientation.average([a,b]) """ if not all(isinstance(item, Orientation) for item in orientations): - print('got a list of {}'.format(orientations)) raise TypeError("Only instances of Orientation can be averaged.") - N = len(orientations) - if multiplicity == [] or not multiplicity: - multiplicity = np.ones(N,dtype='i') + 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 - ref = orientations[0] # take first as reference - for i,(o,n) in enumerate(zip(orientations,multiplicity)): - closest = o.equivalentOrientations(ref.disorientation(o,SST=False,symmetries=True)[2])[0] # select sym orientation with lowest misorientation - M = closest.rotation.asM() * n if i == 0 \ - else M + closest.rotation.asM() * n # noqa add (multiples) of this orientation to average noqa - eig, vec = np.linalg.eig(M/N) - - return Orientation(Rotation.fromQuaternion(np.real(vec.T[eig.argmax()]),acceptHomomorph = True), - ref.lattice) + return Orientation(Rotation.average(closest,weights),ref.lattice) #################################################################################################### @@ -1275,7 +1297,7 @@ def qu2ax(qu): return np.array(ax) -def qu2ro(qu,vector=False): +def qu2ro(qu): """Quaternion to Rodrigues vector""" if iszero(qu[0]): ro = [qu[1], qu[2], qu[3], np.inf] @@ -1284,7 +1306,7 @@ def qu2ro(qu,vector=False): 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[:3])*ro[3] if vector else np.array(ro) + return np.array(ro) def qu2ho(qu): @@ -1354,9 +1376,9 @@ def om2qu(om): return ax2qu(om2ax(om)) -def om2ro(om,vector=False): +def om2ro(om): """Orientation matrix to Rodriques vector""" - return eu2ro(om2eu(om,vector)) + return eu2ro(om2eu(om)) def om2cu(om): @@ -1404,7 +1426,7 @@ def eu2ax(eu): return ax -def eu2ro(eu,vector=False): +def eu2ro(eu): """Euler angles to Rodrigues vector""" ro = eu2ax(eu) # convert to axis angle representation if ro[3] >= np.pi: # Differs from original implementation. check convention 5 @@ -1414,7 +1436,7 @@ def eu2ro(eu,vector=False): else: ro[3] = np.tan(ro[3]*0.5) - return ro[:3] * ro[3] if vector else ro + return ro def eu2ho(eu): @@ -1463,7 +1485,7 @@ def ax2om(ax): return om if P < 0.0 else om.T -def ax2ro(ax,vector=False): +def ax2ro(ax): """Axis angle to Rodrigues vector""" if iszero(ax[3]): ro = [ 0.0, 0.0, P, 0.0 ] @@ -1473,7 +1495,7 @@ def ax2ro(ax,vector=False): 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[:3])*ro[3] if vector else np.array(ro) + return np.array(ro) def ax2qu(ax): @@ -1594,9 +1616,9 @@ def ho2om(ho): return ax2om(ho2ax(ho)) -def ho2ro(ho,vector=False): +def ho2ro(ho): """Axis angle to Rodriques vector""" - return ax2ro(ho2ax(ho,vector)) + return ax2ro(ho2ax(ho)) def ho2qu(ho): @@ -1619,9 +1641,9 @@ def cu2ax(cu): return ho2ax(cu2ho(cu)) -def cu2ro(cu,vector=False): +def cu2ro(cu): """Cubochoric to Rodrigues vector""" - return ho2ro(cu2ho(cu,vector)) + return ho2ro(cu2ho(cu)) def cu2qu(cu): From d842902a76b41c19a0ccffacc721617c9c470739 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 16 Apr 2019 21:47:29 +0200 Subject: [PATCH 05/29] consistent ordering - qu, om, eu, ax, ro, ho, cu --- python/damask/orientation.py | 198 +++++++++++++++++------------------ 1 file changed, 97 insertions(+), 101 deletions(-) diff --git a/python/damask/orientation.py b/python/damask/orientation.py index fde168d27..007aecfd9 100644 --- a/python/damask/orientation.py +++ b/python/damask/orientation.py @@ -1247,6 +1247,20 @@ def iszero(a): #---------- 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 @@ -1265,19 +1279,6 @@ def qu2eu(qu): 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 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): """ @@ -1330,6 +1331,14 @@ def qu2cu(qu): #---------- orientation matrix ---------- +def om2qu(om): + """ + Orientation matrix to quaternion + + The original formulation (direct conversion) had numerical issues + """ + return ax2qu(om2ax(om)) + def om2eu(om): """Euler angles to orientation matrix""" @@ -1365,25 +1374,11 @@ def om2ax(om): ax[0:3] = np.where(iszero(diagDelta), ax[0:3],np.abs(ax[0:3])*np.sign(-P*diagDelta)) return np.array(ax) - - -def om2qu(om): - """ - Orientation matrix to quaternion - - The original formulation (direct conversion) had numerical issues - """ - return ax2qu(om2ax(om)) def om2ro(om): """Orientation matrix to Rodriques vector""" return eu2ro(om2eu(om)) - - -def om2cu(om): - """Orientation matrix to cubochoric""" - return ho2cu(om2ho(om)) def om2ho(om): @@ -1391,9 +1386,26 @@ def om2ho(om): 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.homomorph() !ToDo: Check with original + return qu + + def eu2om(eu): """Euler angles to orientation matrix""" c = np.cos(eu) @@ -1444,19 +1456,6 @@ def eu2ho(eu): return ax2ho(eu2ax(eu)) -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 eu2cu(eu): """Euler angles to cubochoric""" return ho2cu(eu2ho(eu)) @@ -1464,10 +1463,16 @@ def eu2cu(eu): #---------- axis angle ---------- - -def ax2eu(ax): - """Orientation matrix to Euler angles""" - return om2eu(ax2om(ax)) +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 ax2om(ax): @@ -1484,6 +1489,11 @@ 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 ax2ro(ax): """Axis angle to Rodrigues vector""" @@ -1498,18 +1508,6 @@ def ax2ro(ax): 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 ax2ho(ax): """Axis angle to homochoric""" f = (0.75 * ( ax[3] - np.sin(ax[3]) ))**(1.0/3.0) @@ -1524,6 +1522,20 @@ def ax2cu(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""" @@ -1541,21 +1553,6 @@ def ro2ax(ro): return np.array(ax) -def ro2eu(ro): - """Rodrigues vector to orientation matrix""" - return om2eu(ro2om(ro)) - - -def ro2om(ro): - """Rodgrigues vector to orientation matrix""" - return ax2om(ro2ax(ro)) - - -def ro2qu(ro): - """Rodrigues vector to quaternion""" - return ax2qu(ro2ax(ro)) - - def ro2ho(ro): """Rodrigues vector to homochoric""" if iszero(np.sum(ro[0:3]**2.0)): @@ -1574,7 +1571,21 @@ def ro2cu(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, @@ -1601,34 +1612,21 @@ def ho2ax(ho): return ax -def ho2cu(ho): - """Homochoric to cubochoric""" - return Lambert.BallToCube(ho) - - -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 ho2cu(ho): + """Homochoric to cubochoric""" + return Lambert.BallToCube(ho) -def cu2eu(cu): - """Cubochoric to Euler angles""" - return ho2eu(cu2ho(cu)) +#---------- cubochoric ---------- + +def cu2qu(cu): + """Cubochoric to quaternion""" + return ho2qu(cu2ho(cu)) def cu2om(cu): @@ -1636,6 +1634,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)) @@ -1646,13 +1649,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) - - From 5d23a61fb0a381a6a9a3192c1b9be634b43213b2 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 16 Apr 2019 22:17:56 +0200 Subject: [PATCH 06/29] same order as in the python module ... and only one "implicit none" --- src/rotations.f90 | 1307 ++++++++++++++++++++++----------------------- 1 file changed, 626 insertions(+), 681 deletions(-) diff --git a/src/rotations.f90 b/src/rotations.f90 index b899adacb..aed4c5773 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,7 +301,6 @@ pure function qu2ax(qu) result(ax) PI, & math_clip - implicit none type(quaternion), intent(in) :: qu real(pReal), dimension(4) :: ax @@ -835,7 +364,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 +384,129 @@ 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 ho2cu +end function qu2cu !--------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University -!> @brief convert cubochoric to homochoric +!> @brief convert rotation matrix to a unit quaternion !--------------------------------------------------------------------------------------------------- -function cu2ho(cu) result(ho) - use Lambert, only: & - LambertCubeToBall +function om2qu(om) result(qu) + use prec, only: & + dEq - implicit none - real(pReal), intent(in), dimension(3) :: cu - real(pReal), dimension(3) :: ho + real(pReal), intent(in), dimension(3,3) :: om + type(quaternion) :: qu + + real(pReal), dimension(4) :: qu_A + real(pReal), dimension(4) :: s - ho = LambertCubeToBall(cu) + 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] -end function cu2ho + 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) + + 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' + +end function om2qu !--------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University -!> @brief convert Rodrigues vector to Euler angles +!> @brief orientation matrix to Euler angles !--------------------------------------------------------------------------------------------------- -pure function ro2eu(ro) result(eu) +pure function om2eu(om) result(eu) + use math, only: & + PI - implicit none - real(pReal), intent(in), dimension(4) :: ro - real(pReal), dimension(3) :: eu - - eu = om2eu(ro2om(ro)) - -end function ro2eu + 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 Euler angles to homochoric +!> @brief convert orientation matrix to axis angle pair !--------------------------------------------------------------------------------------------------- -pure function eu2ho(eu) result(ho) +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(3) :: eu - real(pReal), dimension(3) :: ho + 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 - ho = ax2ho(eu2ax(eu)) - -end function eu2ho +end function om2ax !--------------------------------------------------------------------------------------------------- @@ -924,7 +515,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 +529,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 +543,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 +551,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 +832,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 +949,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 +1081,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 +1089,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 +1123,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 +1133,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 From fa182004470eb67f0e7f8a4a253b4b43787b6dce Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 17 Apr 2019 08:46:03 +0200 Subject: [PATCH 07/29] quaternion is now in separate module avoid long modules with multiple, only loosely related classes --- python/damask/orientation.py | 212 +---------------------------------- python/damask/quaternion.py | 210 ++++++++++++++++++++++++++++++++++ 2 files changed, 213 insertions(+), 209 deletions(-) create mode 100644 python/damask/quaternion.py diff --git a/python/damask/orientation.py b/python/damask/orientation.py index 007aecfd9..349c368fa 100644 --- a/python/damask/orientation.py +++ b/python/damask/orientation.py @@ -3,214 +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 __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 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): - return [self.q]+list(self.p) - - 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 as P #################################################################################################### @@ -488,7 +282,7 @@ class Rotation: 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 + else M + r.asM() * n # noqa add (multiples) of this rotation to average noqa eig, vec = np.linalg.eig(M/N) return Rotation.fromQuaternion(np.real(vec.T[eig.argmax()]),acceptHomomorph = True) diff --git a/python/damask/quaternion.py b/python/damask/quaternion.py new file mode 100644 index 000000000..7c98a6d1c --- /dev/null +++ b/python/damask/quaternion.py @@ -0,0 +1,210 @@ +# -*- coding: UTF-8 no BOM -*- + +import numpy as np + +P = -1 # convention (sed 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. + 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 __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 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): + return [self.q]+list(self.p) + + 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 c7c6627dfb7e1871b3dc334338bd22974862c010 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 17 Apr 2019 08:49:26 +0200 Subject: [PATCH 08/29] don't try to figure out the input --- python/damask/orientation.py | 10 ++-------- 1 file changed, 2 insertions(+), 8 deletions(-) diff --git a/python/damask/orientation.py b/python/damask/orientation.py index 349c368fa..f2302e992 100644 --- a/python/damask/orientation.py +++ b/python/damask/orientation.py @@ -446,13 +446,7 @@ 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 R.shape[0]==4: # transition old (length not stored separately) to new - R = (R[0:3]*R[3]) + R = rodrigues epsilon = 0.0 if self.lattice == 'cubic': @@ -910,7 +904,7 @@ class Orientation: for k in range(2): r.inverse() breaker = self.lattice.symmetry.inFZ(r.asRodrigues()) \ - and (not SST or other.lattice.symmetry.inDisorientationSST(r.asRodrigues())) + and (not SST or other.lattice.symmetry.inDisorientationSST(r.asRodrigues(vector=True))) if breaker: break if breaker: break if breaker: break From 8609eb6eb77ec62c2dec93d994d966cc0abc9b2c Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 17 Apr 2019 09:29:49 +0200 Subject: [PATCH 09/29] improved names and layout --- python/damask/orientation.py | 129 ++++++++++++++++++----------------- 1 file changed, 66 insertions(+), 63 deletions(-) diff --git a/python/damask/orientation.py b/python/damask/orientation.py index f2302e992..c989bc7c0 100644 --- a/python/damask/orientation.py +++ b/python/damask/orientation.py @@ -110,6 +110,11 @@ class Rotation: def misorientation(self,other): """Misorientation""" return self.__class__(other.quaternion*self.quaternion.conjugated()) + + + def average(self,other): + """Calculate the average rotation""" + return Rotation.fromAverage([self,other]) ################################################################################################ @@ -262,9 +267,9 @@ class Rotation: @classmethod - def average(cls, - rotations, - weights = []): + def fromAverage(cls, + rotations, + weights = []): """ Average rotation @@ -272,6 +277,8 @@ class Rotation: 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 not all(isinstance(item, Rotation) for item in rotations): raise TypeError("Only instances of Rotation can be averaged.") @@ -335,56 +342,56 @@ class Symmetry: """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 = [ @@ -447,6 +454,8 @@ class Symmetry: Acta Cryst. (1991). A47, 780-789 """ R = rodrigues + if (len(R) != 3): + raise ValueError('Input is not a Rodriques-Frank vector.\n') epsilon = 0.0 if self.lattice == 'cubic': @@ -967,21 +976,10 @@ class Orientation: @classmethod - def average(cls, - orientations, - weights = []): - """ - Average orientation - - ref: F. Landis Markley, Yang Cheng, John Lucas Crassidis, and Yaakov Oshman. - Averaging Quaternions, - Journal of Guidance, Control, and Dynamics, Vol. 30, No. 4 (2007), pp. 1193-1197. - doi: 10.2514/1.28949 - usage: - a = Orientation(Eulers=np.radians([10, 10, 0]), symmetry='hexagonal') - b = Orientation(Eulers=np.radians([20, 0, 0]), symmetry='hexagonal') - avg = Orientation.average([a,b]) - """ + 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.") @@ -993,7 +991,12 @@ class Orientation: SST = False, # select (o[ther]'s) sym orientation symmetries = True)[2]).rotation) # with lowest misorientation - return Orientation(Rotation.average(closest,weights),ref.lattice) + return Orientation(Rotation.fromAverage(closest,weights),ref.lattice) + + + def average(self,other): + """Calculate the average rotation""" + return Orientation.fromAverage([self,other]) #################################################################################################### From 3f8c16262e6ea32bc18a6937c11bbc991cc6b657 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 17 Apr 2019 09:45:18 +0200 Subject: [PATCH 10/29] fixed import statement previous form only worked in python/damask folder --- python/damask/orientation.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/damask/orientation.py b/python/damask/orientation.py index c989bc7c0..4a6d4e05d 100644 --- a/python/damask/orientation.py +++ b/python/damask/orientation.py @@ -3,8 +3,8 @@ import math import numpy as np from . import Lambert -from quaternion import Quaternion -from quaternion import P as P +from .quaternion import Quaternion +from .quaternion import P #################################################################################################### From 079e683dd181be27bf26c6fd7b4860555793899f Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 17 Apr 2019 10:48:18 +0200 Subject: [PATCH 11/29] disorientation returns Orientation, not Rotation --- python/damask/orientation.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/python/damask/orientation.py b/python/damask/orientation.py index 4a6d4e05d..94ce9c416 100644 --- a/python/damask/orientation.py +++ b/python/damask/orientation.py @@ -899,8 +899,8 @@ 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.') mySymEqs = self.equivalentOrientations() if SST else self.equivalentOrientations([0]) # take all or only first sym operation otherSymEqs = other.equivalentOrientations() @@ -918,7 +918,7 @@ class Orientation: if breaker: break if breaker: break - return (r, i,j, k == 1) if symmetries else r # disorientation ... + return (Orientation(r,self.lattice), i,j, k == 1) if symmetries else Orientation(r,self.lattice)# disorientation ... # ... own sym, other sym, # self-->other: True, self<--other: False From a6e6db05590c4d8ebfcc4e19a4a875a7e7fedbbc Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 17 Apr 2019 12:51:00 +0200 Subject: [PATCH 12/29] more stable/robust conversions --- python/damask/orientation.py | 8 ++++---- src/rotations.f90 | 40 ++++++++---------------------------- 2 files changed, 13 insertions(+), 35 deletions(-) diff --git a/python/damask/orientation.py b/python/damask/orientation.py index 94ce9c416..d9472581a 100644 --- a/python/damask/orientation.py +++ b/python/damask/orientation.py @@ -1077,11 +1077,11 @@ def qu2ax(qu): 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 + 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]): - omega = 2.0 * np.arccos(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] @@ -1126,9 +1126,9 @@ def om2qu(om): """ Orientation matrix to quaternion - The original formulation (direct conversion) had numerical issues + The original formulation (direct conversion) had (numerical?) issues """ - return ax2qu(om2ax(om)) + return eu2qu(om2eu(om)) def om2eu(om): diff --git a/src/rotations.f90 b/src/rotations.f90 index aed4c5773..5279b179c 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -306,12 +306,11 @@ pure function qu2ax(qu) result(ax) real(pReal) :: omega, s - omega = 2.0 * acos(math_clip(qu%w,-1.0_pReal,1.0_pReal)) - ! if the angle equals zero, then we return the rotation axis as [001] - if (dEq0(omega)) then - ax = [ 0.0, 0.0, 1.0, 0.0 ] + 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)) ax = [ qu%x*s, qu%y*s, qu%z*s, omega ] else ax = [ qu%x, qu%y, qu%z, PI ] @@ -397,37 +396,16 @@ end function qu2cu !--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert rotation matrix to a unit quaternion +!> @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) - use prec, only: & - dEq - real(pReal), intent(in), dimension(3,3) :: om - type(quaternion) :: qu - - real(pReal), dimension(4) :: qu_A - real(pReal), dimension(4) :: s + real(pReal), intent(in), dimension(3,3) :: om + type(quaternion) :: qu - 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) - - 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' + qu = eu2qu(om2eu(om)) end function om2qu From fdd5b93e7c28ab892d4dc5cc9ed418e30d5f18fd Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Wed, 17 Apr 2019 15:19:41 +0200 Subject: [PATCH 13/29] avoid FPE exceptions --- python/damask/orientation.py | 8 ++++---- src/rotations.f90 | 7 ++++--- 2 files changed, 8 insertions(+), 7 deletions(-) diff --git a/python/damask/orientation.py b/python/damask/orientation.py index d9472581a..89a4638ca 100644 --- a/python/damask/orientation.py +++ b/python/damask/orientation.py @@ -1132,14 +1132,14 @@ def om2qu(om): 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 - else: + """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) diff --git a/src/rotations.f90 b/src/rotations.f90 index 5279b179c..69529ed24 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -422,14 +422,15 @@ pure function om2eu(om) result(eu) 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 + 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 + where(eu<0.0_pReal) eu = mod(eu+2.0_pReal*PI,[2.0_pReal*PI,PI,2.0_pReal*PI]) end function om2eu From 481aac795254e447680ad377ed31a42715c08734 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 18 Apr 2019 13:46:58 +0200 Subject: [PATCH 14/29] misorientation calculation did not always follow convention - return homomorphed quaternion to ensure that real component is > 0 (equivalent to have axis angle limited to [0, pi] --- python/damask/orientation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/orientation.py b/python/damask/orientation.py index 89a4638ca..1f0dc521a 100644 --- a/python/damask/orientation.py +++ b/python/damask/orientation.py @@ -109,7 +109,7 @@ class Rotation: def misorientation(self,other): """Misorientation""" - return self.__class__(other.quaternion*self.quaternion.conjugated()) + return self.__class__((other.quaternion*self.quaternion.conjugated()).homomorph()) def average(self,other): From 71b00363239634477ab4689f0778b65f7bce0807 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 18 Apr 2019 14:26:53 +0200 Subject: [PATCH 15/29] return rotation as misorientation --- python/damask/orientation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/orientation.py b/python/damask/orientation.py index 1f0dc521a..c0c817ca6 100644 --- a/python/damask/orientation.py +++ b/python/damask/orientation.py @@ -918,7 +918,7 @@ class Orientation: if breaker: break if breaker: break - return (Orientation(r,self.lattice), i,j, k == 1) if symmetries else Orientation(r,self.lattice)# disorientation ... + return (r, i,j, k == 1) if symmetries else r # disorientation ... # ... own sym, other sym, # self-->other: True, self<--other: False From c702441bcc9f3a61c628e79ccd7ae31aa7202945 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 18 Apr 2019 14:27:07 +0200 Subject: [PATCH 16/29] now testing averaging and misorientation calculation --- .gitlab-ci.yml | 8 ++++++++ PRIVATE | 2 +- 2 files changed, 9 insertions(+), 1 deletion(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 45fb4e4f4..fa9cab906 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/PRIVATE b/PRIVATE index c7bc54a26..293d04c73 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit c7bc54a26c8b6ed404aabec4653227e93fa028e2 +Subproject commit 293d04c7372ddda2f79d4487795b983b0b552f0f From e8464814ae295d430443cb030efbdf8385871258 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 18 Apr 2019 14:46:12 +0200 Subject: [PATCH 17/29] explicit is better than implicit - if the quaternion does not follow the convention, use "fromQuaternion(acceptHomomorph=True)" --- python/damask/orientation.py | 1 - 1 file changed, 1 deletion(-) diff --git a/python/damask/orientation.py b/python/damask/orientation.py index c0c817ca6..b34875a1f 100644 --- a/python/damask/orientation.py +++ b/python/damask/orientation.py @@ -44,7 +44,6 @@ class Rotation: self.quaternion = quaternion.copy() else: self.quaternion = Quaternion(q=quaternion[0],p=quaternion[1:4]) - self.quaternion.homomorph() # ToDo: Needed? def __repr__(self): """Value in selected representation""" From ce1bb589c7992927b553892d81bf518ab7abb365 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 18 Apr 2019 15:30:54 +0200 Subject: [PATCH 18/29] complicated use statements confuse CMake probably we should adopt a similar "use" policy in Fortran as recommended for Python "import" and have all imports at the beginning https://www.python.org/dev/peps/pep-0008/ --- src/plastic_disloUCLA.f90 | 10 +--------- src/plastic_dislotwin.f90 | 10 +--------- src/plastic_isotropic.f90 | 13 +------------ src/plastic_kinematichardening.f90 | 13 +------------ src/plastic_none.f90 | 8 +------- src/plastic_nonlocal.f90 | 10 +--------- src/plastic_phenopowerlaw.f90 | 10 +--------- 7 files changed, 7 insertions(+), 67 deletions(-) diff --git a/src/plastic_disloUCLA.f90 b/src/plastic_disloUCLA.f90 index 88aa27432..ccb37be48 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 cb13265b4..ed6638c4f 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 05a31ab75..45bc136db 100644 --- a/src/plastic_isotropic.f90 +++ b/src/plastic_isotropic.f90 @@ -91,18 +91,7 @@ subroutine plastic_isotropic_init debug_levelBasic use IO, only: & IO_error - use material, only: & -#ifdef DEBUG - phasememberAt, & -#endif - phase_plasticity, & - phase_plasticityInstance, & - phase_Noutput, & - material_allocatePlasticState, & - PLASTICITY_ISOTROPIC_label, & - PLASTICITY_ISOTROPIC_ID, & - material_phase, & - plasticState + use material use config, only: & config_phase use lattice diff --git a/src/plastic_kinematichardening.f90 b/src/plastic_kinematichardening.f90 index 0a4a6b3bc..6890c7006 100644 --- a/src/plastic_kinematichardening.f90 +++ b/src/plastic_kinematichardening.f90 @@ -112,18 +112,7 @@ subroutine plastic_kinehardening_init math_expand use IO, only: & IO_error - use material, only: & -#ifdef DEBUG - phasememberAt, & -#endif - phase_plasticity, & - phase_plasticityInstance, & - phase_Noutput, & - material_allocatePlasticState, & - PLASTICITY_kinehardening_label, & - PLASTICITY_kinehardening_ID, & - material_phase, & - plasticState + use material use config, only: & config_phase use lattice diff --git a/src/plastic_none.f90 b/src/plastic_none.f90 index f8a64b55b..fb19a9968 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 implicit none integer :: & diff --git a/src/plastic_nonlocal.f90 b/src/plastic_nonlocal.f90 index f5f48ed11..19f6b38cc 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 4124856d1..8e7ee82d0 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 From 7570b187e75536d23ee90219d3a57e8641ebffd0 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 18 Apr 2019 16:26:05 +0200 Subject: [PATCH 19/29] do explicit homomorphing only in cases when the underlying math cannot guarantee that the resulting quaternion follows the convention --- python/damask/orientation.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/python/damask/orientation.py b/python/damask/orientation.py index b34875a1f..1cf570583 100644 --- a/python/damask/orientation.py +++ b/python/damask/orientation.py @@ -45,6 +45,7 @@ class Rotation: else: self.quaternion = Quaternion(q=quaternion[0],p=quaternion[1:4]) + def __repr__(self): """Value in selected representation""" return '\n'.join([ @@ -60,7 +61,7 @@ class Rotation: 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).asArray()) + return self.__class__((self.quaternion * other.quaternion).homomorph()) elif isinstance(other, np.ndarray): if other.shape == (3,): # rotate a single (3)-vector ( x, y, z) = self.quaternion.p @@ -108,7 +109,7 @@ class Rotation: def misorientation(self,other): """Misorientation""" - return self.__class__((other.quaternion*self.quaternion.conjugated()).homomorph()) + return other*self.inversed() def average(self,other): @@ -1192,7 +1193,7 @@ def eu2qu(eu): -P*sPhi*np.cos(ee[0]-ee[2]), -P*sPhi*np.sin(ee[0]-ee[2]), -P*cPhi*np.sin(ee[0]+ee[2]) ]) - #if qu[0] < 0.0: qu.homomorph() !ToDo: Check with original + if qu[0] < 0.0: qu*=-1 return qu From 1904ac211e9a173c8cb7c3fbffff04ce2db85862 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 18 Apr 2019 16:44:39 +0200 Subject: [PATCH 20/29] test working for more explicit rotation class --- PRIVATE | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PRIVATE b/PRIVATE index 293d04c73..f6171a748 160000 --- a/PRIVATE +++ b/PRIVATE @@ -1 +1 @@ -Subproject commit 293d04c7372ddda2f79d4487795b983b0b552f0f +Subproject commit f6171a748e51b994db27c2cc74cc0168b7aea93f From c4f56703a405da3c7609190bfd5af1d933c17c36 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Thu, 18 Apr 2019 17:04:00 +0200 Subject: [PATCH 21/29] proper definition of homomorph --- python/damask/orientation.py | 4 +++- python/damask/quaternion.py | 5 ++--- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/python/damask/orientation.py b/python/damask/orientation.py index 1cf570583..626acc0f2 100644 --- a/python/damask/orientation.py +++ b/python/damask/orientation.py @@ -61,7 +61,9 @@ class Rotation: 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).homomorph()) + qu = self.quaternion * other.quaternion + if qu.q < 0: qu.homomorph() + return self.__class__(qu) elif isinstance(other, np.ndarray): if other.shape == (3,): # rotate a single (3)-vector ( x, y, z) = self.quaternion.p diff --git a/python/damask/quaternion.py b/python/damask/quaternion.py index 7c98a6d1c..1d476e598 100644 --- a/python/damask/quaternion.py +++ b/python/damask/quaternion.py @@ -201,9 +201,8 @@ class Quaternion: def homomorph(self): - if self.q < 0.0: - self.q = -self.q - self.p = -self.p + self.q = -self.q + self.p = -self.p return self def homomorphed(self): From 4183eeaeca1b15a9ec5f66cbb65b071b5baf84ba Mon Sep 17 00:00:00 2001 From: Philip Eisenlohr Date: Thu, 18 Apr 2019 18:05:46 -0400 Subject: [PATCH 22/29] [skip ci] disorientation may return full-fledged orientation obj and aux info, added rotation.standardize(d) --- python/damask/orientation.py | 11 ++++++++++- 1 file changed, 10 insertions(+), 1 deletion(-) diff --git a/python/damask/orientation.py b/python/damask/orientation.py index 626acc0f2..26e926357 100644 --- a/python/damask/orientation.py +++ b/python/damask/orientation.py @@ -109,6 +109,15 @@ class Rotation: return self.__class__(self.quaternion.conjugated()) + def standardize(self): + """Ensure quaternion representation with positive q""" + if self.quaternion.q < 0.0: self.quaternion.homomorph() + + def standardized(self): + """Ensure quaternion representation with positive q""" + return self.__class__(self.quaternion.homomorphed() if self.quaternion.q < 0.0 else self.quaternion) + + def misorientation(self,other): """Misorientation""" return other*self.inversed() @@ -920,7 +929,7 @@ class Orientation: if breaker: break if breaker: break - return (r, i,j, k == 1) if symmetries else r # disorientation ... + return (Orientation(r,self.symmetry), i,j, k == 1) if symmetries else r # disorientation ... # ... own sym, other sym, # self-->other: True, self<--other: False From 25518df80c67f270b3cfbab57e299060d2bcb76b Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 19 Apr 2019 00:42:20 +0200 Subject: [PATCH 23/29] fixes to quaternion class - always store data as float (solves issues with truediv) - fixed in-place multiplication and division - consistent handling of "..ed()" functions --- python/damask/quaternion.py | 33 +++++++++++++++++++-------------- 1 file changed, 19 insertions(+), 14 deletions(-) diff --git a/python/damask/quaternion.py b/python/damask/quaternion.py index 1d476e598..1bb6f5f0d 100644 --- a/python/damask/quaternion.py +++ b/python/damask/quaternion.py @@ -10,15 +10,15 @@ class Quaternion: 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}. + 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 = q - self.p = np.array(p) + self.q = float(q) + self.p = np.array(p,dtype=float) def __copy__(self): @@ -102,7 +102,8 @@ class Quaternion: 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 + self.q *= other + self.p *= other return self else: return NotImplemented @@ -115,9 +116,8 @@ class Quaternion: return self.__class__(q=self.q * s, p=self.p * s) elif isinstance(other, (int, float)): - self.q /= other - self.p /= other - return self + return self.__class__(q=self.q / other, + p=self.p / other) else: return NotImplemented @@ -129,6 +129,7 @@ class Quaternion: return self elif isinstance(other, (int, float)): self.q /= other + self.p /= other return self else: return NotImplemented @@ -149,6 +150,7 @@ class Quaternion: 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 @@ -179,31 +181,34 @@ class Quaternion: 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 (no return value)""" d = self.magnitude() - if d > 0.0: - self.q /= d - self.p /= d - return self + if d > 0.0: self /= d def normalized(self): + """Returns normalized copy""" return self.copy().normalize() - + def conjugate(self): + """Conjugates in-place (no return value)""" self.p = -self.p - return self def conjugated(self): + """Returns conjugated copy""" return self.copy().conjugate() def homomorph(self): + """Homomorphs in-place (no return value)""" self.q = -self.q self.p = -self.p - return self def homomorphed(self): + """Returns homomorphed copy""" return self.copy().homomorph() From 9a43c2e4c5102101026dfc2e973e167c0b5509b4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 19 Apr 2019 01:05:48 +0200 Subject: [PATCH 24/29] fixed return values - homomorph, standardize, etc. are silent in-place operations (return None) - homomorphed, standardized, etc. are out-of place operations that report --- python/damask/orientation.py | 23 ++++++++++++++++------- python/damask/quaternion.py | 12 +++++++++--- 2 files changed, 25 insertions(+), 10 deletions(-) diff --git a/python/damask/orientation.py b/python/damask/orientation.py index 26e926357..dc3389977 100644 --- a/python/damask/orientation.py +++ b/python/damask/orientation.py @@ -44,6 +44,12 @@ class Rotation: self.quaternion = quaternion.copy() else: self.quaternion = Quaternion(q=quaternion[0],p=quaternion[1:4]) + + def __copy__(self): + """Copy""" + return self.__class__(self.quaternion) + + copy = __copy__ def __repr__(self): @@ -61,9 +67,9 @@ class Rotation: Rotation: Details needed (active/passive), rotation of (3,3,3,3)-matrix should be considered """ if isinstance(other, Rotation): # rotate a rotation - qu = self.quaternion * other.quaternion - if qu.q < 0: qu.homomorph() - return self.__class__(qu) + qu = self.__class__(self.quaternion * other.quaternion) + qu.standardize() + return qu elif isinstance(other, np.ndarray): if other.shape == (3,): # rotate a single (3)-vector ( x, y, z) = self.quaternion.p @@ -102,11 +108,12 @@ class Rotation: def inverse(self): """In-place inverse rotation/backward rotation""" self.quaternion.conjugate() - return self def inversed(self): """Inverse rotation/backward rotation""" - return self.__class__(self.quaternion.conjugated()) + c = self.copy() + c.inverse() + return c def standardize(self): @@ -115,7 +122,9 @@ class Rotation: def standardized(self): """Ensure quaternion representation with positive q""" - return self.__class__(self.quaternion.homomorphed() if self.quaternion.q < 0.0 else self.quaternion) + c = self.copy() + c.standardize() + return c def misorientation(self,other): @@ -929,7 +938,7 @@ class Orientation: if breaker: break if breaker: break - return (Orientation(r,self.symmetry), i,j, k == 1) if symmetries else r # disorientation ... + return (Orientation(r,self.lattice), i,j, k == 1) if symmetries else r # disorientation ... # ... own sym, other sym, # self-->other: True, self<--other: False diff --git a/python/damask/quaternion.py b/python/damask/quaternion.py index 1bb6f5f0d..69f4e6ef0 100644 --- a/python/damask/quaternion.py +++ b/python/damask/quaternion.py @@ -192,7 +192,9 @@ class Quaternion: def normalized(self): """Returns normalized copy""" - return self.copy().normalize() + c = self.copy() + c.normalize() + return c def conjugate(self): @@ -201,7 +203,9 @@ class Quaternion: def conjugated(self): """Returns conjugated copy""" - return self.copy().conjugate() + c = self.copy() + c.conjugate() + return c def homomorph(self): @@ -211,4 +215,6 @@ class Quaternion: def homomorphed(self): """Returns homomorphed copy""" - return self.copy().homomorph() + c = self.copy() + c.homomorph() + return c From feca9fe0a0c51aab090e42fc66002b5699a20226 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 19 Apr 2019 08:04:04 +0200 Subject: [PATCH 25/29] Again changes to the return values. A python function with no return value returns 'None' --- python/damask/orientation.py | 16 +++++++--------- python/damask/quaternion.py | 23 ++++++++++------------- 2 files changed, 17 insertions(+), 22 deletions(-) diff --git a/python/damask/orientation.py b/python/damask/orientation.py index dc3389977..77fa3665c 100644 --- a/python/damask/orientation.py +++ b/python/damask/orientation.py @@ -108,23 +108,21 @@ class Rotation: def inverse(self): """In-place inverse rotation/backward rotation""" self.quaternion.conjugate() + return self def inversed(self): """Inverse rotation/backward rotation""" - c = self.copy() - c.inverse() - return c + return self.copy.inverse() def standardize(self): - """Ensure quaternion representation with positive q""" - if self.quaternion.q < 0.0: self.quaternion.homomorph() + """In-place quaternion representation with positive q""" + if self.quaternion.q < 0.0: self.quaternion.homomorph() + return self def standardized(self): - """Ensure quaternion representation with positive q""" - c = self.copy() - c.standardize() - return c + """Quaternion representation with positive q""" + return self.copy.standardize() def misorientation(self,other): diff --git a/python/damask/quaternion.py b/python/damask/quaternion.py index 69f4e6ef0..9a361cb3d 100644 --- a/python/damask/quaternion.py +++ b/python/damask/quaternion.py @@ -2,7 +2,7 @@ import numpy as np -P = -1 # convention (sed DOI:10.1088/0965-0393/23/8/083501) +P = -1 # convention (see DOI:10.1088/0965-0393/23/8/083501) #################################################################################################### class Quaternion: @@ -186,35 +186,32 @@ class Quaternion: def normalize(self): - """Normalizes in-place (no return value)""" + """Normalizes in-place""" d = self.magnitude() if d > 0.0: self /= d + return self def normalized(self): """Returns normalized copy""" - c = self.copy() - c.normalize() - return c + return self.copy.normalize() def conjugate(self): - """Conjugates in-place (no return value)""" + """Conjugates in-place""" self.p = -self.p + return self def conjugated(self): """Returns conjugated copy""" - c = self.copy() - c.conjugate() - return c + return self.copy.conjugate() def homomorph(self): - """Homomorphs in-place (no return value)""" + """Homomorphs in-place""" self.q = -self.q self.p = -self.p + return self def homomorphed(self): """Returns homomorphed copy""" - c = self.copy() - c.homomorph() - return c + return self.copy.homomorph() From b97f10b6ff7322c97c07dc5c3e0e5fcdc7296be4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 19 Apr 2019 08:17:02 +0200 Subject: [PATCH 26/29] forgotten changes in last commit + fromRandom --- python/damask/orientation.py | 42 ++++++++++++++++++++++++------------ python/damask/quaternion.py | 6 +++--- 2 files changed, 31 insertions(+), 17 deletions(-) diff --git a/python/damask/orientation.py b/python/damask/orientation.py index 77fa3665c..e584415b0 100644 --- a/python/damask/orientation.py +++ b/python/damask/orientation.py @@ -67,9 +67,7 @@ class Rotation: Rotation: Details needed (active/passive), rotation of (3,3,3,3)-matrix should be considered """ if isinstance(other, Rotation): # rotate a rotation - qu = self.__class__(self.quaternion * other.quaternion) - qu.standardize() - return qu + 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 @@ -112,7 +110,7 @@ class Rotation: def inversed(self): """Inverse rotation/backward rotation""" - return self.copy.inverse() + return self.copy().inverse() def standardize(self): @@ -122,7 +120,7 @@ class Rotation: def standardized(self): """Quaternion representation with positive q""" - return self.copy.standardize() + return self.copy().standardize() def misorientation(self,other): @@ -310,7 +308,20 @@ class Rotation: else M + r.asM() * n # noqa add (multiples) of this rotation to average noqa eig, vec = np.linalg.eig(M/N) - return Rotation.fromQuaternion(np.real(vec.T[eig.argmax()]),acceptHomomorph = True) + 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) + # ****************************************************************************************** @@ -433,15 +444,18 @@ class Symmetry: 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. """ - if np.any(R == np.inf): return False + if (len(rodrigues) != 3): + raise ValueError('Input is not a Rodriques-Frank vector.\n') + + if np.any(rodrigues == np.inf): return False - Rabs = abs(R[0:3]*R[3]) + Rabs = abs(rodrigues) if self.lattice == 'cubic': return math.sqrt(2.0)-1.0 >= Rabs[0] \ @@ -471,10 +485,10 @@ class Symmetry: Representation of Orientation and Disorientation Data for Cubic, Hexagonal, Tetragonal and Orthorhombic Crystals Acta Cryst. (1991). A47, 780-789 """ - R = rodrigues - if (len(R) != 3): + if (len(rodrigues) != 3): raise ValueError('Input is not a Rodriques-Frank vector.\n') - + R = rodrigues + epsilon = 0.0 if self.lattice == 'cubic': return R[0] >= R[1]+epsilon and R[1] >= R[2]+epsilon and R[2] >= epsilon @@ -930,7 +944,7 @@ class Orientation: r = b*aInv for k in range(2): r.inverse() - breaker = self.lattice.symmetry.inFZ(r.asRodrigues()) \ + 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 @@ -942,7 +956,7 @@ class Orientation: def inFZ(self): - return self.lattice.symmetry.inFZ(self.rotation.asRodrigues()) + return self.lattice.symmetry.inFZ(self.rotation.asRodrigues(vector=True)) def equivalentOrientations(self,members=[]): """List of orientations which are symmetrically equivalent""" diff --git a/python/damask/quaternion.py b/python/damask/quaternion.py index 9a361cb3d..3a5fa9a89 100644 --- a/python/damask/quaternion.py +++ b/python/damask/quaternion.py @@ -193,7 +193,7 @@ class Quaternion: def normalized(self): """Returns normalized copy""" - return self.copy.normalize() + return self.copy().normalize() def conjugate(self): @@ -203,7 +203,7 @@ class Quaternion: def conjugated(self): """Returns conjugated copy""" - return self.copy.conjugate() + return self.copy().conjugate() def homomorph(self): @@ -214,4 +214,4 @@ class Quaternion: def homomorphed(self): """Returns homomorphed copy""" - return self.copy.homomorph() + return self.copy().homomorph() From 13307307174276952a207da4d3c8660f8e3ef2c4 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Fri, 19 Apr 2019 09:04:54 +0200 Subject: [PATCH 27/29] inFZ needs 'proper' Rodriques--Franck vector --- python/damask/orientation.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/damask/orientation.py b/python/damask/orientation.py index e584415b0..00550d986 100644 --- a/python/damask/orientation.py +++ b/python/damask/orientation.py @@ -976,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) From 9b0d3def1d2fafce03d7d833224e027549bdfd09 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 5 May 2019 22:35:12 +0200 Subject: [PATCH 28/29] small changes suggested by Franz --- python/damask/orientation.py | 2 +- python/damask/quaternion.py | 9 +++------ 2 files changed, 4 insertions(+), 7 deletions(-) diff --git a/python/damask/orientation.py b/python/damask/orientation.py index 51b6f849f..b19037b0c 100644 --- a/python/damask/orientation.py +++ b/python/damask/orientation.py @@ -83,7 +83,7 @@ class Rotation: ]) 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): + elif other.shape == (3,3,3,3,): raise NotImplementedError else: return NotImplemented diff --git a/python/damask/quaternion.py b/python/damask/quaternion.py index 3a5fa9a89..fd681fe9b 100644 --- a/python/damask/quaternion.py +++ b/python/damask/quaternion.py @@ -78,10 +78,8 @@ class Quaternion: return NotImplemented def __neg__(self): - """Unary positive operator""" - self.q *= -1.0 - self.p *= -1.0 - return self + """Unary negative operator""" + return self * -1.0 def __mul__(self, other): @@ -208,8 +206,7 @@ class Quaternion: def homomorph(self): """Homomorphs in-place""" - self.q = -self.q - self.p = -self.p + self *= -1.0 return self def homomorphed(self): From 8dd2a131fefcc25a2a62ba58cefdf07987fffbc1 Mon Sep 17 00:00:00 2001 From: Test User Date: Mon, 6 May 2019 14:30:22 +0200 Subject: [PATCH 29/29] [skip ci] updated version information after successful test of v2.0.3-229-g90fa809f --- VERSION | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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