functionality to rotate 4-tensor

This commit is contained in:
Martin Diehl 2019-09-20 16:39:33 -07:00
parent 8d66136f9c
commit 6c0e92d5c1
3 changed files with 68 additions and 37 deletions

View File

@ -68,7 +68,7 @@ contains
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief map from 3D cubic grid to 3D ball
!--------------------------------------------------------------------------
function Lambert_CubeToBall(cube) result(ball)
pure function Lambert_CubeToBall(cube) result(ball)
real(pReal), intent(in), dimension(3) :: cube
real(pReal), dimension(3) :: ball, LamXYZ, XYZ

View File

@ -1428,12 +1428,12 @@ end function math_areaTriangle
!--------------------------------------------------------------------------------------------------
!> @brief rotate 33 tensor forward
!--------------------------------------------------------------------------------------------------
pure function math_rotate_forward33(tensor,rot_tensor)
pure function math_rotate_forward33(tensor,R)
real(pReal), dimension(3,3) :: math_rotate_forward33
real(pReal), dimension(3,3), intent(in) :: tensor, rot_tensor
real(pReal), dimension(3,3) :: math_rotate_forward33
real(pReal), dimension(3,3), intent(in) :: tensor, R
math_rotate_forward33 = matmul(rot_tensor,matmul(tensor,transpose(rot_tensor)))
math_rotate_forward33 = matmul(R,matmul(tensor,transpose(R)))
end function math_rotate_forward33
@ -1441,12 +1441,12 @@ end function math_rotate_forward33
!--------------------------------------------------------------------------------------------------
!> @brief rotate 33 tensor backward
!--------------------------------------------------------------------------------------------------
pure function math_rotate_backward33(tensor,rot_tensor)
pure function math_rotate_backward33(tensor,R)
real(pReal), dimension(3,3) :: math_rotate_backward33
real(pReal), dimension(3,3), intent(in) :: tensor, rot_tensor
real(pReal), dimension(3,3) :: math_rotate_backward33
real(pReal), dimension(3,3), intent(in) :: tensor, R
math_rotate_backward33 = matmul(transpose(rot_tensor),matmul(tensor,rot_tensor))
math_rotate_backward33 = matmul(transpose(R),matmul(tensor,R))
end function math_rotate_backward33
@ -1454,19 +1454,18 @@ end function math_rotate_backward33
!--------------------------------------------------------------------------------------------------
!> @brief rotate 3333 tensor C'_ijkl=g_im*g_jn*g_ko*g_lp*C_mnop
!--------------------------------------------------------------------------------------------------
pure function math_rotate_forward3333(tensor,rot_tensor)
pure function math_rotate_forward3333(tensor,R)
real(pReal), dimension(3,3,3,3) :: math_rotate_forward3333
real(pReal), dimension(3,3), intent(in) :: rot_tensor
real(pReal), dimension(3,3,3,3) :: math_rotate_forward3333
real(pReal), dimension(3,3), intent(in) :: R
real(pReal), dimension(3,3,3,3), intent(in) :: tensor
integer :: i,j,k,l,m,n,o,p
math_rotate_forward3333 = 0.0_pReal
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
math_rotate_forward3333(i,j,k,l) &
= math_rotate_forward3333(i,j,k,l) &
+ rot_tensor(i,m) * rot_tensor(j,n) * rot_tensor(k,o) * rot_tensor(l,p) * tensor(m,n,o,p)
math_rotate_forward3333(i,j,k,l) = math_rotate_forward3333(i,j,k,l) &
+ R(i,m) * R(j,n) * R(k,o) * R(l,p) * tensor(m,n,o,p)
enddo; enddo; enddo; enddo; enddo; enddo; enddo; enddo
end function math_rotate_forward3333
@ -1482,8 +1481,8 @@ real(pReal) pure elemental function math_clip(a, left, right)
real(pReal), intent(in), optional :: left, right
math_clip = a
if (present(left)) math_clip = max(left,math_clip)
if (present(right)) math_clip = min(right,math_clip)
if (present(left)) math_clip = max(left,math_clip)
if (present(right)) math_clip = min(right,math_clip)
if (present(left) .and. present(right)) &
math_clip = merge (IEEE_value(1.0_pReal,IEEE_quiet_NaN),math_clip, left>right)

View File

