Merge branch 'orientationClass_with_negative_P' into development

lib/damask/orientation.py

@@ -62,8 +62,7 @@ class Quaternion:
 
     def __copy__(self):
       """Copy"""
-      Q = Quaternion(q=self.q,p=self.p)
+      return self.__class__(q=self.q,p=self.p.copy())
-      return Q
 
     copy = __copy__
 
@@ -73,15 +72,13 @@ class Quaternion:
 
     def __pow__(self, exponent):
       """Power"""
-      Q = Quaternion()
       omega = math.acos(self.q)
-      Q.q =          math.cos(exponent*omega)
+      return self.__class__(q=         math.cos(exponent*omega),
-      Q.p = self.p * math.sin(exponent*omega)/math.sin(omega)
+                            p=self.p * math.sin(exponent*omega)/math.sin(omega))
-      return Q
 
     def __ipow__(self, exponent):
       """In-place power"""
-      omega = math.acos(self.q[0])
+      omega = math.acos(self.q)
       self.q  = math.cos(exponent*omega)
       self.p *= math.sin(exponent*omega)/math.sin(omega)
       return self
@@ -91,21 +88,25 @@ class Quaternion:
       # Rowenhorst_etal2015 MSMSE: value of P is selected as -1
       P = -1.0
       try:                                                          # quaternion
-          Q = Quaternion()
+          return self.__class__(q=self.q*other.q - np.dot(self.p,other.p),
-          Q.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))
-          Q.p = self.q*other.p + other.q*self.p + P * np.cross(self.p,other.p)
-          return Q
       except: pass
       try:                                                         # vector (perform passive rotation)
-          return  (self.q*self.q - np.dot(self.p,self.p)) * np.array(other[:3]) \
+          ( x, y, z) = self.p
-                 + 2.0*np.dot(self.p,other[:3])           * self.p \
+          (Vx,Vy,Vz) = other[0:3]
-                 + 2.0*P*self.q                           * np.cross(self.p,other[:3])
+          A = self.q*self.q - np.dot(self.p,self.p)
+          B = 2.0 * (x*Vx + y*Vy + z*Vz)
+          C = 2.0 * P*self.q
+
+          return np.array([
+            A*Vx + B*x + C*(y*Vz - z*Vy),
+            A*Vy + B*y + C*(z*Vx - x*Vz),
+            A*Vz + B*z + C*(x*Vy - y*Vx),
+            ])
       except: pass
       try:                                                        # scalar
-          Q = self.copy()
+          return self.__class__(q=self.q*other,
-          Q.q *= other
+                                p=self.p*other)
-          Q.p *= other
-          return Q
       except:
           return self.copy()
 
@@ -122,9 +123,8 @@ class Quaternion:
     def __div__(self, other):
       """Division"""
       if isinstance(other, (int,float)):
-        q = self.q / other
+        return self.__class__(q=self.q / other,
-        p = self.p / other
+                              p=self.p / other)
-        return self.__class__(q=q,p=p)
       else:
           return NotImplemented
 
@@ -138,9 +138,8 @@ class Quaternion:
     def __add__(self, other):
       """Addition"""
       if isinstance(other, Quaternion):
-        q = self.q + other.q
+        return self.__class__(q=self.q + other.q,
-        p = self.p + other.p
+                              p=self.p + other.p)
-        return self.__class__(q=q,p=p)
       else:
           return NotImplemented
 
@@ -154,12 +153,10 @@ class Quaternion:
     def __sub__(self, other):
       """Subtraction"""
       if isinstance(other, Quaternion):
-          Q = self.copy()
+          return self.__class__(q=self.q - other.q,
-          Q.q -= other.q
+                                p=self.p - other.p)
-          Q.p -= other.p
-          return Q
       else:
-          return self.copy()
+          return NotImplemented
 
     def __isub__(self, other):
       """In-place subtraction"""
@@ -190,7 +187,8 @@ class Quaternion:
 
     def __cmp__(self,other):
       """Linear ordering"""
-      return (self.Rodrigues()>other.Rodrigues()) - (self.Rodrigues()<other.Rodrigues())
+      return (1 if np.linalg.norm(self.asRodrigues()) > np.linalg.norm(other.asRodrigues()) else 0) \
+           - (1 if np.linalg.norm(self.asRodrigues()) < np.linalg.norm(other.asRodrigues()) else 0)
 
     def magnitude_squared(self):
       return self.q ** 2 + np.dot(self.p,self.p)
