diff --git a/code/constitutive_nonlocal.f90 b/code/constitutive_nonlocal.f90 index 592355d0e..d004d7d6a 100644 --- a/code/constitutive_nonlocal.f90 +++ b/code/constitutive_nonlocal.f90 @@ -1451,7 +1451,7 @@ do n = 1,FE_NipNeighbors(mesh_element(2,el)) area = mesh_ipArea(n,ip,el) * math_norm3(surfaceNormal) surfaceNormal = surfaceNormal / math_norm3(surfaceNormal) ! normalize the surface normal to unit length - transmissivity = constitutive_nonlocal_transmissivity(misorientation(4,n), misorientation(1:3,n)) + transmissivity = constitutive_nonlocal_transmissivity(misorientation(:,n)) highOrderScheme = .false. if ( neighboring_el > 0 .and. neighboring_ip > 0 ) then ! if neighbor exists... @@ -1637,24 +1637,34 @@ endsubroutine !********************************************************************* !* transmissivity of IP interface * !********************************************************************* -function constitutive_nonlocal_transmissivity(misorientationAngle, misorientationAxis) +function constitutive_nonlocal_transmissivity(misorientation) use prec, only: pReal, & pInt +use math, only: inDeg, & + math_norm3 implicit none !* input variables -real(pReal), dimension(3), intent(in) :: misorientationAxis ! misorientation axis -real(pReal), intent(in) :: misorientationAngle ! misorientation angle +real(pReal), dimension(4), intent(in) :: misorientation ! misorientation as quaternion !* output variables real(pReal) constitutive_nonlocal_transmissivity ! transmissivity of an IP interface for dislocations !* local variables +real(pReal) misorientationAngle, & + axisNorm +real(pReal), dimension(3) :: misorientationAxis -! transmissivity depends on misorientation angle +misorientationAngle = 2.0_pReal * dacos(min(1.0_pReal, max(-1.0_pReal, misorientation(1)))) * inDeg + +misorientationAxis = misorientation(2:4) +axisNorm = math_norm3(misorientationAxis) +if (axisNorm > tiny(axisNorm)) & + misorientationAxis = misorientationAxis / axisNorm + if (misorientationAngle < 3.0_pReal) then constitutive_nonlocal_transmissivity = 1.0_pReal elseif (misorientationAngle < 10.0_pReal) then diff --git a/code/crystallite.f90 b/code/crystallite.f90 index c943aae90..9ad2549d3 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -38,7 +38,7 @@ real(pReal), dimension (:,:,:,:), allocatable :: & crystallite_Tstar0_v, & ! 2nd Piola-Kirchhoff stress vector at start of FE inc crystallite_partionedTstar0_v, & ! 2nd Piola-Kirchhoff stress vector at start of homog inc crystallite_subTstar0_v, & ! 2nd Piola-Kirchhoff stress vector at start of crystallite inc - crystallite_eulerangles ! euler angles phi1 Phi phi2 + crystallite_orientation ! orientation as quaternion real(pReal), dimension (:,:,:,:,:), allocatable :: & crystallite_Fe, & ! current "elastic" def grad (end of converged time step) crystallite_Fp, & ! current plastic def grad (end of converged time step) @@ -56,7 +56,6 @@ real(pReal), dimension (:,:,:,:,:), allocatable :: & crystallite_partionedLp0,& ! plastic velocity grad at start of homog inc crystallite_subLp0,& ! plastic velocity grad at start of crystallite inc crystallite_P, & ! 1st Piola-Kirchhoff stress per grain - crystallite_R, & ! crystal orientation (rotation matrix current -> lattice conf) crystallite_misorientation ! misorientation between two neighboring ips (only calculated for single grain IPs) real(pReal), dimension (:,:,:,:,:,:,:), allocatable :: & crystallite_dPdF, & ! individual dPdF per grain @@ -145,8 +144,7 @@ subroutine crystallite_init(Temperature) allocate(crystallite_subF0(3,3,gMax,iMax,eMax)); crystallite_subF0 = 0.0_pReal allocate(crystallite_subFp0(3,3,gMax,iMax,eMax)); crystallite_subFp0 = 0.0_pReal allocate(crystallite_subLp0(3,3,gMax,iMax,eMax)); crystallite_subLp0 = 0.0_pReal - allocate(crystallite_R(3,3,gMax,iMax,eMax)); crystallite_R = 0.0_pReal - allocate(crystallite_eulerangles(3,gMax,iMax,eMax)); crystallite_eulerangles = 0.0_pReal + allocate(crystallite_orientation(4,gMax,iMax,eMax)); crystallite_orientation = 0.0_pReal allocate(crystallite_misorientation(4,nMax,gMax,iMax,eMax)); crystallite_misorientation = 0.0_pReal allocate(crystallite_subTstar0_v(6,gMax,iMax,eMax)); crystallite_subTstar0_v = 0.0_pReal allocate(crystallite_dPdF(3,3,3,3,gMax,iMax,eMax)); crystallite_dPdF = 0.0_pReal @@ -207,6 +205,8 @@ subroutine crystallite_init(Temperature) case('volume') mySize = 1 case('orientation') + mySize = 4 + case('eulerangles') mySize = 3 case('defgrad') mySize = 9 @@ -291,8 +291,7 @@ subroutine crystallite_init(Temperature) write(6,'(a35,x,7(i5,x))') 'crystallite_subTstar0_v: ', shape(crystallite_subTstar0_v) write(6,'(a35,x,7(i5,x))') 'crystallite_dPdF: ', shape(crystallite_dPdF) write(6,'(a35,x,7(i5,x))') 'crystallite_fallbackdPdF: ', shape(crystallite_fallbackdPdF) - write(6,'(a35,x,7(i5,x))') 'crystallite_R: ', shape(crystallite_R) - write(6,'(a35,x,7(i5,x))') 'crystallite_eulerangles: ', shape(crystallite_eulerangles) + write(6,'(a35,x,7(i5,x))') 'crystallite_orientation: ', shape(crystallite_orientation) write(6,'(a35,x,7(i5,x))') 'crystallite_misorientation: ', shape(crystallite_misorientation) write(6,'(a35,x,7(i5,x))') 'crystallite_dt: ', shape(crystallite_dt) write(6,'(a35,x,7(i5,x))') 'crystallite_subdt: ', shape(crystallite_subdt) @@ -1499,7 +1498,7 @@ subroutine crystallite_orientations() use prec, only: pInt, & pReal use math, only: math_pDecomposition, & - math_RtoEuler, & + math_RtoQuaternion, & math_misorientation, & inDeg use FEsolving, only: FEsolving_execElem, & @@ -1537,8 +1536,7 @@ integer(pInt) e, & ! element index neighboringPhase, & ! phase of my neighbor neighboringStructure, & ! lattice structure of my neighbor symmetryType ! type of crystal symmetry -real(pReal), dimension(3,3) :: U, R, & ! polar decomposition of Fe - netRotation ! net rotation between two orientations +real(pReal), dimension(3,3) :: U, R logical error @@ -1551,11 +1549,9 @@ logical error call math_pDecomposition(crystallite_Fe(:,:,g,i,e), U, R, error) ! polar decomposition of Fe if (error) then call IO_warning(650, e, i, g) - crystallite_R(:,:,g,i,e) = 0.0_pReal - crystallite_eulerangles(:,g,i,e) = (/400.0, 400.0, 400.0/) ! fake orientation + crystallite_orientation(:,g,i,e) = (/1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal/) ! fake orientation else - crystallite_R(:,:,g,i,e) = transpose(R) - crystallite_eulerangles(:,g,i,e) = math_RtoEuler(crystallite_R(:,:,g,i,e)) * inDeg + crystallite_orientation(:,g,i,e) = math_RtoQuaternion(R) endif enddo @@ -1599,19 +1595,17 @@ logical error symmetryType = 0_pInt end select - call math_misorientation( crystallite_misorientation(1:3,n,1,i,e), & - crystallite_misorientation(4,n,1,i,e), & - netRotation, & - crystallite_R(:,:,1,i,e), & - crystallite_R(:,:,1,neighboring_i,neighboring_e), & - symmetryType) ! calculate misorientation + call math_misorientation( crystallite_misorientation(:,n,1,i,e), & + crystallite_orientation(:,1,i,e), & + crystallite_orientation(:,1,neighboring_i,neighboring_e), & + symmetryType) ! calculate misorientation else ! for neighbor with different phase - crystallite_misorientation(4,n,1,i,e) = 400.0_pReal ! set misorientation angle to 400 + crystallite_misorientation(:,n,1,i,e) = (/-1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal/) ! set misorientation to maximum endif else ! no existing neighbor - crystallite_misorientation(4,n,1,i,e) = 0.0_pReal ! set misorientation angle to zero + crystallite_misorientation(:,n,1,i,e) = (/1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal/) ! set misorientation to zero endif enddo endif @@ -1636,6 +1630,7 @@ function crystallite_postResults(& !*** variables and functions from other modules ***! use prec, only: pInt, & pReal + use math, only: math_QuaternionToEuler use mesh, only: mesh_element use material, only: microstructure_crystallite, & crystallite_Noutput, & @@ -1676,7 +1671,10 @@ function crystallite_postResults(& crystallite_postResults(c+1) = material_volume(g,i,e) ! grain volume (not fraction but absolute, right?) c = c + 1_pInt case ('orientation') - crystallite_postResults(c+1:c+3) = crystallite_eulerangles(:,g,i,e) ! grain orientation + crystallite_postResults(c+1:c+4) = crystallite_orientation(:,g,i,e) ! grain orientation + c = c + 4_pInt + case ('eulerangles') + crystallite_postResults(c+1:c+3) = math_QuaternionToEuler(crystallite_orientation(:,g,i,e)) ! grain orientation c = c + 3_pInt case ('defgrad') forall (k=0:2,l=0:2) crystallite_postResults(c+1+k*3+l) = crystallite_partionedF(k+1,l+1,g,i,e) diff --git a/code/material.config b/code/material.config index 4e870dd0e..60c1e7fc3 100644 --- a/code/material.config +++ b/code/material.config @@ -56,6 +56,7 @@ crystallite 1 (output) phase (output) volume (output) orientation +(output) eulerangles (output) defgrad #-------------------# diff --git a/code/math.f90 b/code/math.f90 index ac75394a5..2c59784a8 100644 --- a/code/math.f90 +++ b/code/math.f90 @@ -65,143 +65,47 @@ -! Symmetry Operations for 3 different materials -! 24 for cubic, 12 for hexagonal, 8 for tetragonal (24+12+8)x3=132 -real(pReal), dimension(3,3,44), parameter :: symOperations = & +! Symmetry operations as quaternions +! 24 for cubic, 12 for hexagonal = 36 +real(pReal), dimension(4,36), parameter :: symOperations = & reshape((/& - 1.0_pReal, 0.0_pReal, 0.0_pReal, & - 0.0_pReal, 1.0_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, 1.0_pReal, & - 0.0_pReal, 0.0_pReal, 1.0_pReal, & - 1.0_pReal, 0.0_pReal, 0.0_pReal, & - 0.0_pReal, 1.0_pReal, 0.0_pReal, & - 0.0_pReal, 1.0_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, 1.0_pReal, & - 1.0_pReal, 0.0_pReal, 0.0_pReal, & - 0.0_pReal, -1.0_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, 1.0_pReal, & - -1.0_pReal, 0.0_pReal, 0.0_pReal, & - 0.0_pReal, -1.0_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, -1.0_pReal, & - 1.0_pReal, 0.0_pReal, 0.0_pReal, & - 0.0_pReal, 1.0_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, -1.0_pReal, & - -1.0_pReal, 0.0_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, -1.0_pReal, & - 1.0_pReal, 0.0_pReal, 0.0_pReal, & - 0.0_pReal, -1.0_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, -1.0_pReal, & - -1.0_pReal, 0.0_pReal, 0.0_pReal, & - 0.0_pReal, 1.0_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, 1.0_pReal, & - -1.0_pReal, 0.0_pReal, 0.0_pReal, & - 0.0_pReal, -1.0_pReal, 0.0_pReal, & - -1.0_pReal, 0.0_pReal, 0.0_pReal, & - 0.0_pReal, 1.0_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, -1.0_pReal, & - -1.0_pReal, 0.0_pReal, 0.0_pReal, & - 0.0_pReal, -1.0_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, 1.0_pReal, & - 1.0_pReal, 0.0_pReal, 0.0_pReal, & - 0.0_pReal, -1.0_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, -1.0_pReal, & - 0.0_pReal, 0.0_pReal, -1.0_pReal, & - 0.0_pReal, -1.0_pReal, 0.0_pReal, & - -1.0_pReal, 0.0_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, 1.0_pReal, & - 0.0_pReal, -1.0_pReal, 0.0_pReal, & - 1.0_pReal, 0.0_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, 1.0_pReal, & - 0.0_pReal, 1.0_pReal, 0.0_pReal, & - -1.0_pReal, 0.0_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, -1.0_pReal, & - 0.0_pReal, 1.0_pReal, 0.0_pReal, & - 1.0_pReal, 0.0_pReal, 0.0_pReal, & - -1.0_pReal, 0.0_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, -1.0_pReal, & - 0.0_pReal, -1.0_pReal, 0.0_pReal, & - 1.0_pReal, 0.0_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, -1.0_pReal, & - 0.0_pReal, 1.0_pReal, 0.0_pReal, & - 1.0_pReal, 0.0_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, 1.0_pReal, & - 0.0_pReal, -1.0_pReal, 0.0_pReal, & - -1.0_pReal, 0.0_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, 1.0_pReal, & - 0.0_pReal, 1.0_pReal, 0.0_pReal, & - 0.0_pReal, -1.0_pReal, 0.0_pReal, & - -1.0_pReal, 0.0_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, -1.0_pReal, & - 0.0_pReal, 1.0_pReal, 0.0_pReal, & - -1.0_pReal, 0.0_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, 1.0_pReal, & - 0.0_pReal, 1.0_pReal, 0.0_pReal, & - 1.0_pReal, 0.0_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, -1.0_pReal, & - 0.0_pReal, -1.0_pReal, 0.0_pReal, & - 1.0_pReal, 0.0_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, 1.0_pReal, & - 1.0_pReal, 0.0_pReal, 0.0_pReal, & - 0.0_pReal, 1.0_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, 1.0_pReal, & - -0.5_pReal, 0.866025403_pReal, 0.0_pReal, & - -0.866025403_pReal, -0.5_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, 1.0_pReal, & - -0.5_pReal, -0.866025403_pReal, 0.0_pReal, & - 0.866025403_pReal, -0.5_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, 1.0_pReal, & - 0.5_pReal, 0.866025403_pReal, 0.0_pReal, & - -0.866025403_pReal, 0.5_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, 1.0_pReal, & - -1.0_pReal, 0.0_pReal, 0.0_pReal, & - 0.0_pReal, -1.0_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, 1.0_pReal, & - 0.5_pReal, -0.866025403_pReal, 0.0_pReal, & - 0.866025403_pReal, 0.5_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, 1.0_pReal, & - -0.5_pReal, -0.866025403_pReal, 0.0_pReal, & - -0.866025403_pReal, 0.5_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, -1.0_pReal, & - 1.0_pReal, 0.0_pReal, 0.0_pReal, & - 0.0_pReal, -1.0_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, -1.0_pReal, & - -0.5_pReal, 0.866025403_pReal, 0.0_pReal, & - 0.866025403_pReal, 0.5_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, -1.0_pReal, & - 0.5_pReal, 0.866025403_pReal, 0.0_pReal, & - 0.866025403_pReal, -0.5_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, -1.0_pReal, & - -1.0_pReal, 0.0_pReal, 0.0_pReal, & - 0.0_pReal, 1.0_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, -1.0_pReal, & - 0.5_pReal, -0.866025403_pReal, 0.0_pReal, & - -0.866025403_pReal, -0.5_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, -1.0_pReal, & - 1.0_pReal, 0.0_pReal, 0.0_pReal, & - 0.0_pReal, 1.0_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, 1.0_pReal, & - -1.0_pReal, 0.0_pReal, 0.0_pReal, & - 0.0_pReal, 1.0_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, -1.0_pReal, & - 1.0_pReal, 0.0_pReal, 0.0_pReal, & - 0.0_pReal, -1.0_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, -1.0_pReal, & - -1.0_pReal, 0.0_pReal, 0.0_pReal, & - 0.0_pReal, -1.0_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, 1.0_pReal, & - 0.0_pReal, 1.0_pReal, 0.0_pReal, & - -1.0_pReal, 0.0_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, 1.0_pReal, & - 0.0_pReal, -1.0_pReal, 0.0_pReal, & - 1.0_pReal, 0.0_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, 1.0_pReal, & - 0.0_pReal, 1.0_pReal, 0.0_pReal, & - 1.0_pReal, 0.0_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, -1.0_pReal, & - 0.0_pReal, -1.0_pReal, 0.0_pReal, & - -1.0_pReal, 0.0_pReal, 0.0_pReal, & - 0.0_pReal, 0.0_pReal, -1.0_pReal & - /),(/3,3,44/)) + 1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal, & ! cubic symmetry operations + 0.0_pReal, 0.0_pReal, 0.7071067811865476_pReal, 0.7071067811865476_pReal, & ! 2-fold symmetry + 0.0_pReal, 0.7071067811865476_pReal, 0.0_pReal, 0.7071067811865476_pReal, & + 0.0_pReal, 0.7071067811865476_pReal, 0.7071067811865476_pReal, 0.0_pReal, & + 0.0_pReal, 0.0_pReal, 0.7071067811865476_pReal, -0.7071067811865476_pReal, & + 0.0_pReal, -0.7071067811865476_pReal, 0.0_pReal, 0.7071067811865476_pReal, & + 0.0_pReal, 0.7071067811865476_pReal, -0.7071067811865476_pReal, 0.0_pReal, & + 0.5_pReal, 0.5_pReal, 0.5_pReal, 0.5_pReal, & ! 3-fold symmetry + -0.5_pReal, 0.5_pReal, 0.5_pReal, 0.5_pReal, & + 0.5_pReal, -0.5_pReal, 0.5_pReal, 0.5_pReal, & + -0.5_pReal, -0.5_pReal, 0.5_pReal, 0.5_pReal, & + 0.5_pReal, 0.5_pReal, -0.5_pReal, 0.5_pReal, & + -0.5_pReal, 0.5_pReal, -0.5_pReal, 0.5_pReal, & + 0.5_pReal, 0.5_pReal, 0.5_pReal, -0.5_pReal, & + -0.5_pReal, 0.5_pReal, 0.5_pReal, -0.5_pReal, & + 0.7071067811865476_pReal, 0.7071067811865476_pReal, 0.0_pReal, 0.0_pReal, & ! 4-fold symmetry + 0.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal, & + -0.7071067811865476_pReal, 0.7071067811865476_pReal, 0.0_pReal, 0.0_pReal, & + 0.7071067811865476_pReal, 0.0_pReal, 0.7071067811865476_pReal, 0.0_pReal, & + 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal, & + -0.7071067811865476_pReal, 0.0_pReal, 0.7071067811865476_pReal, 0.0_pReal, & + 0.7071067811865476_pReal, 0.0_pReal, 0.0_pReal, 0.7071067811865476_pReal, & + 0.0_pReal, 0.0_pReal, 0.0_pReal, 1.0_pReal, & + -0.7071067811865476_pReal, 0.0_pReal, 0.0_pReal, 0.7071067811865476_pReal, & + 1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal, & ! hexagonal symmetry operations + 0.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal, & ! 2-fold symmetry + 0.0_pReal, 0.0_pReal, 1.0_pReal, 0.0_pReal, & + 0.0_pReal, 0.5773502691896259_pReal, 1.154700538379252_pReal, 0.0_pReal, & + 0.0_pReal, -0.5773502691896259_pReal, 1.154700538379252_pReal, 0.0_pReal, & + 0.0_pReal, 1.154700538379252_pReal, 0.5773502691896259_pReal, 0.0_pReal, & + 0.0_pReal, -1.154700538379252_pReal, 0.5773502691896259_pReal, 0.0_pReal, & + 0.866025403784439_pReal, 0.0_pReal, 0.0_pReal, 0.5_pReal, & ! 6-fold symmetry + -0.866025403784439_pReal, 0.0_pReal, 0.0_pReal, 0.5_pReal, & + 0.5_pReal, 0.0_pReal, 0.0_pReal, 0.866025403784439_pReal, & + -0.5_pReal, 0.0_pReal, 0.0_pReal, 0.866025403784439_pReal, & + 0.0_pReal, 0.0_pReal, 0.0_pReal, 1.0_pReal & + /),(/4,36/)) @@ -242,92 +146,45 @@ real(pReal), dimension(3,3,44), parameter :: symOperations = & !************************************************************************** -! calculates the misorientation for 2 orientations C1 and C2 +! calculates the misorientation for 2 orientations Q1 and Q2 (needs quaternions) !************************************************************************** -subroutine math_misorientation(axis, angle, rot, ori1, ori2, symmetryType) +subroutine math_misorientation(dQ, Q1, Q2, symmetryType) use prec, only: pReal, pInt use IO, only: IO_warning implicit none !*** input variables - real(pReal), dimension(3,3), intent(in) :: ori1, & ! 1st orientation - ori2 ! 2nd orientation - integer(pInt), intent(in) :: symmetryType ! Type of crystal symmetry; 1:cubic, 2:hexagonal, 3:tetragonal + real(pReal), dimension(4), intent(in) :: Q1, & ! 1st orientation + Q2 ! 2nd orientation + integer(pInt), intent(in) :: symmetryType ! Type of crystal symmetry; 1:cubic, 2:hexagonal !*** output variables - real(pReal), intent(out) :: angle ! misorientation angle in degrees - real(pReal), dimension(3), intent(out) :: axis ! rotation axis of misorientation - real(pReal), dimension(3,3), intent(out) :: rot ! net rotation of the misorientation + real(pReal), dimension(4), intent(out) :: dQ ! misorientation !*** local variables - real(pReal) this_angle ! candidate for misorientation angle - real(pReal), dimension(3) :: this_axis ! candidate for rotation axis - real(pReal), dimension(3,3) :: this_rot ! candidate for net rotation - real(pReal), dimension(3,3) :: ori1sym, ori2sym ! symmetrical counterpart of 1st and 2nd orientation respectively - real(pReal) invNorm - integer(pInt) NsymOperations, & ! number of possible symmetry operations - startindex, endindex, & - s1, s2, & - i - real(pReal), dimension(:,:,:), allocatable :: mySymOperations ! symmetry Operations for my crystal symmetry + real(pReal), dimension(4) :: this_dQ ! candidate for misorientation + integer(pInt) s + integer(pInt), dimension(2), parameter :: NsymOperations = (/24,12/) ! number of possible symmetry operations + real(pReal), dimension(:,:), allocatable :: mySymOperations ! symmetry Operations for my crystal symmetry - axis = 0.0_pReal - angle = 0.0_pReal - rot = 0.0_pReal - - ! choose my symmetry operations according to my crystal symetry - if (symmetryType == 1_pInt) then - NsymOperations = 24_pInt - startindex = 1_pInt - elseif (symmetryType == 2_pInt) then - NsymOperations = 12_pInt - startindex = 25_pInt - elseif (symmetryType == 3_pInt) then - NsymOperations = 8_pInt - startindex = 37_pInt - else + dQ = 0.0_pReal + + if (symmetryType < 1_pInt .or. symmetryType > 2_pInt) then call IO_warning(700) return endif - allocate(mySymOperations(3,3,NsymOperations)) - endindex = startindex + NsymOperations - 1_pInt - mySymOperations = symOperations(:,:,startindex:endindex) + allocate(mySymOperations(4,NsymOperations(symmetryType))) + mySymOperations = symOperations(:,sum(NsymOperations(1:symmetryType-1))+1:sum(NsymOperations(1:symmetryType))) ! choose symmetry operations according to crystal symmetry - ! Initially set the orientation angle to maximum - angle = 2.0_pReal * pi - - ! apply symmetry operation to 1st orientation - do s1 = 1,NsymOperations - ori1sym = math_mul33x33(mySymOperations(:,:,s1), ori1) - - ! calculate possible net rotation - this_rot = math_mul33x33(ori1sym, transpose(ori2)) - - ! store smallest misorientation for an axis inside standard orientation triangle - ! calculate rotation axis - invNorm = ( (this_rot(1,2) - this_rot(2,1))**2.0_pReal & - + (this_rot(2,3) - this_rot(3,2))**2.0_pReal & - + (this_rot(3,1) - this_rot(1,3))**2.0_pReal ) ** (-0.5_pReal) - this_axis(1) = (this_rot(2,3) - this_rot(3,2)) * invNorm - this_axis(2) = (this_rot(3,1) - this_rot(1,3)) * invNorm - this_axis(3) = (this_rot(1,2) - this_rot(2,1)) * invNorm - - ! calculate rotation angle - this_angle = abs(0.5_pReal * pi - asin(0.4999999_pReal * (this_rot(1,1) + this_rot(2,2) + this_rot(3,3) - 1.0_pReal))) - - if (abs(this_angle) < angle) then - angle = this_angle - rot = this_rot - axis = this_axis - endif + dQ(1) = -1.0_pReal ! start with maximum misorientation angle + do s = 1,NsymOperations(symmetryType) ! loop ver symmetry operations + this_dQ = math_qMul( math_qConj(Q1), math_qMul(mySymOperations(:,s),Q2) ) ! misorientation candidate from Q1^-1*(sym*Q2) + if (this_dQ(1) > dQ(1)) dQ = this_dQ ! store if misorientation angle is smaller (cos is larger) than previous one enddo - ! convert angle to degrees - angle = angle * inDeg - endsubroutine @@ -681,6 +538,124 @@ endsubroutine ENDFUNCTION +!************************************************************************** +! quaternion multiplication q1xq2 = q12 +!************************************************************************** + PURE FUNCTION math_qMul(A,B) + + use prec, only: pReal, pInt + implicit none + + real(pReal), dimension(4), intent(in) :: A, B + real(pReal), dimension(4) :: math_qMul + + math_qMul(1) = A(1)*B(1) - A(2)*B(2) - A(3)*B(3) - A(4)*B(4) + math_qMul(2) = A(1)*B(2) + A(2)*B(1) + A(3)*B(4) - A(4)*B(3) + math_qMul(3) = A(1)*B(3) - A(2)*B(4) + A(3)*B(1) + A(4)*B(2) + math_qMul(4) = A(1)*B(4) + A(2)*B(3) - A(3)*B(2) + A(4)*B(1) + + ENDFUNCTION + + +!************************************************************************** +! quaternion dotproduct +!************************************************************************** + PURE FUNCTION math_qDot(A,B) + + use prec, only: pReal, pInt + implicit none + + real(pReal), dimension(4), intent(in) :: A, B + real(pReal) math_qDot + + math_qDot = A(1)*B(1) + A(2)*B(2) + A(3)*B(3) + A(4)*B(4) + + ENDFUNCTION + + +!************************************************************************** +! quaternion conjugation +!************************************************************************** + PURE FUNCTION math_qConj(Q) + + use prec, only: pReal, pInt + implicit none + + real(pReal), dimension(4), intent(in) :: Q + real(pReal), dimension(4) :: math_qConj + + math_qConj(1) = Q(1) + math_qConj(2:4) = -Q(2:4) + + ENDFUNCTION + + +!************************************************************************** +! quaternion norm +!************************************************************************** + PURE FUNCTION math_qNorm(Q) + + use prec, only: pReal, pInt + implicit none + + real(pReal), dimension(4), intent(in) :: Q + real(pReal) math_qNorm + + math_qNorm = dsqrt(max(0.0_pReal, Q(1)*Q(1) + Q(2)*Q(2) + Q(3)*Q(3) + Q(4)*Q(4))) + + ENDFUNCTION + + +!************************************************************************** +! quaternion inversion +!************************************************************************** + PURE FUNCTION math_qInv(Q) + + use prec, only: pReal, pInt + implicit none + + real(pReal), dimension(4), intent(in) :: Q + real(pReal), dimension(4) :: math_qInv + real(pReal) squareNorm + + math_qInv = 0.0_pReal + + squareNorm = math_qDot(Q,Q) + if (squareNorm > tiny(squareNorm)) & + math_qInv = math_qConj(Q) / squareNorm + + ENDFUNCTION + + +!************************************************************************** +! quaternion inversion +!************************************************************************** + PURE FUNCTION math_qRot(Q,v) + + use prec, only: pReal, pInt + implicit none + + real(pReal), dimension(4), intent(in) :: Q + real(pReal), dimension(3), intent(in) :: v + real(pReal), dimension(3) :: math_qRot + real(pReal), dimension(4,4) :: T + integer(pInt) i, j + + do i = 1,4 + do j = 1,i + T(i,j) = Q(i) * Q(j) + enddo + enddo + + math_qRot(1) = -v(1)*(T(3,3)+T(4,4)) + v(2)*(T(3,2)-T(4,1)) + v(3)*(T(4,2)+T(3,1)) + math_qRot(2) = v(1)*(T(3,2)+T(4,1)) - v(2)*(T(2,2)+T(4,4)) + v(3)*(T(4,3)-T(2,1)) + math_qRot(3) = v(1)*(T(4,2)-T(3,1)) + v(2)*(T(4,3)+T(2,1)) - v(3)*(T(2,2)+T(3,3)) + + math_qRot = 2.0_pReal * math_qRot + v + + ENDFUNCTION + + !************************************************************************** ! transposition of a 3x3 matrix !************************************************************************** @@ -1046,15 +1021,15 @@ pure function math_transpose3x3(A) !******************************************************************** ! euclidic norm of a 3x1 vector !******************************************************************** - pure function math_norm3(v3) + pure function math_norm3(v) use prec, only: pReal,pInt implicit none - real(pReal), dimension(3), intent(in) :: v3 + real(pReal), dimension(3), intent(in) :: v real(pReal) math_norm3 - math_norm3 = sqrt(v3(1)*v3(1)+v3(2)*v3(2)+v3(3)*v3(3)) + math_norm3 = dsqrt(v(1)*v(1) + v(2)*v(2) + v(3)*v(3)) return endfunction @@ -1275,6 +1250,138 @@ pure function math_transpose3x3(A) ENDFUNCTION +!**************************************************************** +! rotation matrix from Euler angles +!**************************************************************** + PURE FUNCTION math_EulerToR (Euler) + + use prec, only: pReal, pInt + implicit none + + real(pReal), dimension(3), intent(in) :: Euler + real(pReal), dimension(3,3) :: math_EulerToR + real(pReal) c1, c, c2, s1, s, s2 + + C1=COS(Euler(1)) + C=COS(Euler(2)) + C2=COS(Euler(3)) + S1=SIN(Euler(1)) + S=SIN(Euler(2)) + S2=SIN(Euler(3)) + math_EulerToR(1,1)=C1*C2-S1*S2*C + math_EulerToR(1,2)=S1*C2+C1*S2*C + math_EulerToR(1,3)=S2*S + math_EulerToR(2,1)=-C1*S2-S1*C2*C + math_EulerToR(2,2)=-S1*S2+C1*C2*C + math_EulerToR(2,3)=C2*S + math_EulerToR(3,1)=S1*S + math_EulerToR(3,2)=-C1*S + math_EulerToR(3,3)=C + return + + ENDFUNCTION + + +!******************************************************************** +! quaternion (w+ix+jy+kz) from orientation matrix +!******************************************************************** + PURE FUNCTION math_RtoQuaternion(R) + + use prec, only: pReal, pInt + implicit none + + real(pReal), dimension (3,3), intent(in) :: R + real(pReal), dimension(4) :: math_RtoQuaternion, T + + T(1) = max(0.0_pReal, 1.0_pReal + R(1,1) + R(2,2) + R(3,3)) + T(2) = max(0.0_pReal, 1.0_pReal + R(1,1) - R(2,2) - R(3,3)) + T(3) = max(0.0_pReal, 1.0_pReal - R(1,1) + R(2,2) - R(3,3)) + T(4) = max(0.0_pReal, 1.0_pReal - R(1,1) - R(2,2) + R(3,3)) + + math_RtoQuaternion = 0.5_pReal * dsqrt(T) + + math_RtoQuaternion(2) = sign(math_RtoQuaternion(2), R(3,2) - R(2,3)) + math_RtoQuaternion(3) = sign(math_RtoQuaternion(3), R(1,3) - R(3,1)) + math_RtoQuaternion(4) = sign(math_RtoQuaternion(4), R(2,1) - R(1,2)) + + ENDFUNCTION + + +!******************************************************************** +! orientation matrix from quaternion (w+ix+jy+kz) +!******************************************************************** + PURE FUNCTION math_QuaternionToR(Q) + + use prec, only: pReal, pInt + implicit none + + real(pReal), dimension(4), intent(in) :: Q + real(pReal), dimension(3,3) :: math_QuaternionToR, T + real(pReal) w2 + integer(pInt) i, j + + forall (i = 1:3, j = 1:3) & + T(i,j) = Q(i+1) * Q(j+1) + + math_QuaternionToR = (2.0_pReal*Q(1)*Q(1) - 1.0_pReal) * math_I3 + 2.0_pReal * T ! symmetrical parts of R + + w2 = 2.0_pReal * Q(1) + math_QuaternionToR(2,1) = math_QuaternionToR(2,1) + w2 * Q(4) ! skew parts of R + math_QuaternionToR(1,2) = math_QuaternionToR(1,2) - w2 * Q(4) + math_QuaternionToR(3,1) = math_QuaternionToR(3,1) - w2 * Q(3) + math_QuaternionToR(1,3) = math_QuaternionToR(1,3) + w2 * Q(3) + math_QuaternionToR(3,2) = math_QuaternionToR(3,2) + w2 * Q(2) + math_QuaternionToR(2,3) = math_QuaternionToR(2,3) - w2 * Q(2) + + ENDFUNCTION + + +!******************************************************************** +! orientation matrix from quaternion (w+ix+jy+kz) +!******************************************************************** + PURE FUNCTION math_EulerToQuaternion(eulerangles) + + use prec, only: pReal, pInt + implicit none + + real(pReal), dimension(3), intent(in) :: eulerangles + real(pReal), dimension(4) :: math_EulerToQuaternion + real(pReal), dimension(3) :: angles + real(pReal) c, s + + angles = 0.5_pReal * eulerangles * inRad + + c = dcos(angles(2)) + s = dsin(angles(2)) + + math_EulerToQuaternion(1) = dcos(angles(1)+angles(3)) * c + math_EulerToQuaternion(2) = dcos(angles(1)-angles(3)) * s + math_EulerToQuaternion(3) = dsin(angles(1)-angles(3)) * s + math_EulerToQuaternion(4) = dsin(angles(1)+angles(3)) * c + + ENDFUNCTION + + +!******************************************************************** +! orientation matrix from quaternion (w+ix+jy+kz) +!******************************************************************** + PURE FUNCTION math_QuaternionToEuler(Q) + + use prec, only: pReal, pInt + implicit none + + real(pReal), dimension(4), intent(in) :: Q + real(pReal), dimension(3) :: math_QuaternionToEuler + + math_QuaternionToEuler(1) = atan2(Q(1)*Q(3)+Q(2)*Q(4), Q(1)*Q(2)-Q(3)*Q(4)) + math_QuaternionToEuler(2) = acos(1.0_pReal-2.0_pReal*(Q(2)*Q(2)+Q(3)*Q(3))) + math_QuaternionToEuler(3) = atan2(-Q(1)*Q(3)+Q(2)*Q(4), Q(1)*Q(2)+Q(3)*Q(4)) + + math_QuaternionToEuler = math_QuaternionToEuler * inDeg + + ENDFUNCTION + + !**************************************************************** ! rotation matrix from axis and angle (in radians) !**************************************************************** @@ -1314,38 +1421,6 @@ pure function math_transpose3x3(A) ENDFUNCTION -!**************************************************************** -! rotation matrix from Euler angles -!**************************************************************** - PURE FUNCTION math_EulerToR (Euler) - - use prec, only: pReal, pInt - implicit none - - real(pReal), dimension(3), intent(in) :: Euler - real(pReal), dimension(3,3) :: math_EulerToR - real(pReal) c1, c, c2, s1, s, s2 - - C1=COS(Euler(1)) - C=COS(Euler(2)) - C2=COS(Euler(3)) - S1=SIN(Euler(1)) - S=SIN(Euler(2)) - S2=SIN(Euler(3)) - math_EulerToR(1,1)=C1*C2-S1*S2*C - math_EulerToR(1,2)=S1*C2+C1*S2*C - math_EulerToR(1,3)=S2*S - math_EulerToR(2,1)=-C1*S2-S1*C2*C - math_EulerToR(2,2)=-S1*S2+C1*C2*C - math_EulerToR(2,3)=C2*S - math_EulerToR(3,1)=S1*S - math_EulerToR(3,2)=-C1*S - math_EulerToR(3,3)=C - return - - ENDFUNCTION - - !************************************************************************** ! disorientation angle between two sets of Euler angles !************************************************************************** @@ -1604,14 +1679,14 @@ math_sampleFiberOri = math_RtoEuler(math_mul33x33(pRot,math_mul33x33(fRot,oRot)) EB2(2,2)=1.0_pReal EB3(3,3)=1.0_pReal ELSE - RHO=SQRT(-3.0_pReal*P**3.0_pReal)/9.0_pReal + RHO=DSQRT(-3.0_pReal*P**3.0_pReal)/9.0_pReal arg=-Q/RHO/2.0_pReal if(arg.GT.1) arg=1 if(arg.LT.-1) arg=-1 - PHI=ACOS(arg) - Y1=2*RHO**(1.0_pReal/3.0_pReal)*COS(PHI/3.0_pReal) - Y2=2*RHO**(1.0_pReal/3.0_pReal)*COS(PHI/3.0_pReal+2.0_pReal/3.0_pReal*PI) - Y3=2*RHO**(1.0_pReal/3.0_pReal)*COS(PHI/3.0_pReal+4.0_pReal/3.0_pReal*PI) + PHI=DACOS(arg) + Y1=2*RHO**(1.0_pReal/3.0_pReal)*DCOS(PHI/3.0_pReal) + Y2=2*RHO**(1.0_pReal/3.0_pReal)*DCOS(PHI/3.0_pReal+2.0_pReal/3.0_pReal*PI) + Y3=2*RHO**(1.0_pReal/3.0_pReal)*DCOS(PHI/3.0_pReal+4.0_pReal/3.0_pReal*PI) EW1=Y1-R/3.0_pReal EW2=Y2-R/3.0_pReal EW3=Y3-R/3.0_pReal