From 48b0307fab9c932b7925bd3fa852dca67bdf1092 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 12 Feb 2019 00:11:22 +0100 Subject: [PATCH] using new rotation class cannot rotate matrices (3,3) and (3,3,3,3) at the moment --- processing/post/addOrientations.py | 47 +++++++---------------------- processing/pre/geom_addPrimitive.py | 11 +++---- python/damask/orientation.py | 32 ++++++++++++++++++-- 3 files changed, 45 insertions(+), 45 deletions(-) diff --git a/processing/post/addOrientations.py b/processing/post/addOrientations.py index a33f96b91..ec824c88c 100755 --- a/processing/post/addOrientations.py +++ b/processing/post/addOrientations.py @@ -9,31 +9,6 @@ import damask scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptID = ' '.join([scriptName,damask.version]) -# -------------------------------------------------------------------- -# convention conformity checks -# -------------------------------------------------------------------- - -def check_Eulers(eulers): - if np.any(eulers < 0.0) or np.any(eulers > 2.0*np.pi) or eulers[1] > np.pi: # Euler angles within valid range? - raise ValueError('Euler angles outside of [0..2π],[0..π],[0..2π].\n{} {} {}.'.format(*eulers)) - return eulers - -def check_quaternion(q): - if q[0] < 0.0: # positive first quaternion component? - raise ValueError('quaternion has negative first component.\n{}'.format(q[0])) - if not np.isclose(np.linalg.norm(q), 1.0): # unit quaternion? - raise ValueError('quaternion is not of unit length.\n{} {} {} {}'.format(*q)) - return q - -def check_matrix(M): - if not np.isclose(np.linalg.det(M),1.0): # proper rotation? - raise ValueError('matrix is not a proper rotation.\n{}'.format(M)) - if not np.isclose(np.dot(M[0],M[1]), 0.0) \ - or not np.isclose(np.dot(M[1],M[2]), 0.0) \ - or not np.isclose(np.dot(M[2],M[0]), 0.0): # all orthogonal? - raise ValueError('matrix is not orthogonal.\n{}'.format(M)) - return M - # -------------------------------------------------------------------- # MAIN # -------------------------------------------------------------------- @@ -133,9 +108,8 @@ if np.sum(input) != 1: parser.error('needs exactly one input format.') (options.quaternion,4,'quaternion'), ][np.where(input)[0][0]] # select input label that was requested -toRadians = np.pi/180.0 if options.degrees else 1.0 # rescale degrees to radians -r = damask.Quaternion.fromAngleAxis(toRadians*options.crystalrotation[0],options.crystalrotation[1:]) # crystal frame rotation -R = damask.Quaternion.fromAngleAxis(toRadians*options. labrotation[0],options. labrotation[1:]) # lab frame rotation +r = damask.Rotation.fromAngleAxis(np.array(options.crystalrotation),options.degrees) # crystal frame rotation +R = damask.Rotation.fromAngleAxis(np.array(options.labrotation),options.degrees) # lab frame rotation # --- loop over input files ------------------------------------------------------------------------ @@ -179,23 +153,24 @@ for name in filenames: outputAlive = True while outputAlive and table.data_read(): # read next data line of ASCII table if inputtype == 'eulers': + o = damask.Rotation.fromEulers(np.array(list(map(float,table.data[column:column+3]))),options.degrees) - o = damask.Orientation(Eulers = check_Eulers(np.array(list(map(float,table.data[column:column+3])))*toRadians)) elif inputtype == 'rodrigues': - o = damask.Orientation(Rodrigues = np.array(list(map(float,table.data[column:column+3])))) - elif inputtype == 'matrix': + o = damask.Rotation.fromRodrigues(np.array(list(map(float,table.data[column:column+3])))) - o = damask.Orientation(matrix = check_matrix(np.array(list(map(float,table.data[column:column+9]))).reshape(3,3))) + elif inputtype == 'matrix': + o = damask.Rotation.fromMatrix(np.array(list(map(float,table.data[column:column+9]))).reshape(3,3)) + elif inputtype == 'frame': M = np.array(list(map(float,table.data[column[0]:column[0]+3] + \ table.data[column[1]:column[1]+3] + \ table.data[column[2]:column[2]+3]))).reshape(3,3).T - o = damask.Orientation(matrix = check_matrix(M/np.linalg.norm(M,axis=0))) - elif inputtype == 'quaternion': + o = damask.Rotation.fromMatrix(M/np.linalg.norm(M,axis=0)) - o = damask.Orientation(quaternion = check_quaternion(np.array(list(map(float,table.data[column:column+4]))))) + elif inputtype == 'quaternion': + o = damask.Rotation.fromQuaternion(np.array(list(map(float,table.data[column:column+4])))) - o.quaternion = r*o.quaternion*R # apply additional lab and crystal frame rotations + o= r*o*R # apply additional lab and crystal frame rotations for output in options.output: if output == 'quaternion': table.data_append(o.asQuaternion()) diff --git a/processing/pre/geom_addPrimitive.py b/processing/pre/geom_addPrimitive.py index 54de558f7..08281bd5c 100755 --- a/processing/pre/geom_addPrimitive.py +++ b/processing/pre/geom_addPrimitive.py @@ -43,7 +43,7 @@ parser.add_option('-f', '--fill', dest='fill', type='int', metavar = 'int' help='grain index to fill primitive. "0" selects maximum microstructure index + 1 [%default]') parser.add_option('-q', '--quaternion', dest='quaternion', type='float', nargs = 4, metavar=' '.join(['float']*4), help = 'rotation of primitive as quaternion') -parser.add_option('-a', '--angleaxis', dest='angleaxis', nargs = 4, metavar=' '.join(['float']*4), +parser.add_option('-a', '--angleaxis', dest='angleaxis', nargs = 4, metavar=' '.join(['float']*4), type=float, help = 'angle,x,y,z clockwise rotation of primitive about axis by angle') parser.add_option( '--degrees', dest='degrees', action='store_true', help = 'angle is given in degrees [%default]') @@ -63,14 +63,11 @@ parser.set_defaults(center = (.0,.0,.0), 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], - options.angleaxis[1:4]) + rotation = damask.Rotation.fromAngleAxis(np.array(options.angleaxis),options.degrees) elif options.quaternion is not None: - options.quaternion = list(map(float,options.quaternion)) - rotation = damask.Quaternion(quat=options.quaternion) + rotation = damask.Rotation.fromQuaternion(np.array(options.quaternion)) else: - rotation = damask.Quaternion() + rotation = damask.Rotation() options.center = np.array(options.center) options.dimension = np.array(options.dimension) diff --git a/python/damask/orientation.py b/python/damask/orientation.py index 9a92c77d4..bba13eb8a 100644 --- a/python/damask/orientation.py +++ b/python/damask/orientation.py @@ -235,7 +235,6 @@ class Rotation: resulting in new coordinates "b" when expressed in system "B". b = Q * a b = np.dot(Q.asMatrix(),a) - ToDo: Denfine how to 3x3 and 3x3x3x3 matrices """ __slots__ = ['quaternion'] @@ -351,9 +350,38 @@ class Rotation: raise ValueError('Rodriques rotation angle not positive.\n'.format(ro[0])) return cls(ro2qu(ro)) - + def __mul__(self, other): + """ + Multiplication + + Rotation: Details needed (active/passive), more cases (3,3), (3,3,3,3) need to be considered + """ + if isinstance(other, Rotation): + return self.__class__((self.quaternion * other.quaternion).asArray()) + elif isinstance(other, np.ndarray): + if other.shape == (3,): + ( 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,): + raise NotImplementedError + elif other.shape == (3,3,3,3): + raise NotImplementedError + else: + return NotImplemented + else: + return NotImplemented + # ****************************************************************************************** class Quaternion: u"""