223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

changed internal representation of orientation and misorientation from euler angles to quaternions (this should also fix some bugs in the math_misorientation subroutine). includes a couple of new functions in math.f90 and some changes in crystallite.f90. 

beware that crystallite output "orientation" now by default returns the orientation as quaternion. if you want euler angles instead, you have to add "eulerangles"  as a crystallite output in your material.config file (see material.config template).

for input of orientations in the texture block of the material.config you still have to specify the rotation in terms of euler angles, quaternions are not yet supported for input.
This commit is contained in:
Christoph Kords 2010-03-18 12:23:17 +00:00
parent e6fdfdfc36
commit 948c119ee9
4 changed files with 353 additions and 269 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
code/constitutive_nonlocal.f90
View File

@ -1451,7 +1451,7 @@ do n = 1,FE_NipNeighbors(mesh_element(2,el))
area = mesh_ipArea(n,ip,el) * math_norm3(surfaceNormal) area = mesh_ipArea(n,ip,el) * math_norm3(surfaceNormal)
surfaceNormal = surfaceNormal / math_norm3(surfaceNormal) ! normalize the surface normal to unit length 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. highOrderScheme = .false.
if ( neighboring_el > 0 .and. neighboring_ip > 0 ) then ! if neighbor exists... if ( neighboring_el > 0 .and. neighboring_ip > 0 ) then ! if neighbor exists...
@ -1637,24 +1637,34 @@ endsubroutine
!********************************************************************* !*********************************************************************
!* transmissivity of IP interface * !* transmissivity of IP interface *
!********************************************************************* !*********************************************************************
function constitutive_nonlocal_transmissivity(misorientationAngle, misorientationAxis) function constitutive_nonlocal_transmissivity(misorientation)
use prec, only: pReal, & use prec, only: pReal, &
pInt pInt
use math, only: inDeg, &
math_norm3
implicit none implicit none
!* input variables !* input variables
real(pReal), dimension(3), intent(in) :: misorientationAxis ! misorientation axis real(pReal), dimension(4), intent(in) :: misorientation ! misorientation as quaternion
real(pReal), intent(in) :: misorientationAngle ! misorientation angle
!* output variables !* output variables
real(pReal) constitutive_nonlocal_transmissivity ! transmissivity of an IP interface for dislocations real(pReal) constitutive_nonlocal_transmissivity ! transmissivity of an IP interface for dislocations
!* local variables !* 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 if (misorientationAngle < 3.0_pReal) then
constitutive_nonlocal_transmissivity = 1.0_pReal constitutive_nonlocal_transmissivity = 1.0_pReal
elseif (misorientationAngle < 10.0_pReal) then elseif (misorientationAngle < 10.0_pReal) then

42
code/crystallite.f90
View File

