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

clearer naming for vector <-> tensor conversion 

and a bunch of other improvements
This commit is contained in:
Martin Diehl 2019-01-14 07:27:18 +01:00
parent c9977e2b78
commit 01a2fffd3c
1 changed files with 266 additions and 186 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

452
src/math.f90
View File

@ -24,24 +24,6 @@ module math
0.0_pReal,0.0_pReal,1.0_pReal &
],[3,3]) !< 3x3 Identity
! ToDo MD: Our naming scheme is a little bit odd: We use essentially the re-ordering according to Nye
! (convenient because Abaqus and Marc want to have 12 on position 4)
! but weight the shear components according to Mandel (convenient for matrix multiplications)
! I suggest to keep Voigt3333to66 (required for reading in elasticity matrices) but rename
! mapMandel to mapNye, math_MandelXtoY to math_XtoY and math_PlainXtoY to math_XtoY.
! It is then clear that math_33to9 just reorders and math_33to6 does the "DAMASK conversion"
! without leaving the impression that it follows any established convention
integer(pInt), dimension (2,6), parameter, private :: &
mapMandel = reshape([&
1_pInt,1_pInt, &
2_pInt,2_pInt, &
3_pInt,3_pInt, &
1_pInt,2_pInt, &
2_pInt,3_pInt, &
1_pInt,3_pInt &
],[2,6]) !< arrangement in Mandel notation. Differs from https://en.wikipedia.org/wiki/Voigt_notation#Mandel_notation
real(pReal), dimension(6), parameter, private :: &
nrmMandel = [&
1.0_pReal, 1.0_pReal, 1.0_pReal, &
@ -52,6 +34,16 @@ module math
1.0_pReal, 1.0_pReal, 1.0_pReal, &
1.0_pReal/sqrt(2.0_pReal), 1.0_pReal/sqrt(2.0_pReal), 1.0_pReal/sqrt(2.0_pReal) ] !< weighting for Mandel notation (backward)
integer(pInt), dimension (2,6), parameter, private :: &
mapNye = reshape([&
1_pInt,1_pInt, &
2_pInt,2_pInt, &
3_pInt,3_pInt, &
1_pInt,2_pInt, &
2_pInt,3_pInt, &
1_pInt,3_pInt &
],[2,6]) !< arrangement in Nye notation.
integer(pInt), dimension (2,6), parameter, private :: &
mapVoigt = reshape([&
1_pInt,1_pInt, &
@ -62,10 +54,6 @@ module math
1_pInt,2_pInt &
],[2,6]) !< arrangement in Voigt notation
real(pReal), dimension(6), parameter, private :: &
nrmVoigt = 1.0_pReal, & !< weighting for Voigt notation (forward)
invnrmVoigt = 1.0_pReal !< weighting for Voigt notation (backward)
integer(pInt), dimension (2,9), parameter, private :: &
mapPlain = reshape([&
1_pInt,1_pInt, &
@ -78,6 +66,56 @@ module math
3_pInt,2_pInt, &
3_pInt,3_pInt &
],[2,9]) !< arrangement in Plain notation
!--------------------------------------------------------------------------------------------------
! Provide deprecated names for compatibility
! ToDo MD: Our naming scheme was a little bit odd: We use essentially the re-ordering according to Nye
! (convenient because Abaqus and Marc want to have 12 on position 4)
! but weight the shear components according to Mandel (convenient for matrix multiplications)
interface math_Plain33to9
module procedure math_33to9
end interface math_Plain33to9
interface math_Plain9to33
module procedure math_9to33
end interface math_Plain9to33
interface math_Mandel33to6
module procedure math_33to6
end interface math_Mandel33to6
interface math_Mandel6to33
module procedure math_6to33
end interface math_Mandel6to33
interface math_Plain3333to99
module procedure math_3333to99
end interface math_Plain3333to99
interface math_Plain99to3333
module procedure math_99to3333
end interface math_Plain99to3333
interface math_Mandel3333to66
module procedure math_3333to66
end interface math_Mandel3333to66
interface math_Mandel66to3333
module procedure math_66to3333
end interface math_Mandel66to3333
public :: &
math_Plain33to9, &
math_Plain9to33, &
math_Mandel33to6, &
math_Mandel6to33, &
math_Plain3333to99, &
math_Plain99to3333, &
math_Mandel3333to66, &
math_Mandel66to3333
!---------------------------------------------------------------------------------------------------
public :: &
math_init, &
@ -116,16 +154,14 @@ module math
math_equivStress33, &
math_trace33, &
math_det33, &
math_Plain33to9, &
math_Plain9to33, &
math_Mandel33to6, &
math_Mandel6to33, &
math_Plain3333to99, &
math_Plain99to3333, &
math_Mandel66toPlain66, &
math_Plain66toMandel66, &
math_Mandel3333to66, &
math_Mandel66to3333, &
math_33to9, &
math_9to33, &
math_33to6, &
math_6to33, &
math_3333to99, &
math_99to3333, &
math_3333to66, &
math_66to3333, &
math_Voigt66to3333, &
math_qRand, &
math_qMul, &
@ -437,9 +473,13 @@ pure function math_identity4th(dimen)
integer(pInt), intent(in) :: dimen !< tensor dimension
integer(pInt) :: i,j,k,l
real(pReal), dimension(dimen,dimen,dimen,dimen) :: math_identity4th
real(pReal), dimension(dimen,dimen) :: identity2nd
forall (i=1_pInt:dimen,j=1_pInt:dimen,k=1_pInt:dimen,l=1_pInt:dimen) math_identity4th(i,j,k,l) = &
0.5_pReal*(math_I3(i,k)*math_I3(j,l)+math_I3(i,l)*math_I3(j,k))
identity2nd = math_identity2nd(dimen)
do concurrent(i=1_pInt:dimen,j=1_pInt:dimen,k=1_pInt:dimen,l=1_pInt:dimen)
math_identity4th(i,j,k,l) &
= 0.5_pReal*(identity2nd(i,k)*identity2nd(j,l)+identity2nd(i,l)*identity2nd(j,k))
enddo
end function math_identity4th
@ -564,7 +604,9 @@ real(pReal) pure function math_mul33xx33(A,B)
integer(pInt) :: i,j
real(pReal), dimension(3,3) :: C
forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt) C(i,j) = A(i,j) * B(i,j)
do concurrent(i=1_pInt:3_pInt,j=1_pInt:3_pInt)
C(i,j) = A(i,j) * B(i,j)
enddo
math_mul33xx33 = sum(C)
end function math_mul33xx33
@ -581,9 +623,10 @@ pure function math_mul3333xx33(A,B)
real(pReal), dimension(3,3), intent(in) :: B
integer(pInt) :: i,j
forall(i = 1_pInt:3_pInt,j = 1_pInt:3_pInt) &
do concurrent(i = 1_pInt:3_pInt,j = 1_pInt:3_pInt)
math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3))
enddo
end function math_mul3333xx33
@ -598,8 +641,9 @@ pure function math_mul3333xx3333(A,B)
real(pReal), dimension(3,3,3,3), intent(in) :: B
real(pReal), dimension(3,3,3,3) :: math_mul3333xx3333
forall(i = 1_pInt:3_pInt,j = 1_pInt:3_pInt, k = 1_pInt:3_pInt, l= 1_pInt:3_pInt) &
do concurrent(i = 1_pInt:3_pInt,j = 1_pInt:3_pInt, k = 1_pInt:3_pInt, l= 1_pInt:3_pInt)
math_mul3333xx3333(i,j,k,l) = sum(A(i,j,1:3,1:3)*B(1:3,1:3,k,l))
enddo
end function math_mul3333xx3333
@ -614,8 +658,9 @@ pure function math_mul33x33(A,B)
real(pReal), dimension(3,3), intent(in) :: A,B
integer(pInt) :: i,j
forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt) &
do concurrent(i=1_pInt:3_pInt,j=1_pInt:3_pInt)
math_mul33x33(i,j) = A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j)
enddo
end function math_mul33x33
@ -630,9 +675,10 @@ pure function math_mul66x66(A,B)
real(pReal), dimension(6,6), intent(in) :: A,B
integer(pInt) :: i,j
forall (i=1_pInt:6_pInt,j=1_pInt:6_pInt) math_mul66x66(i,j) = &
A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j) + &
A(i,4)*B(4,j) + A(i,5)*B(5,j) + A(i,6)*B(6,j)
do concurrent(i=1_pInt:6_pInt,j=1_pInt:6_pInt)
math_mul66x66(i,j) = A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j) &
+ A(i,4)*B(4,j) + A(i,5)*B(5,j) + A(i,6)*B(6,j)
enddo
end function math_mul66x66
@ -647,10 +693,11 @@ pure function math_mul99x99(A,B)
real(pReal), dimension(9,9), intent(in) :: A,B
integer(pInt) i,j
forall (i=1_pInt:9_pInt,j=1_pInt:9_pInt) math_mul99x99(i,j) = &
A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j) + &
A(i,4)*B(4,j) + A(i,5)*B(5,j) + A(i,6)*B(6,j) + &
A(i,7)*B(7,j) + A(i,8)*B(8,j) + A(i,9)*B(9,j)
do concurrent(i=1_pInt:9_pInt,j=1_pInt:9_pInt)
math_mul99x99(i,j) = A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j) &
+ A(i,4)*B(4,j) + A(i,5)*B(5,j) + A(i,6)*B(6,j) &
+ A(i,7)*B(7,j) + A(i,8)*B(8,j) + A(i,9)*B(9,j)
enddo
end function math_mul99x99
@ -747,8 +794,8 @@ end function math_transpose33
!--------------------------------------------------------------------------------------------------
!> @brief Cramer inversion of 33 matrix (function)
! direct Cramer inversion of matrix A.
! returns all zeroes if not possible, i.e. if det close to zero
!> @details Direct Cramer inversion of matrix A. Returns all zeroes if not possible, i.e.
! if determinant is close to zero
!--------------------------------------------------------------------------------------------------
pure function math_inv33(A)
use prec, only: &
@ -784,9 +831,9 @@ end function math_inv33
!--------------------------------------------------------------------------------------------------
!> @brief Cramer inversion of 33 matrix (subroutine)
! direct Cramer inversion of matrix A.
! also returns determinant
! returns error if not possible, i.e. if det close to zero
!> @details Direct Cramer inversion of matrix A. Also returns determinant
! Returns an error if not possible, i.e. if determinant is close to zero
! ToDo: Output arguments should be first
!--------------------------------------------------------------------------------------------------
pure subroutine math_invert33(A, InvA, DetA, error)
use prec, only: &
@ -843,20 +890,38 @@ function math_invSym3333(A)
dgetrf, &
dgetri
temp66_real = math_Mandel3333to66(A)
temp66_real = math_3333to66(A)
call dgetrf(6,6,temp66_real,6,ipiv6,ierr)
call dgetri(6,temp66_real,6,ipiv6,work6,6,ierr)
if (ierr == 0_pInt) then
math_invSym3333 = math_Mandel66to3333(temp66_real)
math_invSym3333 = math_66to3333(temp66_real)
else
call IO_error(400_pInt, ext_msg = 'math_invSym3333')
endif
end function math_invSym3333
!--------------------------------------------------------------------------------------------------
!> @brief invert quadratic matrix of arbitrary dimension
! ToDo: replaces math_invert
!--------------------------------------------------------------------------------------------------
subroutine math_invert2(InvA, error, A)
implicit none
real(pReal), dimension(:,:), intent(in) :: A
real(pReal), dimension(size(A,1),size(A,1)), intent(out) :: invA
logical, intent(out) :: error
call math_invert(size(A,1), A, InvA, error)
end subroutine math_invert2
!--------------------------------------------------------------------------------------------------
!> @brief invert matrix of arbitrary dimension
! ToDo: Wrong order of arguments and superfluous myDim argument.
! Use math_invert2 instead
!--------------------------------------------------------------------------------------------------
subroutine math_invert(myDim,A, InvA, error)
@ -961,15 +1026,14 @@ pure function math_equivStrain33(m)
real(pReal), dimension(3,3), intent(in) :: m
real(pReal), dimension(3) :: e,s
real(pReal) :: math_equivStrain33
real(pReal), parameter :: TWOTHIRD = 2.0_pReal/3.0_pReal
e = [2.0_pReal*m(1,1)-m(2,2)-m(3,3), &
2.0_pReal*m(2,2)-m(3,3)-m(1,1), &
2.0_pReal*m(3,3)-m(1,1)-m(2,2)]/3.0_pReal
s = [m(1,2),m(2,3),m(1,3)]*2.0_pReal
math_equivStrain33 = TWOTHIRD*(1.50_pReal*(sum(e**2.0_pReal)) + &
0.75_pReal*(sum(s**2.0_pReal)))**(0.5_pReal)
math_equivStrain33 = 2.0_pReal/3.0_pReal &
* (1.50_pReal*(sum(e**2.0_pReal))+ 0.75_pReal*(sum(s**2.0_pReal)))**(0.5_pReal)
end function math_equivStrain33
@ -1041,168 +1105,187 @@ end function math_detSym33
!--------------------------------------------------------------------------------------------------
!> @brief convert 33 matrix into vector 9
!--------------------------------------------------------------------------------------------------
pure function math_Plain33to9(m33)
pure function math_33to9(m33)
implicit none
real(pReal), dimension(9) :: math_Plain33to9
real(pReal), dimension(3,3), intent(in) :: m33
integer(pInt) :: i
forall (i=1_pInt:9_pInt) math_Plain33to9(i) = m33(mapPlain(1,i),mapPlain(2,i))
end function math_Plain33to9
!--------------------------------------------------------------------------------------------------
!> @brief convert Plain 9 back to 33 matrix
!--------------------------------------------------------------------------------------------------
pure function math_Plain9to33(v9)
implicit none
real(pReal), dimension(3,3) :: math_Plain9to33
real(pReal), dimension(9), intent(in) :: v9
integer(pInt) :: i
forall (i=1_pInt:9_pInt) math_Plain9to33(mapPlain(1,i),mapPlain(2,i)) = v9(i)
end function math_Plain9to33
!--------------------------------------------------------------------------------------------------
!> @brief convert symmetric 33 matrix into Mandel vector 6
!--------------------------------------------------------------------------------------------------
pure function math_Mandel33to6(m33)
implicit none
real(pReal), dimension(6) :: math_Mandel33to6
real(pReal), dimension(9) :: math_33to9
real(pReal), dimension(3,3), intent(in) :: m33
integer(pInt) :: i
forall (i=1_pInt:6_pInt) math_Mandel33to6(i) = nrmMandel(i)*m33(mapMandel(1,i),mapMandel(2,i))
forall (i=1_pInt:9_pInt) math_33to9(i) = m33(mapPlain(1,i),mapPlain(2,i))
end function math_Mandel33to6
end function math_33to9
!--------------------------------------------------------------------------------------------------
!> @brief convert Mandel 6 back to symmetric 33 matrix
!> @brief convert 9 vector into 33 matrix
!--------------------------------------------------------------------------------------------------
pure function math_Mandel6to33(v6)
pure function math_9to33(v9)
implicit none
real(pReal), dimension(6), intent(in) :: v6
real(pReal), dimension(3,3) :: math_Mandel6to33
real(pReal), dimension(3,3) :: math_9to33
real(pReal), dimension(9), intent(in) :: v9
integer(pInt) :: i
forall (i=1_pInt:6_pInt)
math_Mandel6to33(mapMandel(1,i),mapMandel(2,i)) = invnrmMandel(i)*v6(i)
math_Mandel6to33(mapMandel(2,i),mapMandel(1,i)) = invnrmMandel(i)*v6(i)
end forall
forall (i=1_pInt:9_pInt) math_9to33(mapPlain(1,i),mapPlain(2,i)) = v9(i)
end function math_Mandel6to33
end function math_9to33
!--------------------------------------------------------------------------------------------------
!> @brief convert 3333 tensor into plain matrix 99
!> @brief convert symmetric 33 matrix into 6 vector
!> @details Weighted conversion (default) rearranges according to Nye and weights shear
! components according to Mandel. Advisable for matrix operations.
! Unweighted conversion only changes order according to Nye
!--------------------------------------------------------------------------------------------------
pure function math_Plain3333to99(m3333)
pure function math_33to6(m33,weighted)
implicit none
real(pReal), dimension(6) :: math_33to6
real(pReal), dimension(3,3), intent(in) :: m33
logical, optional, intent(in) :: weighted
real(pReal), dimension(6) :: w
integer(pInt) :: i
if(present(weighted)) then
w = merge(nrmMandel,1.0_pReal,weighted)
else
w = nrmMandel
endif
forall (i=1_pInt:6_pInt) math_33to6(i) = w(i)*m33(mapNye(1,i),mapNye(2,i))
end function math_33to6
!--------------------------------------------------------------------------------------------------
!> @brief convert 6 vector into symmetric 33 matrix
!> @details Weighted conversion (default) rearranges according to Nye and weights shear
! components according to Mandel. Advisable for matrix operations.
! Unweighted conversion only changes order according to Nye
!--------------------------------------------------------------------------------------------------
pure function math_6to33(v6,weighted)
implicit none
real(pReal), dimension(3,3) :: math_6to33
real(pReal), dimension(6), intent(in) :: v6
logical, optional, intent(in) :: weighted
real(pReal), dimension(6) :: w
integer(pInt) :: i
if(present(weighted)) then
w = merge(invnrmMandel,1.0_pReal,weighted)
else
w = invnrmMandel
endif
do i=1_pInt,6_pInt
math_6to33(mapNye(1,i),mapNye(2,i)) = w(i)*v6(i)
math_6to33(mapNye(2,i),mapNye(1,i)) = w(i)*v6(i)
enddo
end function math_6to33
!--------------------------------------------------------------------------------------------------
!> @brief convert 3333 matrix into vector 99
!--------------------------------------------------------------------------------------------------
pure function math_3333to99(m3333)
implicit none
real(pReal), dimension(9,9) :: math_3333to99
real(pReal), dimension(3,3,3,3), intent(in) :: m3333
real(pReal), dimension(9,9) :: math_Plain3333to99
integer(pInt) :: i,j
forall (i=1_pInt:9_pInt,j=1_pInt:9_pInt) math_Plain3333to99(i,j) = &
m3333(mapPlain(1,i),mapPlain(2,i),mapPlain(1,j),mapPlain(2,j))
do concurrent(i=1_pInt:9_pInt,j=1_pInt:9_pInt)
math_3333to99(i,j) = m3333(mapPlain(1,i),mapPlain(2,i),mapPlain(1,j),mapPlain(2,j))
enddo
end function math_3333to99
end function math_Plain3333to99
!--------------------------------------------------------------------------------------------------
!> @brief plain matrix 99 into 3333 tensor
!> @brief convert 99 vector into 3333 matrix
!--------------------------------------------------------------------------------------------------
pure function math_Plain99to3333(m99)
pure function math_99to3333(m99)
implicit none
real(pReal), dimension(3,3,3,3) :: math_99to3333
real(pReal), dimension(9,9), intent(in) :: m99
real(pReal), dimension(3,3,3,3) :: math_Plain99to3333
integer(pInt) :: i,j
forall (i=1_pInt:9_pInt,j=1_pInt:9_pInt) math_Plain99to3333(mapPlain(1,i),mapPlain(2,i),&
mapPlain(1,j),mapPlain(2,j)) = m99(i,j)
do concurrent(i=1_pInt:9_pInt,j=1_pInt:9_pInt)
math_99to3333(mapPlain(1,i),mapPlain(2,i),mapPlain(1,j),mapPlain(2,j)) = m99(i,j)
enddo
end function math_Plain99to3333
end function math_99to3333
!--------------------------------------------------------------------------------------------------
!> @brief convert Mandel matrix 66 into Plain matrix 66
!> @brief convert symmetric 3333 matrix into 66 vector
!> @details Weighted conversion (default) rearranges according to Nye and weights shear
! components according to Mandel. Advisable for matrix operations.
! Unweighted conversion only changes order according to Nye
!--------------------------------------------------------------------------------------------------
pure function math_Mandel66toPlain66(m66)
implicit none
real(pReal), dimension(6,6), intent(in) :: m66
real(pReal), dimension(6,6) :: math_Mandel66toPlain66
integer(pInt) :: i,j
forall (i=1_pInt:6_pInt,j=1_pInt:6_pInt) &
math_Mandel66toPlain66(i,j) = invnrmMandel(i) * invnrmMandel(j) * m66(i,j)
end function math_Mandel66toPlain66
!--------------------------------------------------------------------------------------------------
!> @brief convert Plain matrix 66 into Mandel matrix 66
!--------------------------------------------------------------------------------------------------
pure function math_Plain66toMandel66(m66)
implicit none
real(pReal), dimension(6,6), intent(in) :: m66
real(pReal), dimension(6,6) :: math_Plain66toMandel66
integer(pInt) :: i,j
forall (i=1_pInt:6_pInt,j=1_pInt:6_pInt) &
math_Plain66toMandel66(i,j) = nrmMandel(i) * nrmMandel(j) * m66(i,j)
end function math_Plain66toMandel66
!--------------------------------------------------------------------------------------------------
!> @brief convert symmetric 3333 tensor into Mandel matrix 66
!--------------------------------------------------------------------------------------------------
pure function math_Mandel3333to66(m3333)
pure function math_3333to66(m3333,weighted)
implicit none
real(pReal), dimension(6,6) :: math_3333to66
real(pReal), dimension(3,3,3,3), intent(in) :: m3333
real(pReal), dimension(6,6) :: math_Mandel3333to66
logical, optional, intent(in) :: weighted
real(pReal), dimension(6) :: w
integer(pInt) :: i,j
if(present(weighted)) then
w = merge(nrmMandel,1.0_pReal,weighted)
else
w = nrmMandel
endif
forall (i=1_pInt:6_pInt,j=1_pInt:6_pInt) math_Mandel3333to66(i,j) = &
nrmMandel(i)*nrmMandel(j)*m3333(mapMandel(1,i),mapMandel(2,i),mapMandel(1,j),mapMandel(2,j))
do concurrent(i=1_pInt:6_pInt,j=1_pInt:6_pInt)
math_3333to66(i,j) = w(i)*w(j)*m3333(mapNye(1,i),mapNye(2,i),mapNye(1,j),mapNye(2,j))
enddo
end function math_Mandel3333to66
end function math_3333to66
!--------------------------------------------------------------------------------------------------
!> @brief convert Mandel matrix 66 back to symmetric 3333 tensor
!> @brief convert 66 vector into symmetric 3333 matrix
!> @details Weighted conversion (default) rearranges according to Nye and weights shear
! components according to Mandel. Advisable for matrix operations.
! Unweighted conversion only changes order according to Nye
!--------------------------------------------------------------------------------------------------
pure function math_Mandel66to3333(m66)
pure function math_66to3333(m66,weighted)
implicit none
real(pReal), dimension(3,3,3,3) :: math_Mandel66to3333
real(pReal), dimension(6,6), intent(in) :: m66
real(pReal), dimension(3,3,3,3) :: math_66to3333
real(pReal), dimension(6,6), intent(in) :: m66
logical, optional, intent(in) :: weighted
real(pReal), dimension(6) :: w
integer(pInt) :: i,j
if(present(weighted)) then
w = merge(invnrmMandel,1.0_pReal,weighted)
else
w = invnrmMandel
endif
forall (i=1_pInt:6_pInt,j=1_pInt:6_pInt)
math_Mandel66to3333(mapMandel(1,i),mapMandel(2,i),mapMandel(1,j),mapMandel(2,j)) = &
invnrmMandel(i)*invnrmMandel(j)*m66(i,j)
math_Mandel66to3333(mapMandel(2,i),mapMandel(1,i),mapMandel(1,j),mapMandel(2,j)) = &
invnrmMandel(i)*invnrmMandel(j)*m66(i,j)
math_Mandel66to3333(mapMandel(1,i),mapMandel(2,i),mapMandel(2,j),mapMandel(1,j)) = &
invnrmMandel(i)*invnrmMandel(j)*m66(i,j)
math_Mandel66to3333(mapMandel(2,i),mapMandel(1,i),mapMandel(2,j),mapMandel(1,j)) = &
invnrmMandel(i)*invnrmMandel(j)*m66(i,j)
end forall
do i=1_pInt,6_pInt; do j=1_pInt, 6_pInt
math_66to3333(mapNye(1,i),mapNye(2,i),mapNye(1,j),mapNye(2,j)) = w(i)*w(j)*m66(i,j)
math_66to3333(mapNye(2,i),mapNye(1,i),mapNye(1,j),mapNye(2,j)) = w(i)*w(j)*m66(i,j)
math_66to3333(mapNye(1,i),mapNye(2,i),mapNye(2,j),mapNye(1,j)) = w(i)*w(j)*m66(i,j)
math_66to3333(mapNye(2,i),mapNye(1,i),mapNye(2,j),mapNye(1,j)) = w(i)*w(j)*m66(i,j)
enddo; enddo
end function math_Mandel66to3333
end function math_66to3333
!--------------------------------------------------------------------------------------------------
@ -1215,16 +1298,12 @@ pure function math_Voigt66to3333(m66)
real(pReal), dimension(6,6), intent(in) :: m66
integer(pInt) :: i,j
forall (i=1_pInt:6_pInt,j=1_pInt:6_pInt)
math_Voigt66to3333(mapVoigt(1,i),mapVoigt(2,i),mapVoigt(1,j),mapVoigt(2,j)) = &
invnrmVoigt(i)*invnrmVoigt(j)*m66(i,j)
math_Voigt66to3333(mapVoigt(2,i),mapVoigt(1,i),mapVoigt(1,j),mapVoigt(2,j)) = &
invnrmVoigt(i)*invnrmVoigt(j)*m66(i,j)
math_Voigt66to3333(mapVoigt(1,i),mapVoigt(2,i),mapVoigt(2,j),mapVoigt(1,j)) = &
invnrmVoigt(i)*invnrmVoigt(j)*m66(i,j)
math_Voigt66to3333(mapVoigt(2,i),mapVoigt(1,i),mapVoigt(2,j),mapVoigt(1,j)) = &
invnrmVoigt(i)*invnrmVoigt(j)*m66(i,j)
end forall
do i=1_pInt,6_pInt; do j=1_pInt, 6_pInt
math_Voigt66to3333(mapVoigt(1,i),mapVoigt(2,i),mapVoigt(1,j),mapVoigt(2,j)) = m66(i,j)
math_Voigt66to3333(mapVoigt(2,i),mapVoigt(1,i),mapVoigt(1,j),mapVoigt(2,j)) = m66(i,j)
math_Voigt66to3333(mapVoigt(1,i),mapVoigt(2,i),mapVoigt(2,j),mapVoigt(1,j)) = m66(i,j)
math_Voigt66to3333(mapVoigt(2,i),mapVoigt(1,i),mapVoigt(2,j),mapVoigt(1,j)) = m66(i,j)
enddo; enddo
end function math_Voigt66to3333
@ -1632,8 +1711,9 @@ pure function math_qToR(q)
real(pReal), dimension(3,3) :: math_qToR, T,S
integer(pInt) :: i, j
forall (i = 1_pInt:3_pInt, j = 1_pInt:3_pInt) &
do concurrent (i = 1_pInt:3_pInt, j = 1_pInt:3_pInt)
T(i,j) = q(i+1_pInt) * q(j+1_pInt)
enddo
S = reshape( [0.0_pReal, -q(4), q(3), &
q(4), 0.0_pReal, -q(2), &
@ -2038,7 +2118,7 @@ end function math_eigenvectorBasisSym
!--------------------------------------------------------------------------------------------------
!> @brief eigenvector basis of symmetric 33 matrix m
!--------------------------------------------------------------------------------------------------
function math_eigenvectorBasisSym33(m)
pure function math_eigenvectorBasisSym33(m)
implicit none
real(pReal), dimension(3,3) :: math_eigenvectorBasisSym33
@ -2103,7 +2183,7 @@ end function math_eigenvectorBasisSym33
!--------------------------------------------------------------------------------------------------
!> @brief logarithm eigenvector basis of symmetric 33 matrix m
!--------------------------------------------------------------------------------------------------
function math_eigenvectorBasisSym33_log(m)
pure function math_eigenvectorBasisSym33_log(m)
implicit none
real(pReal), dimension(3,3) :: math_eigenvectorBasisSym33_log
@ -2159,11 +2239,12 @@ function math_eigenvectorBasisSym33_log(m)
endif threeSimilarEigenvalues
math_eigenvectorBasisSym33_log = log(sqrt(values(1))) * EB(1:3,1:3,1) &
+ log(sqrt(values(2))) * EB(1:3,1:3,2) &
+ log(sqrt(values(3))) * EB(1:3,1:3,3)
+ log(sqrt(values(2))) * EB(1:3,1:3,2) &
+ log(sqrt(values(3))) * EB(1:3,1:3,3)
end function math_eigenvectorBasisSym33_log
!--------------------------------------------------------------------------------------------------
!> @brief rotational part from polar decomposition of 33 tensor m
!--------------------------------------------------------------------------------------------------
@ -2608,13 +2689,12 @@ pure function math_rotate_forward3333(tensor,rot_tensor)
real(pReal), dimension(3,3,3,3), intent(in) :: tensor
integer(pInt) :: i,j,k,l,m,n,o,p
math_rotate_forward3333= 0.0_pReal
do i = 1_pInt,3_pInt; do j = 1_pInt,3_pInt; do k = 1_pInt,3_pInt; do l = 1_pInt,3_pInt
do m = 1_pInt,3_pInt; do n = 1_pInt,3_pInt; do o = 1_pInt,3_pInt; do p = 1_pInt,3_pInt
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 = 0.0_pReal
do i = 1_pInt,3_pInt;do j = 1_pInt,3_pInt;do k = 1_pInt,3_pInt;do l = 1_pInt,3_pInt
do m = 1_pInt,3_pInt;do n = 1_pInt,3_pInt;do o = 1_pInt,3_pInt;do p = 1_pInt,3_pInt
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)
enddo; enddo; enddo; enddo; enddo; enddo; enddo; enddo
end function math_rotate_forward3333