using new rotation class

cannot rotate matrices (3,3) and (3,3,3,3) at the moment
This commit is contained in:
Martin Diehl 2019-02-12 00:11:22 +01:00
parent 7ee933b79d
commit 48b0307fab
3 changed files with 45 additions and 45 deletions

View File

@ -9,31 +9,6 @@ import damask
scriptName = os.path.splitext(os.path.basename(__file__))[0] scriptName = os.path.splitext(os.path.basename(__file__))[0]
scriptID = ' '.join([scriptName,damask.version]) 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 # MAIN
# -------------------------------------------------------------------- # --------------------------------------------------------------------
@ -133,9 +108,8 @@ if np.sum(input) != 1: parser.error('needs exactly one input format.')
(options.quaternion,4,'quaternion'), (options.quaternion,4,'quaternion'),
][np.where(input)[0][0]] # select input label that was requested ][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.Rotation.fromAngleAxis(np.array(options.crystalrotation),options.degrees) # crystal frame rotation
r = damask.Quaternion.fromAngleAxis(toRadians*options.crystalrotation[0],options.crystalrotation[1:]) # crystal frame rotation R = damask.Rotation.fromAngleAxis(np.array(options.labrotation),options.degrees) # lab frame rotation
R = damask.Quaternion.fromAngleAxis(toRadians*options. labrotation[0],options. labrotation[1:]) # lab frame rotation
# --- loop over input files ------------------------------------------------------------------------ # --- loop over input files ------------------------------------------------------------------------
@ -179,23 +153,24 @@ for name in filenames:
outputAlive = True outputAlive = True
while outputAlive and table.data_read(): # read next data line of ASCII table while outputAlive and table.data_read(): # read next data line of ASCII table
if inputtype == 'eulers': 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': elif inputtype == 'rodrigues':
o = damask.Orientation(Rodrigues = np.array(list(map(float,table.data[column:column+3])))) o = damask.Rotation.fromRodrigues(np.array(list(map(float,table.data[column:column+3]))))
elif inputtype == 'matrix':
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': elif inputtype == 'frame':
M = np.array(list(map(float,table.data[column[0]:column[0]+3] + \ M = np.array(list(map(float,table.data[column[0]:column[0]+3] + \
table.data[column[1]:column[1]+3] + \ table.data[column[1]:column[1]+3] + \
table.data[column[2]:column[2]+3]))).reshape(3,3).T table.data[column[2]:column[2]+3]))).reshape(3,3).T
o = damask.Orientation(matrix = check_matrix(M/np.linalg.norm(M,axis=0))) o = damask.Rotation.fromMatrix(M/np.linalg.norm(M,axis=0))
elif inputtype == 'quaternion':
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: for output in options.output:
if output == 'quaternion': table.data_append(o.asQuaternion()) if output == 'quaternion': table.data_append(o.asQuaternion())

View File

@ -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]') 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), parser.add_option('-q', '--quaternion', dest='quaternion', type='float', nargs = 4, metavar=' '.join(['float']*4),
help = 'rotation of primitive as quaternion') 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') help = 'angle,x,y,z clockwise rotation of primitive about axis by angle')
parser.add_option( '--degrees', dest='degrees', action='store_true', parser.add_option( '--degrees', dest='degrees', action='store_true',
help = 'angle is given in degrees [%default]') help = 'angle is given in degrees [%default]')
@ -63,14 +63,11 @@ parser.set_defaults(center = (.0,.0,.0),
if options.dimension is None: if options.dimension is None:
parser.error('no dimension specified.') parser.error('no dimension specified.')
if options.angleaxis is not None: if options.angleaxis is not None:
options.angleaxis = list(map(float,options.angleaxis)) rotation = damask.Rotation.fromAngleAxis(np.array(options.angleaxis),options.degrees)
rotation = damask.Quaternion.fromAngleAxis(np.radians(options.angleaxis[0]) if options.degrees else options.angleaxis[0],
options.angleaxis[1:4])
elif options.quaternion is not None: elif options.quaternion is not None:
options.quaternion = list(map(float,options.quaternion)) rotation = damask.Rotation.fromQuaternion(np.array(options.quaternion))
rotation = damask.Quaternion(quat=options.quaternion)
else: else:
rotation = damask.Quaternion() rotation = damask.Rotation()
options.center = np.array(options.center) options.center = np.array(options.center)
options.dimension = np.array(options.dimension) options.dimension = np.array(options.dimension)

View File

@ -235,7 +235,6 @@ class Rotation:
resulting in new coordinates "b" when expressed in system "B". resulting in new coordinates "b" when expressed in system "B".
b = Q * a b = Q * a
b = np.dot(Q.asMatrix(),a) b = np.dot(Q.asMatrix(),a)
ToDo: Denfine how to 3x3 and 3x3x3x3 matrices
""" """
__slots__ = ['quaternion'] __slots__ = ['quaternion']
@ -351,9 +350,38 @@ class Rotation:
raise ValueError('Rodriques rotation angle not positive.\n'.format(ro[0])) raise ValueError('Rodriques rotation angle not positive.\n'.format(ro[0]))
return cls(ro2qu(ro)) 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: class Quaternion:
u""" u"""