@ -38,7 +38,7 @@ real(pReal), dimension (:,:,:,:), allocatable :: &
crystallite_Tstar0_v, & ! 2nd Piola-Kirchhoff stress vector at start of FE inc 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_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_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 :: & real(pReal), dimension (:,:,:,:,:), allocatable :: &
crystallite_Fe, & ! current "elastic" def grad (end of converged time step) crystallite_Fe, & ! current "elastic" def grad (end of converged time step)
crystallite_Fp, & ! current plastic 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_partionedLp0,& ! plastic velocity grad at start of homog inc
crystallite_subLp0,& ! plastic velocity grad at start of crystallite inc crystallite_subLp0,& ! plastic velocity grad at start of crystallite inc
crystallite_P, & ! 1st Piola-Kirchhoff stress per grain 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) crystallite_misorientation ! misorientation between two neighboring ips (only calculated for single grain IPs)
real(pReal), dimension (:,:,:,:,:,:,:), allocatable :: & real(pReal), dimension (:,:,:,:,:,:,:), allocatable :: &
crystallite_dPdF, & ! individual dPdF per grain 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_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_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_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_orientation(4,gMax,iMax,eMax)); crystallite_orientation = 0.0_pReal
allocate(crystallite_eulerangles(3,gMax,iMax,eMax)); crystallite_eulerangles = 0.0_pReal
allocate(crystallite_misorientation(4,nMax,gMax,iMax,eMax)); crystallite_misorientation = 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_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 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') case('volume')
mySize = 1 mySize = 1
case('orientation') case('orientation')
mySize = 4
case('eulerangles')
mySize = 3 mySize = 3
case('defgrad') case('defgrad')
mySize = 9 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_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_dPdF: ', shape(crystallite_dPdF)
write(6,'(a35,x,7(i5,x))') 'crystallite_fallbackdPdF: ', shape(crystallite_fallbackdPdF) 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_orientation: ', shape(crystallite_orientation)
write(6,'(a35,x,7(i5,x))') 'crystallite_eulerangles: ', shape(crystallite_eulerangles)
write(6,'(a35,x,7(i5,x))') 'crystallite_misorientation: ', shape(crystallite_misorientation) 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_dt: ', shape(crystallite_dt)
write(6,'(a35,x,7(i5,x))') 'crystallite_subdt: ', shape(crystallite_subdt) write(6,'(a35,x,7(i5,x))') 'crystallite_subdt: ', shape(crystallite_subdt)
@ -1499,7 +1498,7 @@ subroutine crystallite_orientations()
use prec, only: pInt, & use prec, only: pInt, &
pReal pReal
use math, only: math_pDecomposition, & use math, only: math_pDecomposition, &
math_RtoEuler, & math_RtoQuaternion, &
math_misorientation, & math_misorientation, &
inDeg inDeg
use FEsolving, only: FEsolving_execElem, & use FEsolving, only: FEsolving_execElem, &
@ -1537,8 +1536,7 @@ integer(pInt) e, & ! element index
neighboringPhase, & ! phase of my neighbor neighboringPhase, & ! phase of my neighbor
neighboringStructure, & ! lattice structure of my neighbor neighboringStructure, & ! lattice structure of my neighbor
symmetryType ! type of crystal symmetry symmetryType ! type of crystal symmetry
real(pReal), dimension(3,3) :: U, R, & ! polar decomposition of Fe real(pReal), dimension(3,3) :: U, R
netRotation ! net rotation between two orientations
logical error logical error
@ -1551,11 +1549,9 @@ logical error
call math_pDecomposition(crystallite_Fe(:,:,g,i,e), U, R, error) ! polar decomposition of Fe call math_pDecomposition(crystallite_Fe(:,:,g,i,e), U, R, error) ! polar decomposition of Fe
if (error) then if (error) then
call IO_warning(650, e, i, g) call IO_warning(650, e, i, g)
crystallite_R(:,:,g,i,e) = 0.0_pReal crystallite_orientation(:,g,i,e) = (/1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal/) ! fake orientation
crystallite_eulerangles(:,g,i,e) = (/400.0, 400.0, 400.0/) ! fake orientation
else else
crystallite_R(:,:,g,i,e) = transpose(R) crystallite_orientation(:,g,i,e) = math_RtoQuaternion(R)
crystallite_eulerangles(:,g,i,e) = math_RtoEuler(crystallite_R(:,:,g,i,e)) * inDeg
endif endif
enddo enddo
@ -1599,19 +1595,17 @@ logical error
symmetryType = 0_pInt symmetryType = 0_pInt
end select end select
call math_misorientation( crystallite_misorientation(1:3,n,1,i,e), & call math_misorientation( crystallite_misorientation(:,n,1,i,e), &
crystallite_misorientation(4,n,1,i,e), & crystallite_orientation(:,1,i,e), &
netRotation, & crystallite_orientation(:,1,neighboring_i,neighboring_e), &
crystallite_R(:,:,1,i,e), & symmetryType) ! calculate misorientation
crystallite_R(:,:,1,neighboring_i,neighboring_e), &
symmetryType) ! calculate misorientation
else ! for neighbor with different phase 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 endif
else ! no existing neighbor 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 endif
enddo enddo
endif endif
@ -1636,6 +1630,7 @@ function crystallite_postResults(&
!*** variables and functions from other modules ***! !*** variables and functions from other modules ***!
use prec, only: pInt, & use prec, only: pInt, &
pReal pReal
use math, only: math_QuaternionToEuler
use mesh, only: mesh_element use mesh, only: mesh_element
use material, only: microstructure_crystallite, & use material, only: microstructure_crystallite, &
crystallite_Noutput, & 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?) crystallite_postResults(c+1) = material_volume(g,i,e) ! grain volume (not fraction but absolute, right?)
c = c + 1_pInt c = c + 1_pInt
case ('orientation') 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 c = c + 3_pInt
case ('defgrad') 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) forall (k=0:2,l=0:2) crystallite_postResults(c+1+k*3+l) = crystallite_partionedF(k+1,l+1,g,i,e)

1
code/material.config
View File

@ -56,6 +56,7 @@ crystallite 1
(output) phase (output) phase
(output) volume (output) volume
(output) orientation (output) orientation
(output) eulerangles
(output) defgrad (output) defgrad
#-------------------# #-------------------#

557
code/math.f90
View File

