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

corrected (?) RodrigToR

This commit is contained in:
Philip Eisenlohr 2009-03-18 14:05:27 +00:00
parent 9007fbc73c
commit 121723ffa5
1 changed files with 18 additions and 61 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

79
trunk/math.f90
View File

@ -1213,18 +1213,25 @@ end subroutine
real(pReal) s,c real(pReal) s,c
integer(pInt) i integer(pInt) i
forall (i=1:3) axisNrm(i) = axis(i)/sqrt(dot_product(axis,axis)) forall (i=1:3) axisNrm(i) = axis(i)/dsqrt(math_mul3x3(axis,axis))
s = sin(omega) s = sin(omega)
c = cos(omega) c = cos(omega)
math_RodrigtoR(1,1) = (1.0_pReal - axisNrm(1)**2)*c + axisNrm(1)**2 c1 = 1.0_pReal - c
math_RodrigtoR(1,2) = axisNrm(1)*axisNrm(2)*(1.0_pReal - c) + axisNrm(3)*s
math_RodrigtoR(1,3) = axisNrm(1)*axisNrm(3)*(1.0_pReal - c) - axisNrm(2)*s ! formula taken from http://mathworld.wolfram.com/RodriguesRotationFormula.html
math_RodrigtoR(2,1) = axisNrm(1)*axisNrm(2)*(1.0_pReal - c) - axisNrm(3)*s
math_RodrigtoR(2,2) = (1.0_pReal - axisNrm(2)**2)*c + axisNrm(2)**2 math_RodrigtoR(1,1) = c + c1*axisNrm(1)**2
math_RodrigtoR(2,3) = axisNrm(2)*axisNrm(3)*(1.0_pReal - c) + axisNrm(1)*s math_RodrigtoR(1,2) = -s*axisNrm(3) + c1*axisNrm(1)*axisNrm(2)
math_RodrigtoR(3,1) = axisNrm(1)*axisNrm(3)*(1.0_pReal - c) + axisNrm(2)*s math_RodrigtoR(1,3) = s*axisNrm(2) + c1*axisNrm(1)*axisNrm(3)
math_RodrigtoR(3,2) = axisNrm(2)*axisNrm(3)*(1.0_pReal - c) - axisNrm(1)*s
math_RodrigtoR(3,3) = (1.0_pReal - axisNrm(3)**2)*c + axisNrm(3)**2 math_RodrigtoR(2,1) = s*axisNrm(3) + c1*axisNrm(2)*axisNrm(1)
math_RodrigtoR(2,2) = c + c1*axisNrm(2)**2
math_RodrigtoR(2,3) = -s*axisNrm(1) + c1*axisNrm(2)*axisNrm(3)
math_RodrigtoR(3,1) = -s*axisNrm(2) + c1*axisNrm(3)*axisNrm(1)
math_RodrigtoR(3,2) = s*axisNrm(1) + c1*axisNrm(3)*axisNrm(2)
math_RodrigtoR(3,3) = c + c1*axisNrm(3)**2
return return
END FUNCTION END FUNCTION
@ -1274,10 +1281,7 @@ end subroutine
real(pReal), dimension(3,3) :: r real(pReal), dimension(3,3) :: r
real(pReal) math_disorient, tr real(pReal) math_disorient, tr
!!$OMP CRITICAL (evilmatmul)
r = math_mul33x33(math_EulerToR(EulerB),transpose(math_EulerToR(EulerA))) r = math_mul33x33(math_EulerToR(EulerB),transpose(math_EulerToR(EulerA)))
!!$OMP END CRITICAL (evilmatmul)
tr = (r(1,1)+r(2,2)+r(3,3)-1.0_pReal)*0.4999999_pReal tr = (r(1,1)+r(2,2)+r(3,3)-1.0_pReal)*0.4999999_pReal
math_disorient = abs(0.5_pReal*pi-asin(tr)) math_disorient = abs(0.5_pReal*pi-asin(tr))
@ -1339,10 +1343,7 @@ endif
if (rnd(5) <= exp(-1.0_pReal*(math_disorient(origin,disturb)/scatter)**2)) exit if (rnd(5) <= exp(-1.0_pReal*(math_disorient(origin,disturb)/scatter)**2)) exit
enddo enddo
!!$OMP CRITICAL (evilmatmul)
math_sampleGaussOri = math_RtoEuler(math_mul33x33(math_EulerToR(disturb),math_EulerToR(center))) math_sampleGaussOri = math_RtoEuler(math_mul33x33(math_EulerToR(disturb),math_EulerToR(center)))
!!$OMP END CRITICAL (evilmatmul)
return return
@ -1417,10 +1418,7 @@ endif
pRot = math_RodrigtoR(axis,angle) pRot = math_RodrigtoR(axis,angle)
! ---# apply the three rotations #--- ! ---# apply the three rotations #---
!!$OMP CRITICAL (evilmatmul)
math_sampleFiberOri = math_RtoEuler(math_mul33x33(pRot,math_mul33x33(fRot,oRot))) math_sampleFiberOri = math_RtoEuler(math_mul33x33(pRot,math_mul33x33(fRot,oRot)))
!!$OMP END CRITICAL (evilmatmul)
return return
@ -1428,15 +1426,10 @@ math_sampleFiberOri = math_RtoEuler(math_mul33x33(pRot,math_mul33x33(fRot,oRot))
!******************************************************************** !********************************************************************
! symmetric Euler angles for given symmetry string ! symmetric Euler angles for given symmetry string
! 'triclinic' or '', 'monoclinic', 'orthotropic' ! 'triclinic' or '', 'monoclinic', 'orthotropic'
!******************************************************************** !********************************************************************
PURE FUNCTION math_symmetricEulers(sym,Euler) PURE FUNCTION math_symmetricEulers(sym,Euler)
use prec, only: pReal, pInt use prec, only: pReal, pInt
@ -1496,14 +1489,7 @@ math_sampleFiberOri = math_RtoEuler(math_mul33x33(pRot,math_mul33x33(fRot,oRot))
CALL math_spectral1(CE,EW1,EW2,EW3,EB1,EB2,EB3) CALL math_spectral1(CE,EW1,EW2,EW3,EB1,EB2,EB3)
U=DSQRT(EW1)*EB1+DSQRT(EW2)*EB2+DSQRT(EW3)*EB3 U=DSQRT(EW1)*EB1+DSQRT(EW2)*EB2+DSQRT(EW3)*EB3
call math_invert3x3(U,UI,det,error) call math_invert3x3(U,UI,det,error)
if (.not. error) then if (.not. error) R = math_mul33x33(fe,ui)
!!$OMP CRITICAL (evilmatmul)
R = math_mul33x33(fe,ui)
!!$OMP END CRITICAL (evilmatmul)
endif
return return
@ -1561,10 +1547,7 @@ math_sampleFiberOri = math_RtoEuler(math_mul33x33(pRot,math_mul33x33(fRot,oRot))
D3=1.0_pReal/(EW3-EW1)/(EW3-EW2) D3=1.0_pReal/(EW3-EW1)/(EW3-EW2)
M1=M-EW1*math_I3 M1=M-EW1*math_I3
M2=M-EW2*math_I3 M2=M-EW2*math_I3
!!$OMP CRITICAL (evilmatmul)
EB3=math_mul33x33(M1,M2)*D3 EB3=math_mul33x33(M1,M2)*D3
!!$OMP END CRITICAL (evilmatmul)
EB1=math_I3-EB3 EB1=math_I3-EB3
! both EB2 and EW2 are set to zero so that they do not ! both EB2 and EW2 are set to zero so that they do not
@ -1575,11 +1558,7 @@ math_sampleFiberOri = math_RtoEuler(math_mul33x33(pRot,math_mul33x33(fRot,oRot))
D1=1.0_pReal/(EW1-EW2)/(EW1-EW3) D1=1.0_pReal/(EW1-EW2)/(EW1-EW3)
M2=M-math_I3*EW2 M2=M-math_I3*EW2
M3=M-math_I3*EW3 M3=M-math_I3*EW3
!!$OMP CRITICAL (evilmatmul)
EB1=math_mul33x33(M2,M3)*D1 EB1=math_mul33x33(M2,M3)*D1
!!$OMP END CRITICAL (evilmatmul)
EB2=math_I3-EB1 EB2=math_I3-EB1
! both EB3 and EW3 are set to zero so that they do not ! both EB3 and EW3 are set to zero so that they do not
! contribute to U in PDECOMPOSITION ! contribute to U in PDECOMPOSITION
@ -1589,11 +1568,7 @@ math_sampleFiberOri = math_RtoEuler(math_mul33x33(pRot,math_mul33x33(fRot,oRot))
D2=1.0_pReal/(EW2-EW1)/(EW2-EW3) D2=1.0_pReal/(EW2-EW1)/(EW2-EW3)
M1=M-math_I3*EW1 M1=M-math_I3*EW1
M3=M-math_I3*EW3 M3=M-math_I3*EW3
!!$OMP CRITICAL (evilmatmul)
EB2=math_mul33x33(M1,M3)*D2 EB2=math_mul33x33(M1,M3)*D2
!!$OMP END CRITICAL (evilmatmul)
EB1=math_I3-EB2 EB1=math_I3-EB2
! both EB3 and EW3 are set to zero so that they do not ! both EB3 and EW3 are set to zero so that they do not
! contribute to U in PDECOMPOSITION ! contribute to U in PDECOMPOSITION
@ -1606,12 +1581,9 @@ math_sampleFiberOri = math_RtoEuler(math_mul33x33(pRot,math_mul33x33(fRot,oRot))
M1=M-EW1*math_I3 M1=M-EW1*math_I3
M2=M-EW2*math_I3 M2=M-EW2*math_I3
M3=M-EW3*math_I3 M3=M-EW3*math_I3
!!$OMP CRITICAL (evilmatmul)
EB1=math_mul33x33(M2,M3)*D1 EB1=math_mul33x33(M2,M3)*D1
EB2=math_mul33x33(M1,M3)*D2 EB2=math_mul33x33(M1,M3)*D2
EB3=math_mul33x33(M1,M2)*D3 EB3=math_mul33x33(M1,M2)*D3
!!$OMP END CRITICAL (evilmatmul)
END IF END IF
END IF END IF
@ -1865,19 +1837,14 @@ math_sampleFiberOri = math_RtoEuler(math_mul33x33(pRot,math_mul33x33(fRot,oRot))
else else
ndim_save = value(1) ndim_save = value(1)
end if end if
else if ( name(1:1) == 'S' .or. name(1:1) == 's' ) then else if ( name(1:1) == 'S' .or. name(1:1) == 's' ) then
seed = value(1) seed = value(1)
end if end if
! !
! Get ! Get
! !
else if ( action(1:1) == 'G' .or. action(1:1) == 'g' ) then else if ( action(1:1) == 'G' .or. action(1:1) == 'g' ) then
if ( name(1:1) == 'B' .or. name(1:1) == 'b' ) then if ( name(1:1) == 'B' .or. name(1:1) == 'b' ) then
if ( ndim /= ndim_save ) then if ( ndim /= ndim_save ) then
deallocate ( base ) deallocate ( base )
ndim_save = ndim ndim_save = ndim
@ -1886,27 +1853,19 @@ math_sampleFiberOri = math_RtoEuler(math_mul33x33(pRot,math_mul33x33(fRot,oRot))
base(i) = prime(i) base(i) = prime(i)
end do end do
end if end if
value(1:ndim_save) = base(1:ndim_save) value(1:ndim_save) = base(1:ndim_save)
else if ( name(1:1) == 'N' .or. name(1:1) == 'n' ) then else if ( name(1:1) == 'N' .or. name(1:1) == 'n' ) then
value(1) = ndim_save value(1) = ndim_save
else if ( name(1:1) == 'S' .or. name(1:1) == 's' ) then else if ( name(1:1) == 'S' .or. name(1:1) == 's' ) then
value(1) = seed value(1) = seed
end if end if
! !
! Increment ! Increment
! !
else if ( action(1:1) == 'I' .or. action(1:1) == 'i' ) then else if ( action(1:1) == 'I' .or. action(1:1) == 'i' ) then
if ( name(1:1) == 'S' .or. name(1:1) == 's' ) then if ( name(1:1) == 'S' .or. name(1:1) == 's' ) then
seed = seed + value(1) seed = seed + value(1)
end if end if
end if end if
return return
@ -2078,7 +2037,6 @@ math_sampleFiberOri = math_RtoEuler(math_mul33x33(pRot,math_mul33x33(fRot,oRot))
if ( any ( base(1:ndim) <= 1 ) ) then if ( any ( base(1:ndim) <= 1 ) ) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
write ( *, '(a)' ) ' ' write ( *, '(a)' ) ' '
write ( *, '(a)' ) 'I_TO_HALTON - Fatal error!' write ( *, '(a)' ) 'I_TO_HALTON - Fatal error!'
write ( *, '(a)' ) ' An input base BASE is <= 1!' write ( *, '(a)' ) ' An input base BASE is <= 1!'
@ -2087,7 +2045,6 @@ math_sampleFiberOri = math_RtoEuler(math_mul33x33(pRot,math_mul33x33(fRot,oRot))
end do end do
call flush(6) call flush(6)
!$OMP END CRITICAL (write2out) !$OMP END CRITICAL (write2out)
stop stop
end if end if