diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 02aec8b7a..1df0f0d9f 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -940,8 +940,6 @@ end function crystallite_push33ToRef !-------------------------------------------------------------------------------------------------- function crystallite_postResults(ipc, ip, el) use math, only: & - math_qToEuler, & - math_qToEulerAxisAngle, & math_det33, & math_I3, & inDeg diff --git a/src/material.f90 b/src/material.f90 index ad1227c58..c73058b2e 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -940,7 +940,7 @@ subroutine material_populateGrains material_texture(c,i,e) = microstructure_texture(c,micro) material_EulerAngles(1:3,c,i,e) = texture_Gauss(1:3,1,material_texture(c,i,e)) material_EulerAngles(1:3,c,i,e) = math_RtoEuler( & ! translate back to Euler angles - math_mul33x33( & ! pre-multiply + matmul( & ! pre-multiply math_EulertoR(material_EulerAngles(1:3,c,i,e)), & ! face-value orientation texture_transformation(1:3,1:3,material_texture(c,i,e)) & ! and transformation matrix ) & diff --git a/src/math.f90 b/src/math.f90 index 3338b99e3..855dfe755 100644 --- a/src/math.f90 +++ b/src/math.f90 @@ -96,10 +96,6 @@ module math math_mul33xx33, & math_mul3333xx33, & math_mul3333xx3333, & - math_mul33x33, & - math_mul33x3, & - math_mul33x3_complex, & - math_mul66x6 , & math_exp33 , & math_transpose33, & math_inv33, & @@ -134,20 +130,13 @@ module math math_RtoEuler, & math_RtoQ, & math_EulerToR, & - math_EulerToQ, & - math_EulerAxisAngleToR, & math_axisAngleToR, & math_EulerAxisAngleToQ, & math_axisAngleToQ, & math_qToRodrig, & math_qToEuler, & - math_qToEulerAxisAngle, & math_qToAxisAngle, & - math_qToR, & math_EulerMisorientation, & - math_sampleRandomOri, & - math_sampleGaussOri, & - math_sampleFiberOri, & math_sampleGaussVar, & math_symmetricEulers, & math_eigenvectorBasisSym33, & @@ -226,7 +215,7 @@ subroutine math_check character(len=64) :: error_msg real(pReal), dimension(3,3) :: R,R2 - real(pReal), dimension(3) :: Eulers,v + real(pReal), dimension(3) :: Eulers real(pReal), dimension(4) :: q,q2,axisangle ! --- check rotation dictionary --- @@ -243,26 +232,6 @@ subroutine math_check call IO_error(401,ext_msg=error_msg) endif - ! +++ q -> R -> q +++ - R = math_qToR(q) - q2 = math_RtoQ(R) - if ( any(abs( q-q2) > tol_math_check) .and. & - any(abs(-q-q2) > tol_math_check) ) then - write (error_msg, '(a,e14.6)' ) & - 'quat -> R -> quat maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2))) - call IO_error(401,ext_msg=error_msg) - endif - - ! +++ q -> euler -> q +++ - Eulers = math_qToEuler(q) - q2 = math_EulerToQ(Eulers) - if ( any(abs( q-q2) > tol_math_check) .and. & - any(abs(-q-q2) > tol_math_check) ) then - write (error_msg, '(a,e14.6)' ) & - 'quat -> euler -> quat maximum deviation ',min(maxval(abs( q-q2)),maxval(abs(-q-q2))) - call IO_error(401,ext_msg=error_msg) - endif - ! +++ R -> euler -> R +++ Eulers = math_RtoEuler(R) R2 = math_EulerToR(Eulers) @@ -272,14 +241,6 @@ subroutine math_check call IO_error(401,ext_msg=error_msg) endif - ! +++ check rotation sense of q and R +++ - v = halton([2,8,5]) ! random vector - R = math_qToR(q) - if (any(abs(matmul(R,v) - math_qRot(q,v)) > tol_math_check)) then - write (error_msg, '(a)' ) 'R(q)*v has different sense than q*v' - call IO_error(401,ext_msg=error_msg) - endif - ! +++ check vector expansion +++ if (any(abs([1.0_pReal,2.0_pReal,2.0_pReal,3.0_pReal,3.0_pReal,3.0_pReal] - & math_expand([1.0_pReal,2.0_pReal,3.0_pReal],[1,2,3,0])) > tol_math_check)) then @@ -1486,35 +1447,6 @@ pure function math_EulerToR(Euler) end function math_EulerToR -!-------------------------------------------------------------------------------------------------- -!> @brief quaternion (w+ix+jy+kz) from Bunge-Euler (3-1-3) angles (in radians) -!> @details rotation matrix is meant to represent a PASSIVE rotation, composed of INTRINSIC -!> @details rotations around the axes of the details rotating reference frame. -!> @details similar to eu2qu from "D Rowenhorst et al. Consistent representations of and -!> @details conversions between 3D rotations, Model. Simul. Mater. Sci. Eng. 23-8 (2015)", but -!> @details Q is conjucated and Q is not reversed for Q(0) < 0. -!-------------------------------------------------------------------------------------------------- -pure function math_EulerToQ(eulerangles) - - implicit none - real(pReal), dimension(3), intent(in) :: eulerangles - real(pReal), dimension(4) :: math_EulerToQ - real(pReal) :: c, s, sigma, delta - - c = cos(0.5_pReal * eulerangles(2)) - s = sin(0.5_pReal * eulerangles(2)) - sigma = 0.5_pReal * (eulerangles(1)+eulerangles(3)) - delta = 0.5_pReal * (eulerangles(1)-eulerangles(3)) - - math_EulerToQ= [c * cos(sigma), & - s * cos(delta), & - s * sin(delta), & - c * sin(sigma) ] - math_EulerToQ = math_qConj(math_EulerToQ) ! convert to passive rotation - -end function math_EulerToQ - - !-------------------------------------------------------------------------------------------------- !> @brief rotation matrix from axis and angle (in radians) !> @details rotation matrix is meant to represent a ACTIVE rotation @@ -1558,25 +1490,6 @@ pure function math_axisAngleToR(axis,omega) end function math_axisAngleToR -!-------------------------------------------------------------------------------------------------- -!> @brief rotation matrix from axis and angle (in radians) -!> @details rotation matrix is meant to represent a PASSIVE rotation -!> @details (see http://en.wikipedia.org/wiki/Euler_angles for definitions) -!> @details eq-uivalent to eu2qu (P=+1) from "D Rowenhorst et al. Consistent representations of and -!> @details conversions between 3D rotations, Model. Simul. Mater. Sci. Eng. 23-8 (2015)" -!-------------------------------------------------------------------------------------------------- -pure function math_EulerAxisAngleToR(axis,omega) - - implicit none - real(pReal), dimension(3,3) :: math_EulerAxisAngleToR - real(pReal), dimension(3), intent(in) :: axis - real(pReal), intent(in) :: omega - - math_EulerAxisAngleToR = transpose(math_axisAngleToR(axis,omega)) ! convert to passive rotation - -end function math_EulerAxisAngleToR - - !-------------------------------------------------------------------------------------------------- !> @brief quaternion (w+ix+jy+kz) from Euler axis and angle (in radians) !> @details quaternion is meant to represent a PASSIVE rotation @@ -1625,30 +1538,6 @@ pure function math_axisAngleToQ(axis,omega) end function math_axisAngleToQ -!-------------------------------------------------------------------------------------------------- -!> @brief orientation matrix from quaternion (w+ix+jy+kz) -!> @details taken from http://arxiv.org/pdf/math/0701759v1.pdf -!> @details see also http://en.wikipedia.org/wiki/Rotation_formalisms_in_three_dimensions -!-------------------------------------------------------------------------------------------------- -pure function math_qToR(q) - - implicit none - real(pReal), dimension(4), intent(in) :: q - real(pReal), dimension(3,3) :: math_qToR, T,S - integer :: i, j - - forall(i = 1:3, j = 1:3) T(i,j) = q(i+1) * q(j+1) - - S = reshape( [0.0_pReal, -q(4), q(3), & - q(4), 0.0_pReal, -q(2), & - -q(3), q(2), 0.0_pReal],[3,3]) ! notation is transposed - - math_qToR = (2.0_pReal * q(1)*q(1) - 1.0_pReal) * math_I3 & - + 2.0_pReal * T - 2.0_pReal * q(1) * S - -end function math_qToR - - !-------------------------------------------------------------------------------------------------- !> @brief 3-1-3 Euler angles (in radians) from quaternion (w+ix+jy+kz) !> @details quaternion is meant to represent a PASSIVE rotation, @@ -1707,24 +1596,6 @@ pure function math_qToAxisAngle(Q) end function math_qToAxisAngle -!-------------------------------------------------------------------------------------------------- -!> @brief Euler axis-angle (x, y, z, ang in radians) from quaternion (w+ix+jy+kz) -!> @details quaternion is meant to represent a PASSIVE rotation -!> @details (see http://en.wikipedia.org/wiki/Euler_angles for definitions) -!-------------------------------------------------------------------------------------------------- -pure function math_qToEulerAxisAngle(qPassive) - - implicit none - real(pReal), dimension(4), intent(in) :: qPassive - real(pReal), dimension(4) :: q - real(pReal), dimension(4) :: math_qToEulerAxisAngle - - q = math_qConj(qPassive) ! convert to active rotation - math_qToEulerAxisAngle = math_qToAxisAngle(q) - -end function math_qToEulerAxisAngle - - !-------------------------------------------------------------------------------------------------- !> @brief Rodrigues vector (x, y, z) from unit quaternion (w+ix+jy+kz) !-------------------------------------------------------------------------------------------------- @@ -1760,118 +1631,6 @@ real(pReal) pure function math_EulerMisorientation(EulerA,EulerB) end function math_EulerMisorientation -!-------------------------------------------------------------------------------------------------- -!> @brief draw a random sample from Euler space -!-------------------------------------------------------------------------------------------------- -function math_sampleRandomOri() - - implicit none - real(pReal), dimension(3) :: math_sampleRandomOri, rnd - - rnd = halton([1,7,3]) - math_sampleRandomOri = [rnd(1)*2.0_pReal*PI, & - acos(2.0_pReal*rnd(2)-1.0_pReal), & - rnd(3)*2.0_pReal*PI] - -end function math_sampleRandomOri - - -!-------------------------------------------------------------------------------------------------- -!> @brief draw a sample from an Gaussian distribution around given orientation and Full Width -! at Half Maximum (FWHM) -!> @details: A uniform misorientation (limited to 2*FWHM) is sampled followed by convolution with -! a Gausian distribution -!-------------------------------------------------------------------------------------------------- -function math_sampleGaussOri(center,FWHM) - - implicit none - real(pReal), intent(in) :: FWHM - real(pReal), dimension(3), intent(in) :: center - real(pReal) :: angle - real(pReal), dimension(3) :: math_sampleGaussOri, axis - real(pReal), dimension(4) :: rnd - real(pReal), dimension(3,3) :: R - - if (FWHM < 0.1_pReal*INRAD) then - math_sampleGaussOri = center - else - GaussConvolution: do - rnd = halton([8,3,6,11]) - axis(1) = rnd(1)*2.0_pReal-1.0_pReal ! uniform on [-1,1] - axis(2:3) = [sqrt(1.0-axis(1)**2.0_pReal)*cos(rnd(2)*2.0*PI),& - sqrt(1.0-axis(1)**2.0_pReal)*sin(rnd(2)*2.0*PI)] ! random axis - angle = (rnd(3)-0.5_pReal)*4.0_pReal*FWHM ! rotation by [0, +-2 FWHM] - R = math_axisAngleToR(axis,angle) - angle = math_EulerMisorientation([0.0_pReal,0.0_pReal,0.0_pReal],math_RtoEuler(R)) - if (rnd(4) <= exp(-4.0_pReal*log(2.0_pReal)*(angle/FWHM)**2_pReal)) exit ! rejection sampling (Gaussian) - enddo GaussConvolution - math_sampleGaussOri = math_RtoEuler(matmul(R,math_EulerToR(center))) - endif - -end function math_sampleGaussOri - - -!-------------------------------------------------------------------------------------------------- -!> @brief draw a sample from an Gaussian distribution around given fiber texture and Full Width -! at Half Maximum (FWHM) -!------------------------------------------------------------------------------------------------- -function math_sampleFiberOri(alpha,beta,FWHM) - - implicit none - real(pReal), dimension(2), intent(in) :: alpha,beta - real(pReal), intent(in) :: FWHM - real(pReal), dimension(3) :: math_sampleFiberOri, & - fInC,& !< fiber axis in crystal coordinate system - fInS,& !< fiber axis in sample coordinate system - u - real(pReal), dimension(3) :: rnd - real(pReal), dimension(:),allocatable :: a !< 2D vector to tilt - integer, dimension(:),allocatable :: idx !< components of 2D vector - real(pReal), dimension(3,3) :: R !< Rotation matrix (composed of three components) - real(pReal):: angle,c - integer:: j,& !< index of smallest component - i - - allocate(a(0)) - allocate(idx(0)) - fInC = [sin(alpha(1))*cos(alpha(2)), sin(alpha(1))*sin(alpha(2)), cos(alpha(1))] - fInS = [sin(beta(1))*cos(beta(2)), sin(beta(1))*sin(beta(2)), cos(beta(1))] - - R = math_EulerAxisAngleToR(math_crossproduct(fInC,fInS),-acos(dot_product(fInC,fInS))) !< rotation to align fiber axis in crystal and sample system - - rnd = halton([7,10,3]) - R = matmul(R,math_EulerAxisAngleToR(fInS,rnd(1)*2.0_pReal*PI)) !< additional rotation (0..360deg) perpendicular to fiber axis - - if (FWHM > 0.1_pReal*INRAD) then - reducedTo2D: do i=1,3 - if (i /= minloc(abs(fInS),1)) then - a=[a,fInS(i)] - idx=[idx,i] - else - j = i - endif - enddo reducedTo2D - GaussConvolution: do - angle = (rnd(2)-0.5_pReal)*4.0_pReal*FWHM ! rotation by [0, +-2 FWHM] - ! solve cos(angle) = dot_product(fInS,u) under the assumption that their smallest component is the same - c = cos(angle)-fInS(j)**2 - u(idx(2)) = -(2.0_pReal*c*a(2) + sqrt(4*((c*a(2))**2-sum(a**2)*(c**2-a(1)**2*(1-fInS(j)**2)))))/& - (2*sum(a**2)) - u(idx(1)) = sqrt(1-u(idx(2))**2-fInS(j)**2) - u(j) = fInS(j) - - rejectionSampling: if (rnd(3) <= exp(-4.0_pReal*log(2.0_pReal)*(angle/FWHM)**2_pReal)) then - R = matmul(R,math_EulerAxisAngleToR(math_crossproduct(u,fInS),angle)) ! tilt around direction of smallest component - exit - endif rejectionSampling - rnd = halton([7,10,3]) - enddo GaussConvolution - endif - math_sampleFiberOri = math_RtoEuler(R) - -end function math_sampleFiberOri - - !-------------------------------------------------------------------------------------------------- !> @brief draw a random sample from Gauss variable !--------------------------------------------------------------------------------------------------