@ -65,143 +65,47 @@
! Symmetry Operations for 3 different materials ! Symmetry operations as quaternions
! 24 for cubic, 12 for hexagonal, 8 for tetragonal (24+12+8)x3=132 ! 24 for cubic, 12 for hexagonal = 36
real(pReal), dimension(3,3,44), parameter :: symOperations = & real(pReal), dimension(4,36), parameter :: symOperations = &
reshape((/& reshape((/&
1.0_pReal, 0.0_pReal, 0.0_pReal, & 1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal, & ! cubic symmetry operations
0.0_pReal, 1.0_pReal, 0.0_pReal, & 0.0_pReal, 0.0_pReal, 0.7071067811865476_pReal, 0.7071067811865476_pReal, & ! 2-fold symmetry
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.0_pReal, 1.0_pReal, & 0.0_pReal, 0.7071067811865476_pReal, 0.7071067811865476_pReal, 0.0_pReal, &
1.0_pReal, 0.0_pReal, 0.0_pReal, & 0.0_pReal, 0.0_pReal, 0.7071067811865476_pReal, -0.7071067811865476_pReal, &
0.0_pReal, 1.0_pReal, 0.0_pReal, & 0.0_pReal, -0.7071067811865476_pReal, 0.0_pReal, 0.7071067811865476_pReal, &
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.0_pReal, 1.0_pReal, & 0.5_pReal, 0.5_pReal, 0.5_pReal, 0.5_pReal, & ! 3-fold symmetry
1.0_pReal, 0.0_pReal, 0.0_pReal, & -0.5_pReal, 0.5_pReal, 0.5_pReal, 0.5_pReal, &
0.0_pReal, -1.0_pReal, 0.0_pReal, & 0.5_pReal, -0.5_pReal, 0.5_pReal, 0.5_pReal, &
0.0_pReal, 0.0_pReal, 1.0_pReal, & -0.5_pReal, -0.5_pReal, 0.5_pReal, 0.5_pReal, &
-1.0_pReal, 0.0_pReal, 0.0_pReal, & 0.5_pReal, 0.5_pReal, -0.5_pReal, 0.5_pReal, &
0.0_pReal, -1.0_pReal, 0.0_pReal, & -0.5_pReal, 0.5_pReal, -0.5_pReal, 0.5_pReal, &
0.0_pReal, 0.0_pReal, -1.0_pReal, & 0.5_pReal, 0.5_pReal, 0.5_pReal, -0.5_pReal, &
1.0_pReal, 0.0_pReal, 0.0_pReal, & -0.5_pReal, 0.5_pReal, 0.5_pReal, -0.5_pReal, &
0.0_pReal, 1.0_pReal, 0.0_pReal, & 0.7071067811865476_pReal, 0.7071067811865476_pReal, 0.0_pReal, 0.0_pReal, & ! 4-fold symmetry
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, 0.0_pReal, & -0.7071067811865476_pReal, 0.7071067811865476_pReal, 0.0_pReal, 0.0_pReal, &
0.0_pReal, 0.0_pReal, -1.0_pReal, & 0.7071067811865476_pReal, 0.0_pReal, 0.7071067811865476_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, &
0.0_pReal, -1.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.7071067811865476_pReal, 0.0_pReal, 0.0_pReal, 0.7071067811865476_pReal, &
-1.0_pReal, 0.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.7071067811865476_pReal, 0.0_pReal, 0.0_pReal, 0.7071067811865476_pReal, &
0.0_pReal, 0.0_pReal, 1.0_pReal, & 1.0_pReal, 0.0_pReal, 0.0_pReal, 0.0_pReal, & ! hexagonal symmetry operations
-1.0_pReal, 0.0_pReal, 0.0_pReal, & 0.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal, & ! 2-fold symmetry
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, 0.0_pReal, & 0.0_pReal, 0.5773502691896259_pReal, 1.154700538379252_pReal, 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.0_pReal, -1.0_pReal, & 0.0_pReal, 1.154700538379252_pReal, 0.5773502691896259_pReal, 0.0_pReal, &
-1.0_pReal, 0.0_pReal, 0.0_pReal, & 0.0_pReal, -1.154700538379252_pReal, 0.5773502691896259_pReal, 0.0_pReal, &
0.0_pReal, -1.0_pReal, 0.0_pReal, & 0.866025403784439_pReal, 0.0_pReal, 0.0_pReal, 0.5_pReal, & ! 6-fold symmetry
0.0_pReal, 0.0_pReal, 1.0_pReal, & -0.866025403784439_pReal, 0.0_pReal, 0.0_pReal, 0.5_pReal, &
1.0_pReal, 0.0_pReal, 0.0_pReal, & 0.5_pReal, 0.0_pReal, 0.0_pReal, 0.866025403784439_pReal, &
0.0_pReal, -1.0_pReal, 0.0_pReal, & -0.5_pReal, 0.0_pReal, 0.0_pReal, 0.866025403784439_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, & /),(/4,36/))
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/))
@ -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 prec, only: pReal, pInt
use IO, only: IO_warning use IO, only: IO_warning
implicit none implicit none
!*** input variables !*** input variables
real(pReal), dimension(3,3), intent(in) :: ori1, & ! 1st orientation real(pReal), dimension(4), intent(in) :: Q1, & ! 1st orientation
ori2 ! 2nd orientation Q2 ! 2nd orientation
integer(pInt), intent(in) :: symmetryType ! Type of crystal symmetry; 1:cubic, 2:hexagonal, 3:tetragonal integer(pInt), intent(in) :: symmetryType ! Type of crystal symmetry; 1:cubic, 2:hexagonal
!*** output variables !*** output variables
real(pReal), intent(out) :: angle ! misorientation angle in degrees real(pReal), dimension(4), intent(out) :: dQ ! misorientation
real(pReal), dimension(3), intent(out) :: axis ! rotation axis of misorientation
real(pReal), dimension(3,3), intent(out) :: rot ! net rotation of the misorientation
!*** local variables !*** local variables
real(pReal) this_angle ! candidate for misorientation angle real(pReal), dimension(4) :: this_dQ ! candidate for misorientation
real(pReal), dimension(3) :: this_axis ! candidate for rotation axis integer(pInt) s
real(pReal), dimension(3,3) :: this_rot ! candidate for net rotation integer(pInt), dimension(2), parameter :: NsymOperations = (/24,12/) ! number of possible symmetry operations
real(pReal), dimension(3,3) :: ori1sym, ori2sym ! symmetrical counterpart of 1st and 2nd orientation respectively real(pReal), dimension(:,:), allocatable :: mySymOperations ! symmetry Operations for my crystal symmetry
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
axis = 0.0_pReal dQ = 0.0_pReal
angle = 0.0_pReal
rot = 0.0_pReal
! choose my symmetry operations according to my crystal symetry if (symmetryType < 1_pInt .or. symmetryType > 2_pInt) then
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
call IO_warning(700) call IO_warning(700)
return return
endif 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 dQ(1) = -1.0_pReal ! start with maximum misorientation angle
angle = 2.0_pReal * pi 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)
! apply symmetry operation to 1st orientation if (this_dQ(1) > dQ(1)) dQ = this_dQ ! store if misorientation angle is smaller (cos is larger) than previous one
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
enddo enddo
! convert angle to degrees
angle = angle * inDeg
endsubroutine endsubroutine
@ -681,6 +538,124 @@ endsubroutine
ENDFUNCTION 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 ! transposition of a 3x3 matrix
!************************************************************************** !**************************************************************************
@ -1046,15 +1021,15 @@ pure function math_transpose3x3(A)
!******************************************************************** !********************************************************************
! euclidic norm of a 3x1 vector ! euclidic norm of a 3x1 vector
!******************************************************************** !********************************************************************
pure function math_norm3(v3) pure function math_norm3(v)
use prec, only: pReal,pInt use prec, only: pReal,pInt
implicit none implicit none
real(pReal), dimension(3), intent(in) :: v3 real(pReal), dimension(3), intent(in) :: v
real(pReal) math_norm3 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 return
endfunction endfunction
@ -1275,6 +1250,138 @@ pure function math_transpose3x3(A)
ENDFUNCTION 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) ! rotation matrix from axis and angle (in radians)
!**************************************************************** !****************************************************************
@ -1314,38 +1421,6 @@ pure function math_transpose3x3(A)
ENDFUNCTION 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 ! 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 EB2(2,2)=1.0_pReal
EB3(3,3)=1.0_pReal EB3(3,3)=1.0_pReal
ELSE 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 arg=-Q/RHO/2.0_pReal
if(arg.GT.1) arg=1 if(arg.GT.1) arg=1
if(arg.LT.-1) arg=-1 if(arg.LT.-1) arg=-1
PHI=ACOS(arg) PHI=DACOS(arg)
Y1=2*RHO**(1.0_pReal/3.0_pReal)*COS(PHI/3.0_pReal) Y1=2*RHO**(1.0_pReal/3.0_pReal)*DCOS(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) 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)*COS(PHI/3.0_pReal+4.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 EW1=Y1-R/3.0_pReal
EW2=Y2-R/3.0_pReal EW2=Y2-R/3.0_pReal
EW3=Y3-R/3.0_pReal EW3=Y3-R/3.0_pReal