From ef3fc0b58ad32eea2fc0c8dbffaf99070f129c1d Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 12 Feb 2019 07:42:46 +0100 Subject: [PATCH] bugfix: wrong array indexing rotation of meshgrid tuple implemented --- processing/post/rotateData.py | 11 +++---- processing/pre/geom_addPrimitive.py | 5 ++- python/damask/orientation.py | 49 ++++++++++++++++++++++------- 3 files changed, 44 insertions(+), 21 deletions(-) diff --git a/processing/post/rotateData.py b/processing/post/rotateData.py index 95102345b..3712b8c73 100755 --- a/processing/post/rotateData.py +++ b/processing/post/rotateData.py @@ -40,9 +40,7 @@ parser.set_defaults(rotation = (0.,1.,1.,1.), if options.data is None: parser.error('no data column specified.') -toRadians = math.pi/180.0 if options.degrees else 1.0 # rescale degrees to radians -q = damask.Quaternion().fromAngleAxis(toRadians*options.rotation[0],options.rotation[1:]) -R = q.asMatrix() +r = damask.Rotation.fromAngleAxis(options.rotation,degrees) # --- loop over input files ------------------------------------------------------------------------- @@ -90,12 +88,11 @@ for name in filenames: while outputAlive and table.data_read(): # read next data line of ASCII table for v in active['vector']: column = table.label_index(v) - table.data[column:column+3] = q * np.array(list(map(float,table.data[column:column+3]))) + table.data[column:column+3] = r * np.array(list(map(float,table.data[column:column+3]))) for t in active['tensor']: column = table.label_index(t) - table.data[column:column+9] = \ - np.dot(R,np.dot(np.array(list(map(float,table.data[column:column+9]))).reshape((3,3)), - R.transpose())).reshape((9)) + table.data[column:column+9] = (r * (np.array(list(map(float,table.data[column:column+9]))))).reshape((3,3)) + outputAlive = table.data_write() # output processed line # ------------------------------------------ output finalization ----------------------------------- diff --git a/processing/pre/geom_addPrimitive.py b/processing/pre/geom_addPrimitive.py index 08281bd5c..f93cdd54f 100755 --- a/processing/pre/geom_addPrimitive.py +++ b/processing/pre/geom_addPrimitive.py @@ -63,7 +63,7 @@ parser.set_defaults(center = (.0,.0,.0), if options.dimension is None: parser.error('no dimension specified.') if options.angleaxis is not None: - rotation = damask.Rotation.fromAngleAxis(np.array(options.angleaxis),options.degrees) + rotation = damask.Rotation.fromAngleAxis(np.array(options.angleaxis),options.degrees,normalise=True) elif options.quaternion is not None: rotation = damask.Rotation.fromQuaternion(np.array(options.quaternion)) else: @@ -156,8 +156,7 @@ for name in filenames: X -= options.center[0] - 0.5 Y -= options.center[1] - 0.5 Z -= options.center[2] - 0.5 - # and then by applying the quaternion - # this should be rotation.conjugate() * (X,Y,Z), but it is this way for backwards compatibility with the older version of this script + # and then by applying the rotation (X, Y, Z) = rotation * (X, Y, Z) # and finally by scaling (we don't worry about options.dimension being negative, np.abs occurs on the microstructure = np.where... line) X /= options.dimension[0] * 0.5 diff --git a/python/damask/orientation.py b/python/damask/orientation.py index cb3b61ee2..54482332b 100644 --- a/python/damask/orientation.py +++ b/python/damask/orientation.py @@ -284,7 +284,7 @@ class Rotation: P = -1): qu = quaternion - if P > 0: qu[1:3] *= -1 # convert from P=1 to P=-1 + if P > 0: qu[1:4] *= -1 # convert from P=1 to P=-1 if qu[0] < 0.0: raise ValueError('Quaternion has negative first component.\n{}'.format(qu[0])) if not np.isclose(np.linalg.norm(qu), 1.0): @@ -307,15 +307,17 @@ class Rotation: def fromAngleAxis(cls, angleAxis, degrees = False, + normalise = False, P = -1): ax = angleAxis - if P > 0: ax[1:3] *= -1 # convert from P=1 to P=-1 - if degrees: ax[0] = np.degrees(ax[0]) + if P > 0: ax[1:4] *= -1 # convert from P=1 to P=-1 + if degrees: ax[0] = np.degrees(ax[0]) + if normalise: ax[1:4] /=np.linalg.norm(ax[1:4]) if ax[0] < 0.0 or ax[0] > np.pi: raise ValueError('Axis angle rotation angle outside of [0..π].\n'.format(ax[0])) - if not np.isclose(np.linalg.norm(ax[1:3]), 1.0): - raise ValueError('Axis angle rotation axis is not of unit length.\n{} {} {}'.format(*ax[1:3])) + if not np.isclose(np.linalg.norm(ax[1:4]), 1.0): + raise ValueError('Axis angle rotation axis is not of unit length.\n{} {} {}'.format(*ax[1:4])) return cls(ax2qu(ax)) @@ -336,12 +338,14 @@ class Rotation: @classmethod def fromRodrigues(cls, rodrigues, + normalise = False, P = -1): ro = rodrigues - if P > 0: ro[1:3] *= -1 # convert from P=1 to P=-1 - if not np.isclose(np.linalg.norm(ro[1:3]), 1.0): - raise ValueError('Rodrigues rotation axis is not of unit length.\n{} {} {}'.format(*ro[1:3])) + if P > 0: ro[1:4] *= -1 # convert from P=1 to P=-1 + if normalise: ro[1:4] /=np.linalg.norm(ro[1:4]) + if not np.isclose(np.linalg.norm(ro[1:4]), 1.0): + raise ValueError('Rodrigues rotation axis is not of unit length.\n{} {} {}'.format(*ro[1:4])) if ro[0] < 0.0: raise ValueError('Rodriques rotation angle not positive.\n'.format(ro[0])) @@ -354,10 +358,10 @@ class Rotation: Rotation: Details needed (active/passive), more cases (3,3), (3,3,3,3) need to be considered """ - if isinstance(other, Rotation): + if isinstance(other, Rotation): # rotate a rotation return self.__class__((self.quaternion * other.quaternion).asArray()) elif isinstance(other, np.ndarray): - if other.shape == (3,): + 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) @@ -369,14 +373,37 @@ class Rotation: A*Vy + B*y + C*(z*Vx - x*Vz), A*Vz + B*z + C*(x*Vy - y*Vx), ]) - elif other.shape == (3,3,): + 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): + """Inverse rotation/backward rotation""" + self.quaternion.conjugate() + return self + + def inversed(self): + """In-place inverse rotation/backward rotation""" + return self.__class__(self.quaternion.conjugated().asArray()) + # ****************************************************************************************** class Quaternion: