From 9d1c1fdb92eae4b18bc60b17823ee1c8ade93d6a Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 14 Jan 2020 11:33:18 +0100 Subject: [PATCH 1/2] enforce Re(q) > 0 --- python/damask/orientation.py | 3 ++- src/rotations.f90 | 22 ++++++++++++++++++---- 2 files changed, 20 insertions(+), 5 deletions(-) diff --git a/python/damask/orientation.py b/python/damask/orientation.py index a86ba9331..dc752a94a 100644 --- a/python/damask/orientation.py +++ b/python/damask/orientation.py @@ -26,7 +26,8 @@ class Rotation: Convention 4: Euler angle triplets are implemented using the Bunge convention, with the angular ranges as [0, 2π],[0, π],[0, 2π]. Convention 5: The rotation angle ω is limited to the interval [0, π]. - Convention 6: P = -1 (as default). + Convention 6: the real part of a quaternion is positive, Re(q) > 0 + Convention 7: P = -1 (as default). Usage ----- diff --git a/src/rotations.f90 b/src/rotations.f90 index 5deb02a20..caec1cc74 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -42,7 +42,8 @@ ! Convention 4: Euler angle triplets are implemented using the Bunge convention, ! with the angular ranges as [0, 2π],[0, π],[0, 2π] ! Convention 5: the rotation angle ω is limited to the interval [0, π] -! Convention 6: P = -1 +! Convention 6: the real part of a quaternion is positive, Re(q) > 0 +! Convention 7: P = -1 !--------------------------------------------------------------------------------------------------- module rotations @@ -77,6 +78,7 @@ module rotations procedure, public :: rotTensor4 procedure, public :: rotTensor4sym procedure, public :: misorientation + procedure, public :: standardize end type rotation public :: & @@ -92,7 +94,7 @@ contains subroutine rotations_init call quaternions_init - write(6,'(/,a)') ' <<<+- rotations init -+>>>' + write(6,'(/,a)') ' <<<+- rotations init -+>>>'; flush(6) call unitTest end subroutine rotations_init @@ -237,7 +239,6 @@ end subroutine fromMatrix !--------------------------------------------------------------------------------------------------- !> @brief: Rotate a rotation -!> ToDo: completly untested !--------------------------------------------------------------------------------------------------- pure elemental function rotRot__(self,R) result(rRot) @@ -245,10 +246,23 @@ pure elemental function rotRot__(self,R) result(rRot) class(rotation), intent(in) :: self,R rRot = rotation(self%q*R%q) + call rRot%standardize() end function rotRot__ +!--------------------------------------------------------------------------------------------------- +!> @brief quaternion representation with positive q +!--------------------------------------------------------------------------------------------------- +pure elemental subroutine standardize(self) + + class(rotation), intent(inout) :: self + + if (real(self%q) < 0.0_pReal) self%q = self%q%homomorphed() + +end subroutine standardize + + !--------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University !> @brief rotate a vector passively (default) or actively @@ -375,7 +389,7 @@ pure elemental function misorientation(self,other) type(rotation) :: misorientation class(rotation), intent(in) :: self, other - misorientation%q = conjg(self%q) * other%q !ToDo: this is the convention used in math + misorientation%q = other%q * conjg(self%q) end function misorientation From d722c6db4aa1a328eb24b433e6d16c5af82e7d8e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Tue, 14 Jan 2020 11:52:22 +0100 Subject: [PATCH 2/2] clear separation between OO and imperative arguments --- src/quaternions.f90 | 55 +++++++++++++++++++++------------------------ 1 file changed, 26 insertions(+), 29 deletions(-) diff --git a/src/quaternions.f90 b/src/quaternions.f90 index 0ca404880..0215aca6e 100644 --- a/src/quaternions.f90 +++ b/src/quaternions.f90 @@ -48,10 +48,7 @@ module quaternions procedure, private :: pow_scal__ generic, public :: operator(**) => pow_quat__, pow_scal__ - procedure, public :: abs__ - procedure, public :: dot_product__ - procedure, public :: exp__ - procedure, public :: log__ + procedure, public :: abs => abs__ procedure, public :: conjg => conjg__ procedure, public :: real => real__ procedure, public :: aimag => aimag__ @@ -317,17 +314,17 @@ end function pow_scal__ !--------------------------------------------------------------------------------------------------- !> take exponential !--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function exp__(self) +type(quaternion) elemental pure function exp__(a) - class(quaternion), intent(in) :: self + class(quaternion), intent(in) :: a real(pReal) :: absImag - absImag = norm2(aimag(self)) + absImag = norm2(aimag(a)) - exp__ = merge(exp(self%w) * [ cos(absImag), & - self%x/absImag * sin(absImag), & - self%y/absImag * sin(absImag), & - self%z/absImag * sin(absImag)], & + exp__ = merge(exp(a%w) * [ cos(absImag), & + a%x/absImag * sin(absImag), & + a%y/absImag * sin(absImag), & + a%z/absImag * sin(absImag)], & IEEE_value(1.0_pReal,IEEE_SIGNALING_NAN), & dNeq0(absImag)) @@ -337,17 +334,17 @@ end function exp__ !--------------------------------------------------------------------------------------------------- !> take logarithm !--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function log__(self) +type(quaternion) elemental pure function log__(a) - class(quaternion), intent(in) :: self + class(quaternion), intent(in) :: a real(pReal) :: absImag - absImag = norm2(aimag(self)) + absImag = norm2(aimag(a)) - log__ = merge([log(abs(self)), & - self%x/absImag * acos(self%w/abs(self)), & - self%y/absImag * acos(self%w/abs(self)), & - self%z/absImag * acos(self%w/abs(self))], & + log__ = merge([log(abs(a)), & + a%x/absImag * acos(a%w/abs(a)), & + a%y/absImag * acos(a%w/abs(a)), & + a%z/absImag * acos(a%w/abs(a))], & IEEE_value(1.0_pReal,IEEE_SIGNALING_NAN), & dNeq0(absImag)) @@ -357,11 +354,11 @@ end function log__ !--------------------------------------------------------------------------------------------------- !> return norm !--------------------------------------------------------------------------------------------------- -real(pReal) elemental pure function abs__(a) +real(pReal) elemental pure function abs__(self) - class(quaternion), intent(in) :: a + class(quaternion), intent(in) :: self - abs__ = norm2([a%w,a%x,a%y,a%z]) + abs__ = norm2([self%w,self%x,self%y,self%z]) end function abs__ @@ -381,11 +378,11 @@ end function dot_product__ !--------------------------------------------------------------------------------------------------- !> take conjugate complex !--------------------------------------------------------------------------------------------------- -type(quaternion) elemental pure function conjg__(a) +type(quaternion) elemental pure function conjg__(self) - class(quaternion), intent(in) :: a + class(quaternion), intent(in) :: self - conjg__ = [a%w, -a%x, -a%y, -a%z] + conjg__ = [self%w,-self%x,-self%y,-self%z] end function conjg__ @@ -464,7 +461,7 @@ subroutine unitTest call random_number(qu) qu = (qu-0.5_pReal) * 2.0_pReal q = quaternion(qu) - + q_2= qu if(any(dNeq(q%asArray(),q_2%asArray()))) call IO_error(401,ext_msg='assign_vec__') @@ -479,10 +476,10 @@ subroutine unitTest q_2 = q / 0.5_pReal if(any(dNeq(q_2%asArray(),2.0_pReal*qu))) call IO_error(401,ext_msg='div__') - + q_2 = q * 0.3_pReal if(dNeq0(abs(q)) .and. q_2 == q) call IO_error(401,ext_msg='eq__') - + q_2 = q if(q_2 /= q) call IO_error(401,ext_msg='neq__') @@ -504,12 +501,12 @@ subroutine unitTest if(q /= conjg(q_2)) call IO_error(401,ext_msg='conjg/involution') if(dNeq(q_2%real(), q%real())) call IO_error(401,ext_msg='conjg/real') if(any(dNeq(q_2%aimag(),q%aimag()*(-1.0_pReal)))) call IO_error(401,ext_msg='conjg/aimag') - + if(abs(q) > 0.0_pReal) then q_2 = q * q%inverse() if( dNeq(real(q_2), 1.0_pReal,1.0e-15_pReal)) call IO_error(401,ext_msg='inverse/real') if(any(dNeq0(aimag(q_2), 1.0e-15_pReal))) call IO_error(401,ext_msg='inverse/aimag') - + q_2 = q/abs(q) q_2 = conjg(q_2) - inverse(q_2) if(any(dNeq0(q_2%asArray(),1.0e-15_pReal))) call IO_error(401,ext_msg='inverse/conjg')