From 9e0a0ee166b85a13e580448b1f2c0c6fad4f82c8 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 11 Jun 2022 23:04:21 +0200 Subject: [PATCH 1/4] not needed --- src/rotations.f90 | 676 +--------------------------------------------- 1 file changed, 10 insertions(+), 666 deletions(-) diff --git a/src/rotations.f90 b/src/rotations.f90 index 5367f2c07..7e0d33cff 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -61,7 +61,6 @@ module rotations procedure, public :: asQuaternion procedure, public :: asEulers procedure, public :: asAxisAngle - procedure, public :: asRodrigues procedure, public :: asMatrix !------------------------------------------ procedure, public :: fromQuaternion @@ -150,24 +149,7 @@ pure function asMatrix(self) asMatrix = qu2om(self%q) end function asMatrix -!--------------------------------------------------------------------------------------------------- -pure function asRodrigues(self) - class(tRotation), intent(in) :: self - real(pReal), dimension(4) :: asRodrigues - - asRodrigues = qu2ro(self%q) - -end function asRodrigues -!--------------------------------------------------------------------------------------------------- -pure function asHomochoric(self) - - class(tRotation), intent(in) :: self - real(pReal), dimension(3) :: asHomochoric - - asHomochoric = qu2ho(self%q) - -end function asHomochoric !--------------------------------------------------------------------------------------------------- ! Initialize rotation from different representations @@ -508,71 +490,6 @@ pure function qu2ax(qu) result(ax) end function qu2ax -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert unit quaternion to Rodrigues vector -!--------------------------------------------------------------------------------------------------- -pure function qu2ro(qu) result(ro) - - real(pReal), intent(in), dimension(4) :: qu - real(pReal), dimension(4) :: ro - - real(pReal) :: s - real(pReal), parameter :: thr = 1.0e-8_pReal - - if (abs(qu(1)) < thr) then - ro = [qu(2), qu(3), qu(4), IEEE_value(1.0_pReal,IEEE_positive_inf)] - else - s = norm2(qu(2:4)) - if (s < thr) then - ro = [0.0_pReal, 0.0_pReal, P, 0.0_pReal] - else - ro = [qu(2)/s,qu(3)/s,qu(4)/s, tan(acos(math_clip(qu(1),-1.0_pReal,1.0_pReal)))] - endif - - end if - -end function qu2ro - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert unit quaternion to homochoric -!--------------------------------------------------------------------------------------------------- -pure function qu2ho(qu) result(ho) - - real(pReal), intent(in), dimension(4) :: qu - real(pReal), dimension(3) :: ho - - real(pReal) :: omega, f - - omega = 2.0 * acos(math_clip(qu(1),-1.0_pReal,1.0_pReal)) - - if (dEq0(omega,tol=1.e-5_pReal)) then - ho = [ 0.0_pReal, 0.0_pReal, 0.0_pReal ] - else - ho = qu(2:4) - f = 0.75_pReal * ( omega - sin(omega) ) - ho = ho/norm2(ho)* f**(1.0_pReal/3.0_pReal) - end if - -end function qu2ho - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert unit quaternion to cubochoric -!--------------------------------------------------------------------------------------------------- -pure function qu2cu(qu) result(cu) - - real(pReal), intent(in), dimension(4) :: qu - real(pReal), dimension(3) :: cu - - cu = ho2cu(qu2ho(qu)) - -end function qu2cu - - !--------------------------------------------------------------------------------------------------- !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief convert rotation matrix to unit quaternion @@ -669,48 +586,6 @@ function om2ax(om) result(ax) end function om2ax -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert rotation matrix to Rodrigues vector -!--------------------------------------------------------------------------------------------------- -pure function om2ro(om) result(ro) - - real(pReal), intent(in), dimension(3,3) :: om - real(pReal), dimension(4) :: ro - - ro = eu2ro(om2eu(om)) - -end function om2ro - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert rotation matrix to homochoric -!--------------------------------------------------------------------------------------------------- -function om2ho(om) result(ho) - - real(pReal), intent(in), dimension(3,3) :: om - real(pReal), dimension(3) :: ho - - ho = ax2ho(om2ax(om)) - -end function om2ho - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert rotation matrix to cubochoric -!--------------------------------------------------------------------------------------------------- -function om2cu(om) result(cu) - - real(pReal), intent(in), dimension(3,3) :: om - real(pReal), dimension(3) :: cu - - cu = ho2cu(om2ho(om)) - -end function om2cu - - !--------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University !> @brief Euler angles to unit quaternion @@ -794,55 +669,6 @@ pure function eu2ax(eu) result(ax) end function eu2ax -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief Euler angles to Rodrigues vector -!--------------------------------------------------------------------------------------------------- -pure function eu2ro(eu) result(ro) - - real(pReal), intent(in), dimension(3) :: eu - real(pReal), dimension(4) :: ro - - ro = eu2ax(eu) - if (ro(4) >= PI) then - ro(4) = IEEE_value(ro(4),IEEE_positive_inf) - elseif(dEq0(ro(4))) then - ro = [ 0.0_pReal, 0.0_pReal, P, 0.0_pReal ] - else - ro(4) = tan(ro(4)*0.5_pReal) - end if - -end function eu2ro - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert Euler angles to homochoric -!--------------------------------------------------------------------------------------------------- -pure function eu2ho(eu) result(ho) - - real(pReal), intent(in), dimension(3) :: eu - real(pReal), dimension(3) :: ho - - ho = ax2ho(eu2ax(eu)) - -end function eu2ho - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert Euler angles to cubochoric -!--------------------------------------------------------------------------------------------------- -function eu2cu(eu) result(cu) - - real(pReal), intent(in), dimension(3) :: eu - real(pReal), dimension(3) :: cu - - cu = ho2cu(eu2ho(eu)) - -end function eu2cu - - !--------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University !> @brief convert axis angle pair to quaternion @@ -916,465 +742,6 @@ pure function ax2eu(ax) result(eu) end function ax2eu -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert axis angle pair to Rodrigues vector -!--------------------------------------------------------------------------------------------------- -pure function ax2ro(ax) result(ro) - - real(pReal), intent(in), dimension(4) :: ax - real(pReal), dimension(4) :: ro - - real(pReal), parameter :: thr = 1.0e-7_pReal - - if (dEq0(ax(4))) then - ro = [ 0.0_pReal, 0.0_pReal, P, 0.0_pReal ] - else - ro(1:3) = ax(1:3) - ! we need to deal with the 180 degree case - ro(4) = merge(IEEE_value(ro(4),IEEE_positive_inf),tan(ax(4)*0.5_pReal),abs(ax(4)-PI) < thr) - end if - -end function ax2ro - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert axis angle pair to homochoric -!--------------------------------------------------------------------------------------------------- -pure function ax2ho(ax) result(ho) - - real(pReal), intent(in), dimension(4) :: ax - real(pReal), dimension(3) :: ho - - real(pReal) :: f - - f = 0.75_pReal * ( ax(4) - sin(ax(4)) ) - f = f**(1.0_pReal/3.0_pReal) - ho = ax(1:3) * f - -end function ax2ho - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert axis angle pair to cubochoric -!--------------------------------------------------------------------------------------------------- -function ax2cu(ax) result(cu) - - real(pReal), intent(in), dimension(4) :: ax - real(pReal), dimension(3) :: cu - - cu = ho2cu(ax2ho(ax)) - -end function ax2cu - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert Rodrigues vector to unit quaternion -!--------------------------------------------------------------------------------------------------- -pure function ro2qu(ro) result(qu) - - real(pReal), intent(in), dimension(4) :: ro - real(pReal), dimension(4) :: qu - - qu = ax2qu(ro2ax(ro)) - -end function ro2qu - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert Rodrigues vector to rotation matrix -!--------------------------------------------------------------------------------------------------- -pure function ro2om(ro) result(om) - - real(pReal), intent(in), dimension(4) :: ro - real(pReal), dimension(3,3) :: om - - om = ax2om(ro2ax(ro)) - -end function ro2om - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert Rodrigues vector to Euler angles -!--------------------------------------------------------------------------------------------------- -pure function ro2eu(ro) result(eu) - - real(pReal), intent(in), dimension(4) :: ro - real(pReal), dimension(3) :: eu - - eu = om2eu(ro2om(ro)) - -end function ro2eu - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert Rodrigues vector to axis angle pair -!--------------------------------------------------------------------------------------------------- -pure function ro2ax(ro) result(ax) - - real(pReal), intent(in), dimension(4) :: ro - real(pReal), dimension(4) :: ax - - real(pReal) :: ta, angle - - ta = ro(4) - - if (.not. IEEE_is_finite(ta)) then - ax = [ ro(1), ro(2), ro(3), PI ] - elseif (dEq0(ta)) then - ax = [ 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal ] - else - angle = 2.0_pReal*atan(ta) - ta = 1.0_pReal/norm2(ro(1:3)) - ax = [ ro(1)/ta, ro(2)/ta, ro(3)/ta, angle ] - end if - -end function ro2ax - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert Rodrigues vector to homochoric -!--------------------------------------------------------------------------------------------------- -pure function ro2ho(ro) result(ho) - - real(pReal), intent(in), dimension(4) :: ro - real(pReal), dimension(3) :: ho - - real(pReal) :: f - - if (dEq0(norm2(ro(1:3)))) then - ho = [ 0.0_pReal, 0.0_pReal, 0.0_pReal ] - else - f = merge(2.0_pReal*atan(ro(4)) - sin(2.0_pReal*atan(ro(4))),PI, IEEE_is_finite(ro(4))) - ho = ro(1:3) * (0.75_pReal*f)**(1.0_pReal/3.0_pReal) - end if - -end function ro2ho - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert Rodrigues vector to cubochoric -!--------------------------------------------------------------------------------------------------- -pure function ro2cu(ro) result(cu) - - real(pReal), intent(in), dimension(4) :: ro - real(pReal), dimension(3) :: cu - - cu = ho2cu(ro2ho(ro)) - -end function ro2cu - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert homochoric to unit quaternion -!--------------------------------------------------------------------------------------------------- -pure function ho2qu(ho) result(qu) - - real(pReal), intent(in), dimension(3) :: ho - real(pReal), dimension(4) :: qu - - qu = ax2qu(ho2ax(ho)) - -end function ho2qu - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert homochoric to rotation matrix -!--------------------------------------------------------------------------------------------------- -pure function ho2om(ho) result(om) - - real(pReal), intent(in), dimension(3) :: ho - real(pReal), dimension(3,3) :: om - - om = ax2om(ho2ax(ho)) - -end function ho2om - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert homochoric to Euler angles -!--------------------------------------------------------------------------------------------------- -pure function ho2eu(ho) result(eu) - - real(pReal), intent(in), dimension(3) :: ho - real(pReal), dimension(3) :: eu - - eu = ax2eu(ho2ax(ho)) - -end function ho2eu - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert homochoric to axis angle pair -!--------------------------------------------------------------------------------------------------- -pure function ho2ax(ho) result(ax) - - real(pReal), intent(in), dimension(3) :: ho - real(pReal), dimension(4) :: ax - - integer :: i - real(pReal) :: hmag_squared, s, hm - real(pReal), parameter, dimension(16) :: & - tfit = [ 1.0000000000018852_pReal, -0.5000000002194847_pReal, & - -0.024999992127593126_pReal, -0.003928701544781374_pReal, & - -0.0008152701535450438_pReal, -0.0002009500426119712_pReal, & - -0.00002397986776071756_pReal, -0.00008202868926605841_pReal, & - +0.00012448715042090092_pReal, -0.0001749114214822577_pReal, & - +0.0001703481934140054_pReal, -0.00012062065004116828_pReal, & - +0.000059719705868660826_pReal, -0.00001980756723965647_pReal, & - +0.000003953714684212874_pReal, -0.00000036555001439719544_pReal ] - - ! normalize h and store the magnitude - hmag_squared = sum(ho**2) - if (dEq0(hmag_squared)) then - ax = [ 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal ] - else - hm = hmag_squared - - ! convert the magnitude to the rotation angle - s = tfit(1) + tfit(2) * hmag_squared - do i=3,16 - hm = hm*hmag_squared - s = s + tfit(i) * hm - enddo - ax = [ho/sqrt(hmag_squared), 2.0_pReal*acos(s)] - end if - -end function ho2ax - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert homochoric to Rodrigues vector -!--------------------------------------------------------------------------------------------------- -pure function ho2ro(ho) result(ro) - - real(pReal), intent(in), dimension(3) :: ho - real(pReal), dimension(4) :: ro - - ro = ax2ro(ho2ax(ho)) - -end function ho2ro - - -!-------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH -!> @brief convert homochoric to cubochoric -!-------------------------------------------------------------------------- -pure function ho2cu(ho) result(cu) - - real(pReal), intent(in), dimension(3) :: ho - real(pReal), dimension(3) :: cu, xyz1, xyz3 - real(pReal), dimension(2) :: Tinv, xyz2 - real(pReal) :: rs, qxy, q2, sq2, q, tt - integer, dimension(3,2) :: p - - rs = norm2(ho) - if (rs > R1+1.e-6_pReal) then - cu = IEEE_value(cu,IEEE_positive_inf) - return - endif - - center: if (all(dEq0(ho))) then - cu = 0.0_pReal - else center - p = GetPyramidOrder(ho) - xyz3 = ho(p(:,1)) - - ! inverse M_3 - xyz2 = xyz3(1:2) * sqrt( 2.0*rs/(rs+abs(xyz3(3))) ) - - ! inverse M_2 - qxy = sum(xyz2**2) - - special: if (dEq0(qxy)) then - Tinv = 0.0_pReal - else special - q2 = qxy + maxval(abs(xyz2))**2 - sq2 = sqrt(q2) - q = (BETA/R2/R1) * sqrt(q2*qxy/(q2-maxval(abs(xyz2))*sq2)) - tt = (minval(abs(xyz2))**2+maxval(abs(xyz2))*sq2)/R2/qxy - Tinv = q * sign(1.0_pReal,xyz2) * merge([ 1.0_pReal, acos(math_clip(tt,-1.0_pReal,1.0_pReal))/PI12], & - [ acos(math_clip(tt,-1.0_pReal,1.0_pReal))/PI12, 1.0_pReal], & - abs(xyz2(2)) <= abs(xyz2(1))) - endif special - - ! inverse M_1 - xyz1 = [ Tinv(1), Tinv(2), sign(1.0_pReal,xyz3(3)) * rs / PREF ]/SC - - ! reverse the coordinates back to order according to the original pyramid number - cu = xyz1(p(:,2)) - - endif center - -end function ho2cu - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert cubochoric to unit quaternion -!--------------------------------------------------------------------------------------------------- -pure function cu2qu(cu) result(qu) - - real(pReal), intent(in), dimension(3) :: cu - real(pReal), dimension(4) :: qu - - qu = ho2qu(cu2ho(cu)) - -end function cu2qu - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert cubochoric to rotation matrix -!--------------------------------------------------------------------------------------------------- -pure function cu2om(cu) result(om) - - real(pReal), intent(in), dimension(3) :: cu - real(pReal), dimension(3,3) :: om - - om = ho2om(cu2ho(cu)) - -end function cu2om - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert cubochoric to Euler angles -!--------------------------------------------------------------------------------------------------- -pure function cu2eu(cu) result(eu) - - real(pReal), intent(in), dimension(3) :: cu - real(pReal), dimension(3) :: eu - - eu = ho2eu(cu2ho(cu)) - -end function cu2eu - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert cubochoric to axis angle pair -!--------------------------------------------------------------------------------------------------- -function cu2ax(cu) result(ax) - - real(pReal), intent(in), dimension(3) :: cu - real(pReal), dimension(4) :: ax - - ax = ho2ax(cu2ho(cu)) - -end function cu2ax - - -!--------------------------------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @brief convert cubochoric to Rodrigues vector -!--------------------------------------------------------------------------------------------------- -pure function cu2ro(cu) result(ro) - - real(pReal), intent(in), dimension(3) :: cu - real(pReal), dimension(4) :: ro - - ro = ho2ro(cu2ho(cu)) - -end function cu2ro - - -!-------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH -!> @brief map from 3D cubic grid to 3D ball -!-------------------------------------------------------------------------- -pure function cu2ho(cu) result(ho) - - real(pReal), intent(in), dimension(3) :: cu - real(pReal), dimension(3) :: ho, LamXYZ, XYZ - real(pReal), dimension(2) :: T - real(pReal) :: c, s, q - real(pReal), parameter :: eps = 1.0e-8_pReal - integer, dimension(3,2) :: p - integer, dimension(2) :: order - - if (maxval(abs(cu)) > AP/2.0+eps) then - ho = IEEE_value(cu,IEEE_positive_inf) - return - end if - - ! transform to the sphere grid via the curved square, and intercept the zero point - center: if (all(dEq0(cu))) then - ho = 0.0_pReal - else center - ! get pyramide and scale by grid parameter ratio - p = GetPyramidOrder(cu) - XYZ = cu(p(:,1)) * SC - - ! intercept all the points along the z-axis - special: if (all(dEq0(XYZ(1:2)))) then - LamXYZ = [ 0.0_pReal, 0.0_pReal, PREF * XYZ(3) ] - else special - order = merge( [2,1], [1,2], abs(XYZ(2)) <= abs(XYZ(1))) ! order of absolute values of XYZ - q = PI12 * XYZ(order(1))/XYZ(order(2)) ! smaller by larger - c = cos(q) - s = sin(q) - q = PREK * XYZ(order(2))/ sqrt(R2-c) - T = [ (R2*c - 1.0), R2 * s] * q - - ! transform to sphere grid (inverse Lambert) - ! [note that there is no need to worry about dividing by zero, since XYZ(3) can not become zero] - c = sum(T**2) - s = c * PI/(24.0*XYZ(3)**2) - c = c * sqrt(PI/24.0_pReal) / XYZ(3) - q = sqrt( 1.0 - s ) - LamXYZ = [ T(order(2)) * q, T(order(1)) * q, PREF * XYZ(3) - c ] - end if special - - ! reverse the coordinates back to order according to the original pyramid number - ho = LamXYZ(p(:,2)) - - end if center - -end function cu2ho - - -!-------------------------------------------------------------------------- -!> @author Marc De Graef, Carnegie Mellon University -!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH -!> @brief Determine to which pyramid a point in a cubic grid belongs. -!-------------------------------------------------------------------------- -pure function GetPyramidOrder(xyz) - - real(pReal),intent(in),dimension(3) :: xyz - integer, dimension(3,2) :: GetPyramidOrder - - if (((abs(xyz(1)) <= xyz(3)).and.(abs(xyz(2)) <= xyz(3))) .or. & - ((abs(xyz(1)) <= -xyz(3)).and.(abs(xyz(2)) <= -xyz(3)))) then - GetPyramidOrder = reshape([[1,2,3],[1,2,3]],[3,2]) - else if (((abs(xyz(3)) <= xyz(1)).and.(abs(xyz(2)) <= xyz(1))) .or. & - ((abs(xyz(3)) <= -xyz(1)).and.(abs(xyz(2)) <= -xyz(1)))) then - GetPyramidOrder = reshape([[2,3,1],[3,1,2]],[3,2]) - else if (((abs(xyz(1)) <= xyz(2)).and.(abs(xyz(3)) <= xyz(2))) .or. & - ((abs(xyz(1)) <= -xyz(2)).and.(abs(xyz(3)) <= -xyz(2)))) then - GetPyramidOrder = reshape([[3,1,2],[2,3,1]],[3,2]) - else - GetPyramidOrder = -1 ! should be impossible, but might simplify debugging - end if - -end function GetPyramidOrder - - !-------------------------------------------------------------------------------------------------- !> @brief Multiply two quaternions. !-------------------------------------------------------------------------------------------------- @@ -1425,17 +792,15 @@ subroutine selfTest() do i = 1, 20 if(i==1) then - qu = om2qu(math_I3) + qu = [1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal] elseif(i==2) then - qu = eu2qu([0.0_pReal,0.0_pReal,0.0_pReal]) + qu = [1.0_pReal,-0.0_pReal,-0.0_pReal,-0.0_pReal] elseif(i==3) then - qu = eu2qu([TAU,PI,TAU]) + qu = [0.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal] elseif(i==4) then qu = [0.0_pReal,0.0_pReal,1.0_pReal,0.0_pReal] elseif(i==5) then - qu = ro2qu([1.0_pReal,0.0_pReal,0.0_pReal,IEEE_value(1.0_pReal, IEEE_positive_inf)]) - elseif(i==6) then - qu = ax2qu([1.0_pReal,0.0_pReal,0.0_pReal,0.0_pReal]) + qu = [0.0_pReal, 0.0_pReal, 0.0_pReal, 1.0_pReal] else call random_number(x) A = sqrt(x(3)) @@ -1448,37 +813,16 @@ subroutine selfTest() endif - if(.not. quaternion_equal(om2qu(qu2om(qu)),qu)) error stop 'om2qu/qu2om' - if(.not. quaternion_equal(eu2qu(qu2eu(qu)),qu)) error stop 'eu2qu/qu2eu' - if(.not. quaternion_equal(ax2qu(qu2ax(qu)),qu)) error stop 'ax2qu/qu2ax' - if(.not. quaternion_equal(ro2qu(qu2ro(qu)),qu)) error stop 'ro2qu/qu2ro' - if(.not. quaternion_equal(ho2qu(qu2ho(qu)),qu)) error stop 'ho2qu/qu2ho' - if(.not. quaternion_equal(cu2qu(qu2cu(qu)),qu)) error stop 'cu2qu/qu2cu' + if(.not. quaternion_equal(om2qu(qu2om(qu)),qu)) error stop 'om2qu2om' + if(.not. quaternion_equal(eu2qu(qu2eu(qu)),qu)) error stop 'eu2qu2eu' + if(.not. quaternion_equal(ax2qu(qu2ax(qu)),qu)) error stop 'ax2qu2ax' om = qu2om(qu) - if(.not. quaternion_equal(om2qu(eu2om(om2eu(om))),qu)) error stop 'eu2om/om2eu' - if(.not. quaternion_equal(om2qu(ax2om(om2ax(om))),qu)) error stop 'ax2om/om2ax' - if(.not. quaternion_equal(om2qu(ro2om(om2ro(om))),qu)) error stop 'ro2om/om2ro' - if(.not. quaternion_equal(om2qu(ho2om(om2ho(om))),qu)) error stop 'ho2om/om2ho' - if(.not. quaternion_equal(om2qu(cu2om(om2cu(om))),qu)) error stop 'cu2om/om2cu' + if(.not. quaternion_equal(om2qu(eu2om(om2eu(om))),qu)) error stop 'eu2om2eu' + if(.not. quaternion_equal(om2qu(ax2om(om2ax(om))),qu)) error stop 'ax2om2ax' eu = qu2eu(qu) - if(.not. quaternion_equal(eu2qu(ax2eu(eu2ax(eu))),qu)) error stop 'ax2eu/eu2ax' - if(.not. quaternion_equal(eu2qu(ro2eu(eu2ro(eu))),qu)) error stop 'ro2eu/eu2ro' - if(.not. quaternion_equal(eu2qu(ho2eu(eu2ho(eu))),qu)) error stop 'ho2eu/eu2ho' - if(.not. quaternion_equal(eu2qu(cu2eu(eu2cu(eu))),qu)) error stop 'cu2eu/eu2cu' - - ax = qu2ax(qu) - if(.not. quaternion_equal(ax2qu(ro2ax(ax2ro(ax))),qu)) error stop 'ro2ax/ax2ro' - if(.not. quaternion_equal(ax2qu(ho2ax(ax2ho(ax))),qu)) error stop 'ho2ax/ax2ho' - if(.not. quaternion_equal(ax2qu(cu2ax(ax2cu(ax))),qu)) error stop 'cu2ax/ax2cu' - - ro = qu2ro(qu) - if(.not. quaternion_equal(ro2qu(ho2ro(ro2ho(ro))),qu)) error stop 'ho2ro/ro2ho' - if(.not. quaternion_equal(ro2qu(cu2ro(ro2cu(ro))),qu)) error stop 'cu2ro/ro2cu' - - ho = qu2ho(qu) - if(.not. quaternion_equal(ho2qu(cu2ho(ho2cu(ho))),qu)) error stop 'cu2ho/ho2cu' + if(.not. quaternion_equal(eu2qu(ax2eu(eu2ax(eu))),qu)) error stop 'ax2eu2ax' call R%fromMatrix(om) From 9f4e354b12eb1353c46b0a3885472f8af7644f54 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 11 Jun 2022 23:12:30 +0200 Subject: [PATCH 2/4] polishing --- src/rotations.f90 | 194 +++++++++++++++++++++++++--------------------- 1 file changed, 105 insertions(+), 89 deletions(-) diff --git a/src/rotations.f90 b/src/rotations.f90 index 7e0d33cff..d548404c9 100644 --- a/src/rotations.f90 +++ b/src/rotations.f90 @@ -27,7 +27,7 @@ ! USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ! ################################################################### -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief rotation storage and conversion @@ -44,7 +44,7 @@ ! Convention 5: the rotation angle ω is limited to the interval [0, π] ! Convention 6: the real part of a quaternion is positive, Re(q) > 0 ! Convention 7: P = -1 -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- module rotations use IO @@ -111,60 +111,64 @@ subroutine rotations_init end subroutine rotations_init -!--------------------------------------------------------------------------------------------------- -! Return rotation in different representations -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- +! Return rotation in different representations. +!-------------------------------------------------------------------------------------------------- pure function asQuaternion(self) class(tRotation), intent(in) :: self real(pReal), dimension(4) :: asQuaternion + asQuaternion = self%q end function asQuaternion -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- pure function asEulers(self) class(tRotation), intent(in) :: self real(pReal), dimension(3) :: asEulers + asEulers = qu2eu(self%q) end function asEulers -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- pure function asAxisAngle(self) class(tRotation), intent(in) :: self real(pReal), dimension(4) :: asAxisAngle + asAxisAngle = qu2ax(self%q) end function asAxisAngle -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- pure function asMatrix(self) class(tRotation), intent(in) :: self real(pReal), dimension(3,3) :: asMatrix + asMatrix = qu2om(self%q) end function asMatrix - -!--------------------------------------------------------------------------------------------------- -! Initialize rotation from different representations -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- +! Initialize rotation from different representations. +!-------------------------------------------------------------------------------------------------- subroutine fromQuaternion(self,qu) class(tRotation), intent(out) :: self real(pReal), dimension(4), intent(in) :: qu + if (dNeq(norm2(qu),1.0_pReal,1.0e-8_pReal)) call IO_error(402,ext_msg='fromQuaternion') self%q = qu end subroutine fromQuaternion -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- subroutine fromEulers(self,eu,degrees) class(tRotation), intent(out) :: self @@ -173,6 +177,7 @@ subroutine fromEulers(self,eu,degrees) real(pReal), dimension(3) :: Eulers + if (.not. present(degrees)) then Eulers = eu else @@ -185,7 +190,7 @@ subroutine fromEulers(self,eu,degrees) self%q = eu2qu(Eulers) end subroutine fromEulers -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- subroutine fromAxisAngle(self,ax,degrees,P) class(tRotation), intent(out) :: self @@ -196,6 +201,7 @@ subroutine fromAxisAngle(self,ax,degrees,P) real(pReal) :: angle real(pReal),dimension(3) :: axis + if (.not. present(degrees)) then angle = ax(4) else @@ -215,51 +221,54 @@ subroutine fromAxisAngle(self,ax,degrees,P) self%q = ax2qu([axis,angle]) end subroutine fromAxisAngle -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- subroutine fromMatrix(self,om) class(tRotation), intent(out) :: self real(pReal), dimension(3,3), intent(in) :: om + if (dNeq(math_det33(om),1.0_pReal,tol=1.0e-5_pReal)) & call IO_error(402,ext_msg='fromMatrix') self%q = om2qu(om) end subroutine fromMatrix -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- -!--------------------------------------------------------------------------------------------------- -!> @brief: Rotate a rotation -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- +!> @brief: Compose rotations. +!-------------------------------------------------------------------------------------------------- pure elemental function rotRot__(self,R) result(rRot) type(tRotation) :: rRot class(tRotation), intent(in) :: self,R - rRot = tRotation(multiply_quaternion(self%q,R%q)) + + rRot = tRotation(multiplyQuaternion(self%q,R%q)) call rRot%standardize() end function rotRot__ -!--------------------------------------------------------------------------------------------------- -!> @brief quaternion representation with positive q -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- +!> @brief Convert to quaternion representation with positive q(1). +!-------------------------------------------------------------------------------------------------- pure elemental subroutine standardize(self) class(tRotation), intent(inout) :: self + if (sign(1.0_pReal,self%q(1)) < 0.0_pReal) self%q = - self%q end subroutine standardize -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University !> @brief Rotate a vector passively (default) or actively. -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- pure function rotVector(self,v,active) result(vRot) real(pReal), dimension(3) :: vRot @@ -270,6 +279,7 @@ pure function rotVector(self,v,active) result(vRot) real(pReal), dimension(4) :: v_normed, q logical :: passive + if (present(active)) then passive = .not. active else @@ -280,22 +290,20 @@ pure function rotVector(self,v,active) result(vRot) vRot = v else v_normed = [0.0_pReal,v]/norm2(v) - if (passive) then - q = multiply_quaternion(self%q, multiply_quaternion(v_normed, conjugate_quaternion(self%q))) - else - q = multiply_quaternion(conjugate_quaternion(self%q), multiply_quaternion(v_normed, self%q)) - endif + q = merge(multiplyQuaternion(self%q, multiplyQuaternion(v_normed, conjugateQuaternion(self%q))), & + multiplyQuaternion(conjugateQuaternion(self%q), multiplyQuaternion(v_normed, self%q)), & + passive) vRot = q(2:4)*norm2(v) endif end function rotVector -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University !> @brief Rotate a rank-2 tensor passively (default) or actively. !> @details: Rotation is based on rotation matrix -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- pure function rotTensor2(self,T,active) result(tRot) real(pReal), dimension(3,3) :: tRot @@ -319,11 +327,11 @@ pure function rotTensor2(self,T,active) result(tRot) end function rotTensor2 -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- !> @brief Rotate a rank-4 tensor passively (default) or actively. !> @details: rotation is based on rotation matrix !! ToDo: Need to check active/passive !!! -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- pure function rotTensor4(self,T,active) result(tRot) real(pReal), dimension(3,3,3,3) :: tRot @@ -351,11 +359,11 @@ pure function rotTensor4(self,T,active) result(tRot) end function rotTensor4 -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- !> @brief Rotate a rank-4 stiffness tensor in Voigt 6x6 notation passively (default) or actively. !> @details: https://scicomp.stackexchange.com/questions/35600 !! ToDo: Need to check active/passive !!! -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- pure function rotStiffness(self,C,active) result(cRot) real(pReal), dimension(6,6) :: cRot @@ -391,24 +399,24 @@ pure function rotStiffness(self,C,active) result(cRot) end function rotStiffness -!--------------------------------------------------------------------------------------------------- -!> @brief Misorientation. -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- +!> @brief Calculate misorientation. +!-------------------------------------------------------------------------------------------------- pure elemental function misorientation(self,other) type(tRotation) :: misorientation class(tRotation), intent(in) :: self, other - misorientation%q = multiply_quaternion(other%q, conjugate_quaternion(self%q)) + misorientation%q = multiplyQuaternion(other%q, conjugateQuaternion(self%q)) end function misorientation -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University !> @brief Convert unit quaternion to rotation matrix. -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- pure function qu2om(qu) result(om) real(pReal), intent(in), dimension(4) :: qu @@ -436,10 +444,10 @@ pure function qu2om(qu) result(om) end function qu2om -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University -!> @brief Convert unit quaternion to Euler angles. -!--------------------------------------------------------------------------------------------------- +!> @brief Convert unit quaternion to Bunge Euler angles. +!-------------------------------------------------------------------------------------------------- pure function qu2eu(qu) result(eu) real(pReal), intent(in), dimension(4) :: qu @@ -466,10 +474,10 @@ pure function qu2eu(qu) result(eu) end function qu2eu -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University -!> @brief convert unit quaternion to axis angle pair -!--------------------------------------------------------------------------------------------------- +!> @brief Convert unit quaternion to axis-angle pair. +!-------------------------------------------------------------------------------------------------- pure function qu2ax(qu) result(ax) real(pReal), intent(in), dimension(4) :: qu @@ -477,6 +485,7 @@ pure function qu2ax(qu) result(ax) real(pReal) :: omega, s + if (dEq0(sum(qu(2:4)**2))) then ax = [ 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal ] ! axis = [001] elseif (dNeq0(qu(1))) then @@ -490,11 +499,11 @@ pure function qu2ax(qu) result(ax) end function qu2ax -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH -!> @brief convert rotation matrix to unit quaternion +!> @brief Convert rotation matrix to unit quaternion. !> @details the original formulation (direct conversion) had (numerical?) issues -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- pure function om2qu(om) result(qu) real(pReal), intent(in), dimension(3,3) :: om @@ -503,6 +512,7 @@ pure function om2qu(om) result(qu) real(pReal) :: trace,s trace = math_trace33(om) + if(trace > 0.0_pReal) then s = 0.5_pReal / sqrt(trace+1.0_pReal) qu = [0.25_pReal/s, (om(3,2)-om(2,3))*s,(om(1,3)-om(3,1))*s,(om(2,1)-om(1,2))*s] @@ -525,17 +535,18 @@ pure function om2qu(om) result(qu) end function om2qu -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University -!> @brief convert orientation matrix to Euler angles +!> @brief Convert orientation matrix to Bunge Euler angles. !> @details Two step check for special cases to avoid invalid operations (not needed for python) -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- pure function om2eu(om) result(eu) real(pReal), intent(in), dimension(3,3) :: om real(pReal), dimension(3) :: eu real(pReal) :: zeta + if (dNeq(abs(om(3,3)),1.0_pReal,1.e-8_pReal)) then zeta = 1.0_pReal/sqrt(math_clip(1.0_pReal-om(3,3)**2,1e-64_pReal,1.0_pReal)) eu = [atan2(om(3,1)*zeta,-om(3,2)*zeta), & @@ -550,10 +561,10 @@ pure function om2eu(om) result(eu) end function om2eu -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University -!> @brief convert orientation matrix to axis angle pair -!--------------------------------------------------------------------------------------------------- +!> @brief Convert orientation matrix to axis-angle pair. +!-------------------------------------------------------------------------------------------------- function om2ax(om) result(ax) real(pReal), intent(in), dimension(3,3) :: om @@ -565,6 +576,7 @@ function om2ax(om) result(ax) real(pReal), dimension(3,3) :: VR, devNull, om_ integer :: ierr, i + om_ = om ! first get the rotation angle @@ -586,10 +598,10 @@ function om2ax(om) result(ax) end function om2ax -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University -!> @brief Euler angles to unit quaternion -!--------------------------------------------------------------------------------------------------- +!> @brief Convert Bunge Euler angles to unit quaternion. +!-------------------------------------------------------------------------------------------------- pure function eu2qu(eu) result(qu) real(pReal), intent(in), dimension(3) :: eu @@ -597,6 +609,7 @@ pure function eu2qu(eu) result(qu) real(pReal), dimension(3) :: ee real(pReal) :: cPhi, sPhi + ee = 0.5_pReal*eu cPhi = cos(ee(2)) @@ -611,10 +624,10 @@ pure function eu2qu(eu) result(qu) end function eu2qu -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University -!> @brief Euler angles to orientation matrix -!--------------------------------------------------------------------------------------------------- +!> @brief Convert Euler angles to orientation matrix. +!-------------------------------------------------------------------------------------------------- pure function eu2om(eu) result(om) real(pReal), intent(in), dimension(3) :: eu @@ -622,6 +635,7 @@ pure function eu2om(eu) result(om) real(pReal), dimension(3) :: c, s + c = cos(eu) s = sin(eu) @@ -640,10 +654,10 @@ pure function eu2om(eu) result(om) end function eu2om -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University -!> @brief convert euler to axis angle -!--------------------------------------------------------------------------------------------------- +!> @brief Convert Bunge Euler angles to axis-angle pair. +!-------------------------------------------------------------------------------------------------- pure function eu2ax(eu) result(ax) real(pReal), intent(in), dimension(3) :: eu @@ -651,6 +665,7 @@ pure function eu2ax(eu) result(ax) real(pReal) :: t, delta, tau, alpha, sigma + t = tan(eu(2)*0.5_pReal) sigma = 0.5_pReal*(eu(1)+eu(3)) delta = 0.5_pReal*(eu(1)-eu(3)) @@ -669,10 +684,10 @@ pure function eu2ax(eu) result(ax) end function eu2ax -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University -!> @brief convert axis angle pair to quaternion -!--------------------------------------------------------------------------------------------------- +!> @brief Convert axis-angle pair to unit quaternion. +!-------------------------------------------------------------------------------------------------- pure function ax2qu(ax) result(qu) real(pReal), intent(in), dimension(4) :: ax @@ -692,10 +707,10 @@ pure function ax2qu(ax) result(qu) end function ax2qu -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University -!> @brief convert axis angle pair to orientation matrix -!--------------------------------------------------------------------------------------------------- +!> @brief Convert axis-angle pair to orientation matrix. +!-------------------------------------------------------------------------------------------------- pure function ax2om(ax) result(om) real(pReal), intent(in), dimension(4) :: ax @@ -703,6 +718,7 @@ pure function ax2om(ax) result(om) real(pReal) :: q, c, s, omc + c = cos(ax(4)) s = sin(ax(4)) omc = 1.0_pReal-c @@ -728,15 +744,16 @@ pure function ax2om(ax) result(om) end function ax2om -!--------------------------------------------------------------------------------------------------- +!-------------------------------------------------------------------------------------------------- !> @author Marc De Graef, Carnegie Mellon University -!> @brief convert axis angle pair to Euler angles -!--------------------------------------------------------------------------------------------------- +!> @brief Convert axis-angle pair to Bunge Euler angles. +!-------------------------------------------------------------------------------------------------- pure function ax2eu(ax) result(eu) real(pReal), intent(in), dimension(4) :: ax real(pReal), dimension(3) :: eu + eu = om2eu(ax2om(ax)) end function ax2eu @@ -745,33 +762,32 @@ end function ax2eu !-------------------------------------------------------------------------------------------------- !> @brief Multiply two quaternions. !-------------------------------------------------------------------------------------------------- -pure function multiply_quaternion(qu1,qu2) +pure function multiplyQuaternion(qu1,qu2) real(pReal), dimension(4), intent(in) :: qu1, qu2 - real(pReal), dimension(4) :: multiply_quaternion + real(pReal), dimension(4) :: multiplyQuaternion - multiply_quaternion(1) = qu1(1)*qu2(1) - qu1(2)*qu2(2) - qu1(3)*qu2(3) - qu1(4)*qu2(4) - multiply_quaternion(2) = qu1(1)*qu2(2) + qu1(2)*qu2(1) + P * (qu1(3)*qu2(4) - qu1(4)*qu2(3)) - multiply_quaternion(3) = qu1(1)*qu2(3) + qu1(3)*qu2(1) + P * (qu1(4)*qu2(2) - qu1(2)*qu2(4)) - multiply_quaternion(4) = qu1(1)*qu2(4) + qu1(4)*qu2(1) + P * (qu1(2)*qu2(3) - qu1(3)*qu2(2)) + multiplyQuaternion(1) = qu1(1)*qu2(1) - qu1(2)*qu2(2) - qu1(3)*qu2(3) - qu1(4)*qu2(4) + multiplyQuaternion(2) = qu1(1)*qu2(2) + qu1(2)*qu2(1) + P * (qu1(3)*qu2(4) - qu1(4)*qu2(3)) + multiplyQuaternion(3) = qu1(1)*qu2(3) + qu1(3)*qu2(1) + P * (qu1(4)*qu2(2) - qu1(2)*qu2(4)) + multiplyQuaternion(4) = qu1(1)*qu2(4) + qu1(4)*qu2(1) + P * (qu1(2)*qu2(3) - qu1(3)*qu2(2)) -end function multiply_quaternion +end function multiplyQuaternion !-------------------------------------------------------------------------------------------------- !> @brief Calculate conjugate complex of a quaternion. !-------------------------------------------------------------------------------------------------- -pure function conjugate_quaternion(qu) +pure function conjugateQuaternion(qu) real(pReal), dimension(4), intent(in) :: qu - real(pReal), dimension(4) :: conjugate_quaternion + real(pReal), dimension(4) :: conjugateQuaternion - conjugate_quaternion = [qu(1), -qu(2), -qu(3), -qu(4)] + conjugateQuaternion = [qu(1), -qu(2), -qu(3), -qu(4)] - -end function conjugate_quaternion +end function conjugateQuaternion !-------------------------------------------------------------------------------------------------- @@ -780,8 +796,8 @@ end function conjugate_quaternion subroutine selfTest() type(tRotation) :: R - real(pReal), dimension(4) :: qu, ax, ro - real(pReal), dimension(3) :: x, eu, ho, v3 + real(pReal), dimension(4) :: qu, ax + real(pReal), dimension(3) :: x, eu, v3 real(pReal), dimension(3,3) :: om, t33 real(pReal), dimension(3,3,3,3) :: t3333 real(pReal), dimension(6,6) :: C From f60c950a7c52d6f9ffe48e3f70ae6865e94ef6a0 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sat, 11 Jun 2022 23:21:34 +0200 Subject: [PATCH 3/4] numerically better normalization slightly more expensive, but guarantees a det(R) = 1.0 --- src/math.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/math.f90 b/src/math.f90 index 25c90ccf4..b92e1870b 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -1129,7 +1129,7 @@ pure function math_rotationalPart(F) result(R) - I_U(1)*I_F(1) * transpose(F) & + I_U(1) * transpose(matmul(F,F)) & - matmul(F,C) - R = R /(I_U(1)*I_U(2)-I_U(3)) + R = R*math_det33(R)**(-1.0/3.0) end function math_rotationalPart From 0d63284c01e101612f65fb5719172b65c7cc2fc9 Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Sun, 12 Jun 2022 16:16:25 +0200 Subject: [PATCH 4/4] increase test coverage --- src/math.f90 | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/math.f90 b/src/math.f90 index b92e1870b..882706ea7 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -1129,7 +1129,7 @@ pure function math_rotationalPart(F) result(R) - I_U(1)*I_F(1) * transpose(F) & + I_U(1) * transpose(matmul(F,F)) & - matmul(F,C) - R = R*math_det33(R)**(-1.0/3.0) + R = R*math_det33(R)**(-1.0_pReal/3.0_pReal) end function math_rotationalPart @@ -1422,7 +1422,9 @@ subroutine selfTest() t33_2 = math_rotationalPart(transpose(t33)) t33 = math_rotationalPart(t33) if (any(dNeq0(matmul(t33_2,t33) - math_I3,tol=1.0e-10_pReal))) & - error stop 'math_rotationalPart' + error stop 'math_rotationalPart (forward-backward)' + if (dNeq(1.0_pReal,math_det33(math_rotationalPart(t33)),tol=1.0e-10_pReal)) & + error stop 'math_rotationalPart (determinant)' call random_number(r) d = int(r*5.0_pReal) + 1