@@ -203,7 +201,8 @@ class Quaternion:
     def normalize(self):
       d = self.magnitude()
       if d > 0.0:
-          self /= d
+          self.q /= d
+          self.p /= d
       return self
 
     def conjugate(self):
@@ -322,9 +321,11 @@ class Quaternion:
     @classmethod
     def fromRodrigues(cls, rodrigues):
       if not isinstance(rodrigues, np.ndarray): rodrigues = np.array(rodrigues)
-      halfangle = math.atan(np.linalg.norm(rodrigues))
+      norm = np.linalg.norm(rodrigues)
+      halfangle = math.atan(norm)
+      s = math.sin(halfangle)
       c = math.cos(halfangle)
-      return cls(q=c,p=rodrigues/c)
+      return cls(q=c,p=s*rodrigues/norm)
 
 
     @classmethod
@@ -860,7 +861,7 @@ class Orientation:
       M = closest.quaternion.asM() * n if i == 0 else M + closest.quaternion.asM() * n              # noqa add (multiples) of this orientation to average noqa
     eig, vec = np.linalg.eig(M/N)
 
-    return Orientation(quaternion = Quaternion(quatArray = np.real(vec.T[eig.argmax()])),
+    return Orientation(quaternion = Quaternion(quat = np.real(vec.T[eig.argmax()])),
                        symmetry = reference.symmetry.lattice)
 
 

processing/post/addGrainID.py

@@ -200,9 +200,9 @@ for name in filenames:
       if gID != -1 and gID not in alreadyChecked:                                                   # indexed point belonging to a grain not yet tested?
         alreadyChecked[gID] = True                                                                  # remember not to check again
         disorientation = o.disorientation(orientations[gID],SST = False)[0]                         # compare against other orientation
-        if disorientation.quaternion.w >  cos_disorientation:                                       # within threshold ...
+        if disorientation.quaternion.q >  cos_disorientation:                                       # within threshold ...
           candidates.append(gID)                                                                    # remember as potential candidate
-          if disorientation.quaternion.w >= bestDisorientation.w:                                   # ... and better than current best? 
+          if disorientation.quaternion.q >= bestDisorientation.q:                                   # ... and better than current best? 
             matched = True
             matchedID = gID                                                                         # remember that grain
             bestDisorientation = disorientation.quaternion

processing/pre/geom_addPrimitive.py

@@ -64,11 +64,11 @@ if options.dimension is None:
   parser.error('no dimension specified.')   
 if options.angleaxis is not None:
   options.angleaxis = list(map(float,options.angleaxis))
-  rotation = damask.Quaternion().fromAngleAxis(np.radians(options.angleaxis[0]) if options.degrees else options.angleaxis[0],
+  rotation = damask.Quaternion.fromAngleAxis(np.radians(options.angleaxis[0]) if options.degrees else options.angleaxis[0],
-                                               options.angleaxis[1:4])
+                                             options.angleaxis[1:4])
 elif options.quaternion is not None:
-  options.quaternion = map(float,options.quaternion)
+  options.quaternion = list(map(float,options.quaternion))
-  rotation = damask.Quaternion(options.quaternion)
+  rotation = damask.Quaternion(quat=options.quaternion)
 else:
   rotation = damask.Quaternion()
 

processing/pre/geom_fromTable.py

@@ -258,7 +258,7 @@ for name in filenames:
 
             if len(grains) > 0:                                                                     # check immediate neighborhood first
               cos_disorientations = np.array([o.disorientation(orientations[grainID],
-                                                               SST = False)[0].quaternion.w \
+                                                               SST = False)[0].quaternion.q \
                                               for grainID in grains])                               # store disorientation per grainID
               closest_grain = np.argmax(cos_disorientations)                                        # grain among grains with closest orientation to myself
               match = 'local'
@@ -269,7 +269,7 @@ for name in filenames:
 
               if len(grains) > 0:
                 cos_disorientations = np.array([o.disorientation(orientations[grainID],
-                                                                 SST = False)[0].quaternion.w \
+                                                                 SST = False)[0].quaternion.q \
                                                 for grainID in grains])                             # store disorientation per grainID
                 closest_grain = np.argmax(cos_disorientations)                                      # grain among grains with closest orientation to myself
                 match = 'global'