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

following default style

This commit is contained in:
Martin Diehl 2022-06-15 20:15:35 +02:00
parent 0b8473d7ec
commit 074cdf5dff
1 changed files with 14 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

28
src/rotations.f90
View File

@ -182,7 +182,7 @@ subroutine fromEulers(self,eu,degrees)
Eulers = eu Eulers = eu
else else
Eulers = merge(eu*INRAD,eu,degrees) Eulers = merge(eu*INRAD,eu,degrees)
endif end if
if (any(Eulers<0.0_pReal) .or. any(Eulers>TAU) .or. Eulers(2) > PI) & if (any(Eulers<0.0_pReal) .or. any(Eulers>TAU) .or. Eulers(2) > PI) &
call IO_error(402,ext_msg='fromEulers') call IO_error(402,ext_msg='fromEulers')
@ -206,14 +206,14 @@ subroutine fromAxisAngle(self,ax,degrees,P)
angle = ax(4) angle = ax(4)
else else
angle = merge(ax(4)*INRAD,ax(4),degrees) angle = merge(ax(4)*INRAD,ax(4),degrees)
endif end if
if (.not. present(P)) then if (.not. present(P)) then
axis = ax(1:3) axis = ax(1:3)
else else
axis = ax(1:3) * merge(-1.0_pReal,1.0_pReal,P == 1) axis = ax(1:3) * merge(-1.0_pReal,1.0_pReal,P == 1)
if(abs(P) /= 1) call IO_error(402,ext_msg='fromAxisAngle (P)') if(abs(P) /= 1) call IO_error(402,ext_msg='fromAxisAngle (P)')
endif end if
if(dNeq(norm2(axis),1.0_pReal) .or. angle < 0.0_pReal .or. angle > PI) & if(dNeq(norm2(axis),1.0_pReal) .or. angle < 0.0_pReal .or. angle > PI) &
call IO_error(402,ext_msg='fromAxisAngle') call IO_error(402,ext_msg='fromAxisAngle')
@ -284,7 +284,7 @@ pure function rotVector(self,v,active) result(vRot)
passive = .not. active passive = .not. active
else else
passive = .true. passive = .true.
endif end if
if (dEq0(norm2(v))) then if (dEq0(norm2(v))) then
vRot = v vRot = v
@ -294,7 +294,7 @@ pure function rotVector(self,v,active) result(vRot)
multiplyQuaternion(conjugateQuaternion(self%q), multiplyQuaternion(v_normed, self%q)), & multiplyQuaternion(conjugateQuaternion(self%q), multiplyQuaternion(v_normed, self%q)), &
passive) passive)
vRot = q(2:4)*norm2(v) vRot = q(2:4)*norm2(v)
endif end if
end function rotVector end function rotVector
@ -318,7 +318,7 @@ pure function rotTensor2(self,T,active) result(tRot)
passive = .not. active passive = .not. active
else else
passive = .true. passive = .true.
endif end if
tRot = merge(matmul(matmul(self%asMatrix(),T),transpose(self%asMatrix())), & tRot = merge(matmul(matmul(self%asMatrix(),T),transpose(self%asMatrix())), &
matmul(matmul(transpose(self%asMatrix()),T),self%asMatrix()), & matmul(matmul(transpose(self%asMatrix()),T),self%asMatrix()), &
@ -347,14 +347,14 @@ pure function rotTensor4(self,T,active) result(tRot)
R = merge(transpose(self%asMatrix()),self%asMatrix(),active) R = merge(transpose(self%asMatrix()),self%asMatrix(),active)
else else
R = self%asMatrix() R = self%asMatrix()
endif end if
tRot = 0.0_pReal tRot = 0.0_pReal
do i = 1,3;do j = 1,3;do k = 1,3;do l = 1,3 do i = 1,3;do j = 1,3;do k = 1,3;do l = 1,3
do m = 1,3;do n = 1,3;do o = 1,3;do p = 1,3 do m = 1,3;do n = 1,3;do o = 1,3;do p = 1,3
tRot(i,j,k,l) = tRot(i,j,k,l) & tRot(i,j,k,l) = tRot(i,j,k,l) &
+ R(i,m) * R(j,n) * R(k,o) * R(l,p) * T(m,n,o,p) + R(i,m) * R(j,n) * R(k,o) * R(l,p) * T(m,n,o,p)
enddo; enddo; enddo; enddo; enddo; enddo; enddo; enddo end do; end do; end do; end do; end do; end do; end do; end do
end function rotTensor4 end function rotTensor4
@ -379,7 +379,7 @@ pure function rotStiffness(self,C,active) result(cRot)
R = merge(transpose(self%asMatrix()),self%asMatrix(),active) R = merge(transpose(self%asMatrix()),self%asMatrix(),active)
else else
R = self%asMatrix() R = self%asMatrix()
endif end if
M = reshape([R(1,1)**2, R(2,1)**2, R(3,1)**2, & M = reshape([R(1,1)**2, R(2,1)**2, R(3,1)**2, &
R(2,1)*R(3,1), R(1,1)*R(3,1), R(1,1)*R(2,1), & R(2,1)*R(3,1), R(1,1)*R(3,1), R(1,1)*R(2,1), &
@ -468,7 +468,7 @@ pure function qu2eu(qu) result(eu)
eu = [atan2((-P*qu(1)*qu(3)+qu(2)*qu(4))*chi, (-P*qu(1)*qu(2)-qu(3)*qu(4))*chi ), & eu = [atan2((-P*qu(1)*qu(3)+qu(2)*qu(4))*chi, (-P*qu(1)*qu(2)-qu(3)*qu(4))*chi ), &
atan2( 2.0_pReal*chi, q03-q12 ), & atan2( 2.0_pReal*chi, q03-q12 ), &
atan2(( P*qu(1)*qu(3)+qu(2)*qu(4))*chi, (-P*qu(1)*qu(2)+qu(3)*qu(4))*chi )] atan2(( P*qu(1)*qu(3)+qu(2)*qu(4))*chi, (-P*qu(1)*qu(2)+qu(3)*qu(4))*chi )]
endif degenerated end if degenerated
where(sign(1.0_pReal,eu)<0.0_pReal) eu = mod(eu+TAU,[TAU,PI,TAU]) where(sign(1.0_pReal,eu)<0.0_pReal) eu = mod(eu+TAU,[TAU,PI,TAU])
end function qu2eu end function qu2eu
@ -526,8 +526,8 @@ pure function om2qu(om) result(qu)
else else
s = 2.0_pReal * sqrt( 1.0_pReal + om(3,3) - om(1,1) - om(2,2) ) s = 2.0_pReal * sqrt( 1.0_pReal + om(3,3) - om(1,1) - om(2,2) )
qu = [ (om(2,1) - om(1,2)) /s,(om(1,3) + om(3,1)) / s,(om(2,3) + om(3,2)) / s,0.25_pReal * s] qu = [ (om(2,1) - om(1,2)) /s,(om(1,3) + om(3,1)) / s,(om(2,3) + om(3,2)) / s,0.25_pReal * s]
endif end if
endif end if
if(sign(1.0_pReal,qu(1))<0.0_pReal) qu =-1.0_pReal * qu if(sign(1.0_pReal,qu(1))<0.0_pReal) qu =-1.0_pReal * qu
qu(2:4) = merge(qu(2:4),qu(2:4)*P,dEq0(qu(2:4))) qu(2:4) = merge(qu(2:4),qu(2:4)*P,dEq0(qu(2:4)))
qu = qu/norm2(qu) qu = qu/norm2(qu)
@ -593,7 +593,7 @@ function om2ax(om) result(ax)
ax(1:3) = VR(1:3,i) ax(1:3) = VR(1:3,i)
where ( dNeq0([om(2,3)-om(3,2), om(3,1)-om(1,3), om(1,2)-om(2,1)])) & where ( dNeq0([om(2,3)-om(3,2), om(3,1)-om(1,3), om(1,2)-om(2,1)])) &
ax(1:3) = sign(ax(1:3),-P *[om(2,3)-om(3,2), om(3,1)-om(1,3), om(1,2)-om(2,1)]) ax(1:3) = sign(ax(1:3),-P *[om(2,3)-om(3,2), om(3,1)-om(1,3), om(1,2)-om(2,1)])
endif end if
end function om2ax end function om2ax
@ -826,7 +826,7 @@ subroutine selfTest()
cos(TAU*x(2))*B,& cos(TAU*x(2))*B,&
sin(TAU*x(1))*A] sin(TAU*x(1))*A]
if(qu(1)<0.0_pReal) qu = qu * (-1.0_pReal) if(qu(1)<0.0_pReal) qu = qu * (-1.0_pReal)
endif end if
if(.not. quaternion_equal(om2qu(qu2om(qu)),qu)) error stop 'om2qu2om' if(.not. quaternion_equal(om2qu(qu2om(qu)),qu)) error stop 'om2qu2om'