@ -70,9 +70,10 @@ module rotations
!------------------------------------------
procedure, private :: rotRot__
generic, public :: operator(*) => rotRot__
generic, public :: rotate => rotVector,rotTensor2,rotTensor4
procedure, public :: rotVector
procedure, public :: rotTensor2
!procedure, public :: rotTensor4
procedure, public :: rotTensor4
procedure, public :: misorientation
end type rotation
@ -209,12 +210,12 @@ end subroutine fromRotationMatrix
!> @brief: Rotate a rotation
!> ToDo: completly untested
!---------------------------------------------------------------------------------------------------
function rotRot__(self,r) result(rRot)
pure elemental function rotRot__(self,R) result(rRot)
type(rotation) :: rRot
class(rotation), intent(in) :: self,r
class(rotation), intent(in) :: self,R
rRot = rotation(self%q*r%q)
rRot = rotation(self%q*R%q)
end function rotRot__
@ -223,7 +224,7 @@ end function rotRot__
!> @author Marc De Graef, Carnegie Mellon University
!> @brief rotate a vector passively (default) or actively
!---------------------------------------------------------------------------------------------------
function rotVector(self,v,active) result(vRot)
pure function rotVector(self,v,active) result(vRot)
real(pReal), dimension(3) :: vRot
class(rotation), intent(in) :: self
@ -257,14 +258,14 @@ end function rotVector
!---------------------------------------------------------------------------------------------------
!> @author Marc De Graef, Carnegie Mellon University
!> @brief rotate a second rank tensor passively (default) or actively
!> @brief rotate a rank-2 tensor passively (default) or actively
!> @details: rotation is based on rotation matrix
!---------------------------------------------------------------------------------------------------
function rotTensor2(self,m,active) result(mRot)
pure function rotTensor2(self,T,active) result(tRot)
real(pReal), dimension(3,3) :: mRot
real(pReal), dimension(3,3) :: tRot
class(rotation), intent(in) :: self
real(pReal), intent(in), dimension(3,3) :: m
real(pReal), intent(in), dimension(3,3) :: T
logical, intent(in), optional :: active
logical :: passive
@ -276,18 +277,49 @@ function rotTensor2(self,m,active) result(mRot)
endif
if (passive) then
mRot = matmul(matmul(self%asRotationMatrix(),m),transpose(self%asRotationMatrix()))
tRot = matmul(matmul(self%asRotationMatrix(),T),transpose(self%asRotationMatrix()))
else
mRot = matmul(matmul(transpose(self%asRotationMatrix()),m),self%asRotationMatrix())
tRot = matmul(matmul(transpose(self%asRotationMatrix()),T),self%asRotationMatrix())
endif
end function rotTensor2
!---------------------------------------------------------------------------------------------------
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief rotate a rank-4 tensor passively (default) or actively
!> @details: rotation is based on rotation matrix
!! ToDo: Need to check active/passive !!!
!---------------------------------------------------------------------------------------------------
pure function rotTensor4(self,T,active) result(tRot)
real(pReal), dimension(3,3,3,3) :: tRot
class(rotation), intent(in) :: self
real(pReal), intent(in), dimension(3,3,3,3) :: T
logical, intent(in), optional :: active
real(pReal), dimension(3,3) :: R
integer :: i,j,k,l,m,n,o,p
if (present(active)) then
R = merge(transpose(self%asRotationMatrix()),self%asRotationMatrix(),active)
else
R = self%asRotationMatrix()
endif
tRot = 0.0_pReal
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
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)
enddo; enddo; enddo; enddo; enddo; enddo; enddo; enddo
end function rotTensor4
!---------------------------------------------------------------------------------------------------
!> @brief misorientation
!---------------------------------------------------------------------------------------------------
function misorientation(self,other)
pure elemental function misorientation(self,other)
type(rotation) :: misorientation
class(rotation), intent(in) :: self, other
@ -436,7 +468,7 @@ end function qu2ho
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert unit quaternion to cubochoric
!---------------------------------------------------------------------------------------------------
function qu2cu(qu) result(cu)
pure function qu2cu(qu) result(cu)
real(pReal), intent(in), dimension(4) :: qu
real(pReal), dimension(3) :: cu
@ -922,7 +954,7 @@ end function ro2ho
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert Rodrigues vector to cubochoric
!---------------------------------------------------------------------------------------------------
function ro2cu(ro) result(cu)
pure function ro2cu(ro) result(cu)
real(pReal), intent(in), dimension(4) :: ro
real(pReal), dimension(3) :: cu
@ -1032,7 +1064,7 @@ end function ho2ro
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert homochoric to cubochoric
!---------------------------------------------------------------------------------------------------
function ho2cu(ho) result(cu)
pure function ho2cu(ho) result(cu)
real(pReal), intent(in), dimension(3) :: ho
real(pReal), dimension(3) :: cu
@ -1046,7 +1078,7 @@ end function ho2cu
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert cubochoric to unit quaternion
!---------------------------------------------------------------------------------------------------
function cu2qu(cu) result(qu)
pure function cu2qu(cu) result(qu)
real(pReal), intent(in), dimension(3) :: cu
real(pReal), dimension(4) :: qu
@ -1060,7 +1092,7 @@ end function cu2qu
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert cubochoric to rotation matrix
!---------------------------------------------------------------------------------------------------
function cu2om(cu) result(om)
pure function cu2om(cu) result(om)
real(pReal), intent(in), dimension(3) :: cu
real(pReal), dimension(3,3) :: om
@ -1074,7 +1106,7 @@ end function cu2om
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert cubochoric to Euler angles
!---------------------------------------------------------------------------------------------------
function cu2eu(cu) result(eu)
pure function cu2eu(cu) result(eu)
real(pReal), intent(in), dimension(3) :: cu
real(pReal), dimension(3) :: eu
@ -1102,7 +1134,7 @@ end function cu2ax
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert cubochoric to Rodrigues vector
!---------------------------------------------------------------------------------------------------
function cu2ro(cu) result(ro)
pure function cu2ro(cu) result(ro)
real(pReal), intent(in), dimension(3) :: cu
real(pReal), dimension(4) :: ro
@ -1116,7 +1148,7 @@ end function cu2ro
!> @author Marc De Graef, Carnegie Mellon University
!> @brief convert cubochoric to homochoric
!---------------------------------------------------------------------------------------------------
function cu2ho(cu) result(ho)
pure function cu2ho(cu) result(ho)
real(pReal), intent(in), dimension(3) :: cu
real(pReal), dimension(3) :: ho