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

Merge branch 'tensor-conversion-rename' into development

This commit is contained in:
Martin Diehl 2019-01-15 11:43:45 +01:00
parent 42c86b9711 3c12b7441f
commit c231c808da
4 changed files with 286 additions and 230 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

42
src/CPFEM.f90
View File

@ -286,12 +286,11 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
math_identity2nd, & math_identity2nd, &
math_mul33x33, & math_mul33x33, &
math_det33, & math_det33, &
math_transpose33, & math_delta, &
math_I3, & math_sym3333to66, &
math_Mandel3333to66, & math_66toSym3333, &
math_Mandel66to3333, & math_sym33to6, &
math_Mandel33to6, & math_6toSym33
math_Mandel6to33
use mesh, only: & use mesh, only: &
mesh_FEasCP, & mesh_FEasCP, &
mesh_NcpElems, & mesh_NcpElems, &
@ -353,8 +352,8 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
integer(pInt), intent(in) :: mode !< computation mode 1: regular computation plus aging of results integer(pInt), intent(in) :: mode !< computation mode 1: regular computation plus aging of results
real(pReal), intent(in) :: temperature_inp !< temperature real(pReal), intent(in) :: temperature_inp !< temperature
logical, intent(in) :: parallelExecution !< flag indicating parallel computation of requested IPs logical, intent(in) :: parallelExecution !< flag indicating parallel computation of requested IPs
real(pReal), dimension(6), intent(out) :: cauchyStress !< stress vector in Mandel notation real(pReal), dimension(6), intent(out) :: cauchyStress !< stress as 6 vector
real(pReal), dimension(6,6), intent(out) :: jacobian !< jacobian in Mandel notation (Consistent tangent dcs/dE) real(pReal), dimension(6,6), intent(out) :: jacobian !< jacobian as 66 tensor (Consistent tangent dcs/dE)
real(pReal) J_inverse, & ! inverse of Jacobian real(pReal) J_inverse, & ! inverse of Jacobian
rnd rnd
@ -534,8 +533,8 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) then if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) then
write(6,'(a,1x,i8,1x,i2)') '<< CPFEM >> OUTDATED at elFE ip',elFE,ip write(6,'(a,1x,i8,1x,i2)') '<< CPFEM >> OUTDATED at elFE ip',elFE,ip
write(6,'(a,/,3(12x,3(f10.6,1x),/))') '<< CPFEM >> FFN1 old:',& write(6,'(a,/,3(12x,3(f10.6,1x),/))') '<< CPFEM >> FFN1 old:',&
math_transpose33(materialpoint_F(1:3,1:3,ip,elCP)) transpose(materialpoint_F(1:3,1:3,ip,elCP))
write(6,'(a,/,3(12x,3(f10.6,1x),/))') '<< CPFEM >> FFN1 now:',math_transpose33(ffn1) write(6,'(a,/,3(12x,3(f10.6,1x),/))') '<< CPFEM >> FFN1 now:',transpose(ffn1)
endif endif
outdatedFFN1 = .true. outdatedFFN1 = .true.
endif endif
@ -593,26 +592,25 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
endif endif
! translate from P to CS ! translate from P to CS
Kirchhoff = math_mul33x33(materialpoint_P(1:3,1:3,ip,elCP), math_transpose33(materialpoint_F(1:3,1:3,ip,elCP))) Kirchhoff = math_mul33x33(materialpoint_P(1:3,1:3,ip,elCP), transpose(materialpoint_F(1:3,1:3,ip,elCP)))
J_inverse = 1.0_pReal / math_det33(materialpoint_F(1:3,1:3,ip,elCP)) J_inverse = 1.0_pReal / math_det33(materialpoint_F(1:3,1:3,ip,elCP))
CPFEM_cs(1:6,ip,elCP) = math_Mandel33to6(J_inverse * Kirchhoff) CPFEM_cs(1:6,ip,elCP) = math_sym33to6(J_inverse * Kirchhoff,weighted=.false.)
! translate from dP/dF to dCS/dE ! translate from dP/dF to dCS/dE
H = 0.0_pReal H = 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 i=1,3; do j=1,3; do k=1,3; do l=1,3; do m=1,3; do n=1,3
H(i,j,k,l) = H(i,j,k,l) + & H(i,j,k,l) = H(i,j,k,l) &
materialpoint_F(j,m,ip,elCP) * & + materialpoint_F(j,m,ip,elCP) * materialpoint_F(l,n,ip,elCP) &
materialpoint_F(l,n,ip,elCP) * & * materialpoint_dPdF(i,m,k,n,ip,elCP) &
materialpoint_dPdF(i,m,k,n,ip,elCP) - & - math_delta(j,l) * materialpoint_F(i,m,ip,elCP) * materialpoint_P(k,m,ip,elCP) &
math_I3(j,l) * materialpoint_F(i,m,ip,elCP) * materialpoint_P(k,m,ip,elCP) + & + 0.5_pReal * ( Kirchhoff(j,l)*math_delta(i,k) + Kirchhoff(i,k)*math_delta(j,l) &
0.5_pReal * (math_I3(i,k) * Kirchhoff(j,l) + math_I3(j,l) * Kirchhoff(i,k) + & + Kirchhoff(j,k)*math_delta(i,l) + Kirchhoff(i,l)*math_delta(j,k))
math_I3(i,l) * Kirchhoff(j,k) + math_I3(j,k) * Kirchhoff(i,l))
enddo; enddo; enddo; enddo; enddo; enddo enddo; enddo; enddo; enddo; enddo; enddo
forall(i=1:3, j=1:3,k=1:3,l=1:3) & forall(i=1:3, j=1:3,k=1:3,l=1:3) &
H_sym(i,j,k,l) = 0.25_pReal * (H(i,j,k,l) + H(j,i,k,l) + H(i,j,l,k) + H(j,i,l,k)) H_sym(i,j,k,l) = 0.25_pReal * (H(i,j,k,l) + H(j,i,k,l) + H(i,j,l,k) + H(j,i,l,k))
CPFEM_dcsde(1:6,1:6,ip,elCP) = math_Mandel3333to66(J_inverse * H_sym) CPFEM_dcsde(1:6,1:6,ip,elCP) = math_sym3333to66(J_inverse * H_sym,weighted=.false.)
endif terminalIllness endif terminalIllness
endif validCalculation endif validCalculation
@ -639,7 +637,7 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
!*** remember extreme values of stress ... !*** remember extreme values of stress ...
cauchyStress33 = math_Mandel6to33(CPFEM_cs(1:6,ip,elCP)) cauchyStress33 = math_6toSym33(CPFEM_cs(1:6,ip,elCP),weighted=.false.)
if (maxval(cauchyStress33) > debug_stressMax) then if (maxval(cauchyStress33) > debug_stressMax) then
debug_stressMaxLocation = [elCP, ip] debug_stressMaxLocation = [elCP, ip]
debug_stressMax = maxval(cauchyStress33) debug_stressMax = maxval(cauchyStress33)
@ -649,7 +647,7 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
debug_stressMin = minval(cauchyStress33) debug_stressMin = minval(cauchyStress33)
endif endif
!*** ... and Jacobian !*** ... and Jacobian
jacobian3333 = math_Mandel66to3333(CPFEM_dcsdE(1:6,1:6,ip,elCP)) jacobian3333 = math_66toSym3333(CPFEM_dcsdE(1:6,1:6,ip,elCP),weighted=.false.)
if (maxval(jacobian3333) > debug_jacobianMax) then if (maxval(jacobian3333) > debug_jacobianMax) then
debug_jacobianMaxLocation = [elCP, ip] debug_jacobianMaxLocation = [elCP, ip]
debug_jacobianMax = maxval(jacobian3333) debug_jacobianMax = maxval(jacobian3333)

8
src/DAMASK_abaqus.f
View File

@ -102,8 +102,6 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,&
calcMode, & calcMode, &
terminallyIll, & terminallyIll, &
symmetricSolver symmetricSolver
use math, only: &
invnrmMandel
use debug, only: & use debug, only: &
debug_info, & debug_info, &
debug_reset, & debug_reset, &
@ -305,9 +303,9 @@ subroutine UMAT(STRESS,STATEV,DDSDDE,SSE,SPD,SCD,&
! ABAQUS implicit: 11, 22, 33, 12, 13, 23 ! ABAQUS implicit: 11, 22, 33, 12, 13, 23
! ABAQUS implicit: 11, 22, 33, 12 ! ABAQUS implicit: 11, 22, 33, 12
forall(i=1:ntens) ddsdde(1:ntens,i) = invnrmMandel(i)*ddsdde_h(1:ntens,i)*invnrmMandel(1:ntens) ddsdde = ddsdde_h(1:ntens,1:ntens)
stress(1:ntens) = stress_h(1:ntens)*invnrmMandel(1:ntens) stress = stress_h(1:ntens)
if(symmetricSolver) ddsdde(1:ntens,1:ntens) = 0.5_pReal*(ddsdde(1:ntens,1:ntens) + transpose(ddsdde(1:ntens,1:ntens))) if(symmetricSolver) ddsdde = 0.5_pReal*(ddsdde + transpose(ddsdde))
if(ntens == 6) then if(ntens == 6) then
stress_h = stress stress_h = stress
stress(5) = stress_h(6) stress(5) = stress_h(6)

11
src/DAMASK_marc.f90
View File

@ -127,9 +127,6 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
calcMode, & calcMode, &
terminallyIll, & terminallyIll, &
symmetricSolver symmetricSolver
use math, only: &
math_transpose33,&
invnrmMandel
use debug, only: & use debug, only: &
debug_level, & debug_level, &
debug_LEVELBASIC, & debug_LEVELBASIC, &
@ -235,9 +232,9 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
write(6,'(a,i12)') ' Nodes: ', nnode write(6,'(a,i12)') ' Nodes: ', nnode
write(6,'(a,i1)') ' Deformation gradient: ', itel write(6,'(a,i1)') ' Deformation gradient: ', itel
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' Deformation gradient at t=n:', & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' Deformation gradient at t=n:', &
math_transpose33(ffn) transpose(ffn)
write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' Deformation gradient at t=n+1:', & write(6,'(/,a,/,3(3(f12.7,1x)/))',advance='no') ' Deformation gradient at t=n+1:', &
math_transpose33(ffn1) transpose(ffn1)
endif endif
!$ defaultNumThreadsInt = omp_get_num_threads() ! remember number of threads set by Marc !$ defaultNumThreadsInt = omp_get_num_threads() ! remember number of threads set by Marc
@ -357,8 +354,8 @@ subroutine hypela2(d,g,e,de,s,t,dt,ngens,m,nn,kcus,matus,ndi,nshear,disp, &
! Marc: 11, 22, 33, 12, 23, 13 ! Marc: 11, 22, 33, 12, 23, 13
! Marc: 11, 22, 33, 12 ! Marc: 11, 22, 33, 12
forall(i=1:ngens) d(1:ngens,i) = invnrmMandel(i)*ddsdde(1:ngens,i)*invnrmMandel(1:ngens) d = ddsdde(1:ngens,1:ngens)
s(1:ndi+nshear) = stress(1:ndi+nshear)*invnrmMandel(1:ndi+nshear) s = stress(1:ndi+nshear)
g = 0.0_pReal g = 0.0_pReal
if(symmetricSolver) d = 0.5_pReal*(d+transpose(d)) if(symmetricSolver) d = 0.5_pReal*(d+transpose(d))

405
src/math.f90
View File

@ -24,33 +24,25 @@ module math
0.0_pReal,0.0_pReal,1.0_pReal & 0.0_pReal,0.0_pReal,1.0_pReal &
],[3,3]) !< 3x3 Identity ],[3,3]) !< 3x3 Identity
! ToDo MD: Our naming scheme is a little bit odd: We use essentially the re-ordering according to Nye real(pReal), dimension(6), parameter, private :: &
! (convenient because Abaqus and Marc want to have 12 on position 4) nrmMandel = [&
! but weight the shear components according to Mandel (convenient for matrix multiplications) 1.0_pReal, 1.0_pReal, 1.0_pReal, &
! I suggest to keep Voigt3333to66 (required for reading in elasticity matrices) but rename sqrt(2.0_pReal), sqrt(2.0_pReal), sqrt(2.0_pReal) ] !< weighting for Mandel notation (forward)
! 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" real(pReal), dimension(6), parameter , private :: &
! without leaving the impression that it follows any established convention invnrmMandel = [&
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 :: & integer(pInt), dimension (2,6), parameter, private :: &
mapMandel = reshape([& mapNye = reshape([&
1_pInt,1_pInt, & 1_pInt,1_pInt, &
2_pInt,2_pInt, & 2_pInt,2_pInt, &
3_pInt,3_pInt, & 3_pInt,3_pInt, &
1_pInt,2_pInt, & 1_pInt,2_pInt, &
2_pInt,3_pInt, & 2_pInt,3_pInt, &
1_pInt,3_pInt & 1_pInt,3_pInt &
],[2,6]) !< arrangement in Mandel notation. Differs from https://en.wikipedia.org/wiki/Voigt_notation#Mandel_notation ],[2,6]) !< arrangement in Nye notation.
real(pReal), dimension(6), parameter, private :: &
nrmMandel = [&
1.0_pReal, 1.0_pReal, 1.0_pReal, &
sqrt(2.0_pReal), sqrt(2.0_pReal), sqrt(2.0_pReal) ] !< weighting for Mandel notation (forward)
real(pReal), dimension(6), parameter , public :: &
invnrmMandel = [&
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 :: & integer(pInt), dimension (2,6), parameter, private :: &
mapVoigt = reshape([& mapVoigt = reshape([&
@ -62,10 +54,6 @@ module math
1_pInt,2_pInt & 1_pInt,2_pInt &
],[2,6]) !< arrangement in Voigt notation ],[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 :: & integer(pInt), dimension (2,9), parameter, private :: &
mapPlain = reshape([& mapPlain = reshape([&
1_pInt,1_pInt, & 1_pInt,1_pInt, &
@ -79,6 +67,56 @@ module math
3_pInt,3_pInt & 3_pInt,3_pInt &
],[2,9]) !< arrangement in Plain notation ],[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_sym33to6
end interface math_Mandel33to6
interface math_Mandel6to33
module procedure math_6toSym33
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_sym3333to66
end interface math_Mandel3333to66
interface math_Mandel66to3333
module procedure math_66toSym3333
end interface math_Mandel66to3333
public :: &
math_Plain33to9, &
math_Plain9to33, &
math_Mandel33to6, &
math_Mandel6to33, &
math_Plain3333to99, &
math_Plain99to3333, &
math_Mandel3333to66, &
math_Mandel66to3333
!---------------------------------------------------------------------------------------------------
public :: & public :: &
math_init, & math_init, &
math_qsort, & math_qsort, &
@ -116,16 +154,14 @@ module math
math_equivStress33, & math_equivStress33, &
math_trace33, & math_trace33, &
math_det33, & math_det33, &
math_Plain33to9, & math_33to9, &
math_Plain9to33, & math_9to33, &
math_Mandel33to6, & math_sym33to6, &
math_Mandel6to33, & math_6toSym33, &
math_Plain3333to99, & math_3333to99, &
math_Plain99to3333, & math_99to3333, &
math_Mandel66toPlain66, & math_sym3333to66, &
math_Plain66toMandel66, & math_66toSym3333, &
math_Mandel3333to66, &
math_Mandel66to3333, &
math_Voigt66to3333, & math_Voigt66to3333, &
math_qRand, & math_qRand, &
math_qMul, & math_qMul, &
@ -423,7 +459,7 @@ pure function math_identity2nd(dimen)
real(pReal), dimension(dimen,dimen) :: math_identity2nd real(pReal), dimension(dimen,dimen) :: math_identity2nd
math_identity2nd = 0.0_pReal math_identity2nd = 0.0_pReal
forall (i=1_pInt:dimen) math_identity2nd(i,i) = 1.0_pReal forall(i=1_pInt:dimen) math_identity2nd(i,i) = 1.0_pReal
end function math_identity2nd end function math_identity2nd
@ -437,9 +473,11 @@ pure function math_identity4th(dimen)
integer(pInt), intent(in) :: dimen !< tensor dimension integer(pInt), intent(in) :: dimen !< tensor dimension
integer(pInt) :: i,j,k,l integer(pInt) :: i,j,k,l
real(pReal), dimension(dimen,dimen,dimen,dimen) :: math_identity4th 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) = & identity2nd = math_identity2nd(dimen)
0.5_pReal*(math_I3(i,k)*math_I3(j,l)+math_I3(i,l)*math_I3(j,k)) 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*(identity2nd(i,k)*identity2nd(j,l)+identity2nd(i,l)*identity2nd(j,k))
end function math_identity4th end function math_identity4th
@ -508,7 +546,7 @@ pure function math_tensorproduct(A,B)
real(pReal), dimension(size(A,1),size(B,1)) :: math_tensorproduct real(pReal), dimension(size(A,1),size(B,1)) :: math_tensorproduct
integer(pInt) :: i,j integer(pInt) :: i,j
forall (i=1_pInt:size(A,1),j=1_pInt:size(B,1)) math_tensorproduct(i,j) = A(i)*B(j) forall(i=1_pInt:size(A,1),j=1_pInt:size(B,1)) math_tensorproduct(i,j) = A(i)*B(j)
end function math_tensorproduct end function math_tensorproduct
@ -523,7 +561,7 @@ pure function math_tensorproduct33(A,B)
real(pReal), dimension(3), intent(in) :: A,B real(pReal), dimension(3), intent(in) :: A,B
integer(pInt) :: i,j integer(pInt) :: i,j
forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt) math_tensorproduct33(i,j) = A(i)*B(j) forall(i=1_pInt:3_pInt,j=1_pInt:3_pInt) math_tensorproduct33(i,j) = A(i)*B(j)
end function math_tensorproduct33 end function math_tensorproduct33
@ -564,7 +602,7 @@ real(pReal) pure function math_mul33xx33(A,B)
integer(pInt) :: i,j integer(pInt) :: i,j
real(pReal), dimension(3,3) :: C 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) forall(i=1_pInt:3_pInt,j=1_pInt:3_pInt) C(i,j) = A(i,j) * B(i,j)
math_mul33xx33 = sum(C) math_mul33xx33 = sum(C)
end function math_mul33xx33 end function math_mul33xx33
@ -581,8 +619,7 @@ pure function math_mul3333xx33(A,B)
real(pReal), dimension(3,3), intent(in) :: B real(pReal), dimension(3,3), intent(in) :: B
integer(pInt) :: i,j integer(pInt) :: i,j
forall(i = 1_pInt:3_pInt,j = 1_pInt:3_pInt) & forall(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))
math_mul3333xx33(i,j) = sum(A(i,j,1:3,1:3)*B(1:3,1:3))
end function math_mul3333xx33 end function math_mul3333xx33
@ -614,8 +651,7 @@ pure function math_mul33x33(A,B)
real(pReal), dimension(3,3), intent(in) :: A,B real(pReal), dimension(3,3), intent(in) :: A,B
integer(pInt) :: i,j integer(pInt) :: i,j
forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt) & forall(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)
math_mul33x33(i,j) = A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j)
end function math_mul33x33 end function math_mul33x33
@ -630,9 +666,9 @@ pure function math_mul66x66(A,B)
real(pReal), dimension(6,6), intent(in) :: A,B real(pReal), dimension(6,6), intent(in) :: A,B
integer(pInt) :: i,j integer(pInt) :: i,j
forall (i=1_pInt:6_pInt,j=1_pInt:6_pInt) math_mul66x66(i,j) = & forall(i=1_pInt:6_pInt,j=1_pInt:6_pInt) &
A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j) + & 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) + A(i,4)*B(4,j) + A(i,5)*B(5,j) + A(i,6)*B(6,j)
end function math_mul66x66 end function math_mul66x66
@ -647,10 +683,10 @@ pure function math_mul99x99(A,B)
real(pReal), dimension(9,9), intent(in) :: A,B real(pReal), dimension(9,9), intent(in) :: A,B
integer(pInt) i,j integer(pInt) i,j
forall (i=1_pInt:9_pInt,j=1_pInt:9_pInt) math_mul99x99(i,j) = & forall(i=1_pInt:9_pInt,j=1_pInt:9_pInt) &
A(i,1)*B(1,j) + A(i,2)*B(2,j) + A(i,3)*B(3,j) + & 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,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) + A(i,7)*B(7,j) + A(i,8)*B(8,j) + A(i,9)*B(9,j)
end function math_mul99x99 end function math_mul99x99
@ -698,9 +734,8 @@ pure function math_mul66x6(A,B)
real(pReal), dimension(6), intent(in) :: B real(pReal), dimension(6), intent(in) :: B
integer(pInt) :: i integer(pInt) :: i
forall (i=1_pInt:6_pInt) math_mul66x6(i) = & forall (i=1_pInt:6_pInt) math_mul66x6(i) = A(i,1)*B(1) + A(i,2)*B(2) + A(i,3)*B(3) &
A(i,1)*B(1) + A(i,2)*B(2) + A(i,3)*B(3) + & + A(i,4)*B(4) + A(i,5)*B(5) + A(i,6)*B(6)
A(i,4)*B(4) + A(i,5)*B(5) + A(i,6)*B(6)
end function math_mul66x6 end function math_mul66x6
@ -747,8 +782,8 @@ end function math_transpose33
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Cramer inversion of 33 matrix (function) !> @brief Cramer inversion of 33 matrix (function)
! direct Cramer inversion of matrix A. !> @details Direct Cramer inversion of matrix A. Returns all zeroes if not possible, i.e.
! returns all zeroes if not possible, i.e. if det close to zero ! if determinant is close to zero
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function math_inv33(A) pure function math_inv33(A)
use prec, only: & use prec, only: &
@ -784,9 +819,9 @@ end function math_inv33
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Cramer inversion of 33 matrix (subroutine) !> @brief Cramer inversion of 33 matrix (subroutine)
! direct Cramer inversion of matrix A. !> @details Direct Cramer inversion of matrix A. Also returns determinant
! also returns determinant ! Returns an error if not possible, i.e. if determinant is close to zero
! returns error if not possible, i.e. if det close to zero ! ToDo: Output arguments should be first
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure subroutine math_invert33(A, InvA, DetA, error) pure subroutine math_invert33(A, InvA, DetA, error)
use prec, only: & use prec, only: &
@ -843,20 +878,38 @@ function math_invSym3333(A)
dgetrf, & dgetrf, &
dgetri dgetri
temp66_real = math_Mandel3333to66(A) temp66_real = math_sym3333to66(A)
call dgetrf(6,6,temp66_real,6,ipiv6,ierr) call dgetrf(6,6,temp66_real,6,ipiv6,ierr)
call dgetri(6,temp66_real,6,ipiv6,work6,6,ierr) call dgetri(6,temp66_real,6,ipiv6,work6,6,ierr)
if (ierr == 0_pInt) then if (ierr == 0_pInt) then
math_invSym3333 = math_Mandel66to3333(temp66_real) math_invSym3333 = math_66toSym3333(temp66_real)
else else
call IO_error(400_pInt, ext_msg = 'math_invSym3333') call IO_error(400_pInt, ext_msg = 'math_invSym3333')
endif endif
end function math_invSym3333 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 !> @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) subroutine math_invert(myDim,A, InvA, error)
@ -961,15 +1014,14 @@ pure function math_equivStrain33(m)
real(pReal), dimension(3,3), intent(in) :: m real(pReal), dimension(3,3), intent(in) :: m
real(pReal), dimension(3) :: e,s real(pReal), dimension(3) :: e,s
real(pReal) :: math_equivStrain33 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), & 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(2,2)-m(3,3)-m(1,1), &
2.0_pReal*m(3,3)-m(1,1)-m(2,2)]/3.0_pReal 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 s = [m(1,2),m(2,3),m(1,3)]*2.0_pReal
math_equivStrain33 = TWOTHIRD*(1.50_pReal*(sum(e**2.0_pReal)) + & math_equivStrain33 = 2.0_pReal/3.0_pReal &
0.75_pReal*(sum(s**2.0_pReal)))**(0.5_pReal) * (1.50_pReal*(sum(e**2.0_pReal))+ 0.75_pReal*(sum(s**2.0_pReal)))**(0.5_pReal)
end function math_equivStrain33 end function math_equivStrain33
@ -1041,172 +1093,188 @@ end function math_detSym33
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief convert 33 matrix into vector 9 !> @brief convert 33 matrix into vector 9
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function math_Plain33to9(m33) pure function math_33to9(m33)
implicit none implicit none
real(pReal), dimension(9) :: math_Plain33to9 real(pReal), dimension(9) :: math_33to9
real(pReal), dimension(3,3), intent(in) :: m33 real(pReal), dimension(3,3), intent(in) :: m33
integer(pInt) :: i integer(pInt) :: i
forall (i=1_pInt:9_pInt) math_Plain33to9(i) = m33(mapPlain(1,i),mapPlain(2,i)) forall(i=1_pInt:9_pInt) math_33to9(i) = m33(mapPlain(1,i),mapPlain(2,i))
end function math_Plain33to9 end function math_33to9
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief convert Plain 9 back to 33 matrix !> @brief convert 9 vector into 33 matrix
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function math_Plain9to33(v9) pure function math_9to33(v9)
implicit none implicit none
real(pReal), dimension(3,3) :: math_Plain9to33 real(pReal), dimension(3,3) :: math_9to33
real(pReal), dimension(9), intent(in) :: v9 real(pReal), dimension(9), intent(in) :: v9
integer(pInt) :: i integer(pInt) :: i
forall (i=1_pInt:9_pInt) math_Plain9to33(mapPlain(1,i),mapPlain(2,i)) = v9(i) forall(i=1_pInt:9_pInt) math_9to33(mapPlain(1,i),mapPlain(2,i)) = v9(i)
end function math_Plain9to33 end function math_9to33
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief convert symmetric 33 matrix into Mandel vector 6 !> @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_Mandel33to6(m33) pure function math_sym33to6(m33,weighted)
implicit none implicit none
real(pReal), dimension(6) :: math_Mandel33to6 real(pReal), dimension(6) :: math_sym33to6
real(pReal), dimension(3,3), intent(in) :: m33 real(pReal), dimension(3,3), intent(in) :: m33
logical, optional, intent(in) :: weighted
real(pReal), dimension(6) :: w
integer(pInt) :: i integer(pInt) :: i
forall (i=1_pInt:6_pInt) math_Mandel33to6(i) = nrmMandel(i)*m33(mapMandel(1,i),mapMandel(2,i)) if(present(weighted)) then
w = merge(nrmMandel,1.0_pReal,weighted)
else
w = nrmMandel
endif
end function math_Mandel33to6 forall(i=1_pInt:6_pInt) math_sym33to6(i) = w(i)*m33(mapNye(1,i),mapNye(2,i))
end function math_sym33to6
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief convert Mandel 6 back to symmetric 33 matrix !> @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_Mandel6to33(v6) pure function math_6toSym33(v6,weighted)
implicit none implicit none
real(pReal), dimension(3,3) :: math_6toSym33
real(pReal), dimension(6), intent(in) :: v6 real(pReal), dimension(6), intent(in) :: v6
real(pReal), dimension(3,3) :: math_Mandel6to33 logical, optional, intent(in) :: weighted
real(pReal), dimension(6) :: w
integer(pInt) :: i integer(pInt) :: i
forall (i=1_pInt:6_pInt) if(present(weighted)) then
math_Mandel6to33(mapMandel(1,i),mapMandel(2,i)) = invnrmMandel(i)*v6(i) w = merge(invnrmMandel,1.0_pReal,weighted)
math_Mandel6to33(mapMandel(2,i),mapMandel(1,i)) = invnrmMandel(i)*v6(i) else
end forall w = invnrmMandel
endif
end function math_Mandel6to33 do i=1_pInt,6_pInt
math_6toSym33(mapNye(1,i),mapNye(2,i)) = w(i)*v6(i)
math_6toSym33(mapNye(2,i),mapNye(1,i)) = w(i)*v6(i)
enddo
end function math_6toSym33
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief convert 3333 tensor into plain matrix 99 !> @brief convert 3333 matrix into 99 matrix
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function math_Plain3333to99(m3333) pure function math_3333to99(m3333)
implicit none implicit none
real(pReal), dimension(9,9) :: math_3333to99
real(pReal), dimension(3,3,3,3), intent(in) :: m3333 real(pReal), dimension(3,3,3,3), intent(in) :: m3333
real(pReal), dimension(9,9) :: math_Plain3333to99
integer(pInt) :: i,j integer(pInt) :: i,j
forall (i=1_pInt:9_pInt,j=1_pInt:9_pInt) math_Plain3333to99(i,j) = & forall(i=1_pInt:9_pInt,j=1_pInt:9_pInt) &
m3333(mapPlain(1,i),mapPlain(2,i),mapPlain(1,j),mapPlain(2,j)) math_3333to99(i,j) = m3333(mapPlain(1,i),mapPlain(2,i),mapPlain(1,j),mapPlain(2,j))
end function math_3333to99
end function math_Plain3333to99
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief plain matrix 99 into 3333 tensor !> @brief convert 99 matrix into 3333 matrix
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function math_Plain99to3333(m99) pure function math_99to3333(m99)
implicit none implicit none
real(pReal), dimension(3,3,3,3) :: math_99to3333
real(pReal), dimension(9,9), intent(in) :: m99 real(pReal), dimension(9,9), intent(in) :: m99
real(pReal), dimension(3,3,3,3) :: math_Plain99to3333
integer(pInt) :: i,j integer(pInt) :: i,j
forall (i=1_pInt:9_pInt,j=1_pInt:9_pInt) math_Plain99to3333(mapPlain(1,i),mapPlain(2,i),& forall(i=1_pInt:9_pInt,j=1_pInt:9_pInt) &
mapPlain(1,j),mapPlain(2,j)) = m99(i,j) math_99to3333(mapPlain(1,i),mapPlain(2,i),mapPlain(1,j),mapPlain(2,j)) = m99(i,j)
end function math_Plain99to3333 end function math_99to3333
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief convert Mandel matrix 66 into Plain matrix 66 !> @brief convert symmetric 3333 matrix into 66 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_Mandel66toPlain66(m66) pure function math_sym3333to66(m3333,weighted)
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)
implicit none implicit none
real(pReal), dimension(6,6) :: math_sym3333to66
real(pReal), dimension(3,3,3,3), intent(in) :: m3333 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 integer(pInt) :: i,j
forall (i=1_pInt:6_pInt,j=1_pInt:6_pInt) math_Mandel3333to66(i,j) = & if(present(weighted)) then
nrmMandel(i)*nrmMandel(j)*m3333(mapMandel(1,i),mapMandel(2,i),mapMandel(1,j),mapMandel(2,j)) w = merge(nrmMandel,1.0_pReal,weighted)
else
w = nrmMandel
endif
end function math_Mandel3333to66 forall(i=1_pInt:6_pInt,j=1_pInt:6_pInt) &
math_sym3333to66(i,j) = w(i)*w(j)*m3333(mapNye(1,i),mapNye(2,i),mapNye(1,j),mapNye(2,j))
end function math_sym3333to66
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief convert Mandel matrix 66 back to symmetric 3333 tensor !> @brief convert 66 matrix 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_66toSym3333(m66,weighted)
implicit none implicit none
real(pReal), dimension(3,3,3,3) :: math_Mandel66to3333 real(pReal), dimension(3,3,3,3) :: math_66toSym3333
real(pReal), dimension(6,6), intent(in) :: m66 real(pReal), dimension(6,6), intent(in) :: m66
logical, optional, intent(in) :: weighted
real(pReal), dimension(6) :: w
integer(pInt) :: i,j integer(pInt) :: i,j
forall (i=1_pInt:6_pInt,j=1_pInt:6_pInt) if(present(weighted)) then
math_Mandel66to3333(mapMandel(1,i),mapMandel(2,i),mapMandel(1,j),mapMandel(2,j)) = & w = merge(invnrmMandel,1.0_pReal,weighted)
invnrmMandel(i)*invnrmMandel(j)*m66(i,j) else
math_Mandel66to3333(mapMandel(2,i),mapMandel(1,i),mapMandel(1,j),mapMandel(2,j)) = & w = invnrmMandel
invnrmMandel(i)*invnrmMandel(j)*m66(i,j) endif
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
end function math_Mandel66to3333 do i=1_pInt,6_pInt; do j=1_pInt, 6_pInt
math_66toSym3333(mapNye(1,i),mapNye(2,i),mapNye(1,j),mapNye(2,j)) = w(i)*w(j)*m66(i,j)
math_66toSym3333(mapNye(2,i),mapNye(1,i),mapNye(1,j),mapNye(2,j)) = w(i)*w(j)*m66(i,j)
math_66toSym3333(mapNye(1,i),mapNye(2,i),mapNye(2,j),mapNye(1,j)) = w(i)*w(j)*m66(i,j)
math_66toSym3333(mapNye(2,i),mapNye(1,i),mapNye(2,j),mapNye(1,j)) = w(i)*w(j)*m66(i,j)
enddo; enddo
end function math_66toSym3333
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief convert Voigt matrix 66 back to symmetric 3333 tensor !> @brief convert 66 Voigt matrix into symmetric 3333 matrix
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function math_Voigt66to3333(m66) pure function math_Voigt66to3333(m66)
@ -1215,16 +1283,12 @@ pure function math_Voigt66to3333(m66)
real(pReal), dimension(6,6), intent(in) :: m66 real(pReal), dimension(6,6), intent(in) :: m66
integer(pInt) :: i,j integer(pInt) :: i,j
forall (i=1_pInt:6_pInt,j=1_pInt:6_pInt) 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)) = & math_Voigt66to3333(mapVoigt(1,i),mapVoigt(2,i),mapVoigt(1,j),mapVoigt(2,j)) = m66(i,j)
invnrmVoigt(i)*invnrmVoigt(j)*m66(i,j) math_Voigt66to3333(mapVoigt(2,i),mapVoigt(1,i),mapVoigt(1,j),mapVoigt(2,j)) = m66(i,j)
math_Voigt66to3333(mapVoigt(2,i),mapVoigt(1,i),mapVoigt(1,j),mapVoigt(2,j)) = & math_Voigt66to3333(mapVoigt(1,i),mapVoigt(2,i),mapVoigt(2,j),mapVoigt(1,j)) = m66(i,j)
invnrmVoigt(i)*invnrmVoigt(j)*m66(i,j) math_Voigt66to3333(mapVoigt(2,i),mapVoigt(1,i),mapVoigt(2,j),mapVoigt(1,j)) = m66(i,j)
math_Voigt66to3333(mapVoigt(1,i),mapVoigt(2,i),mapVoigt(2,j),mapVoigt(1,j)) = & enddo; enddo
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
end function math_Voigt66to3333 end function math_Voigt66to3333
@ -1632,8 +1696,7 @@ pure function math_qToR(q)
real(pReal), dimension(3,3) :: math_qToR, T,S real(pReal), dimension(3,3) :: math_qToR, T,S
integer(pInt) :: i, j integer(pInt) :: i, j
forall (i = 1_pInt:3_pInt, j = 1_pInt:3_pInt) & forall(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt) T(i,j) = q(i+1_pInt) * q(j+1_pInt)
T(i,j) = q(i+1_pInt) * q(j+1_pInt)
S = reshape( [0.0_pReal, -q(4), q(3), & S = reshape( [0.0_pReal, -q(4), q(3), &
q(4), 0.0_pReal, -q(2), & q(4), 0.0_pReal, -q(2), &
@ -2038,7 +2101,7 @@ end function math_eigenvectorBasisSym
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief eigenvector basis of symmetric 33 matrix m !> @brief eigenvector basis of symmetric 33 matrix m
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function math_eigenvectorBasisSym33(m) pure function math_eigenvectorBasisSym33(m)
implicit none implicit none
real(pReal), dimension(3,3) :: math_eigenvectorBasisSym33 real(pReal), dimension(3,3) :: math_eigenvectorBasisSym33
@ -2103,7 +2166,7 @@ end function math_eigenvectorBasisSym33
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief logarithm eigenvector basis of symmetric 33 matrix m !> @brief logarithm eigenvector basis of symmetric 33 matrix m
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function math_eigenvectorBasisSym33_log(m) pure function math_eigenvectorBasisSym33_log(m)
implicit none implicit none
real(pReal), dimension(3,3) :: math_eigenvectorBasisSym33_log real(pReal), dimension(3,3) :: math_eigenvectorBasisSym33_log
@ -2164,6 +2227,7 @@ function math_eigenvectorBasisSym33_log(m)
end function math_eigenvectorBasisSym33_log end function math_eigenvectorBasisSym33_log
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief rotational part from polar decomposition of 33 tensor m !> @brief rotational part from polar decomposition of 33 tensor m
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -2608,13 +2672,12 @@ pure function math_rotate_forward3333(tensor,rot_tensor)
real(pReal), dimension(3,3,3,3), intent(in) :: tensor real(pReal), dimension(3,3,3,3), intent(in) :: tensor
integer(pInt) :: i,j,k,l,m,n,o,p integer(pInt) :: i,j,k,l,m,n,o,p
math_rotate_forward3333= 0.0_pReal 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 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
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) = math_rotate_forward3333(i,j,k,l) & = math_rotate_forward3333(i,j,k,l) &
+ rot_tensor(i,m) * rot_tensor(j,n) & + rot_tensor(i,m) * rot_tensor(j,n) * rot_tensor(k,o) * rot_tensor(l,p) * tensor(m,n,o,p)
* rot_tensor(k,o) * rot_tensor(l,p) * tensor(m,n,o,p)
enddo; enddo; enddo; enddo; enddo; enddo; enddo; enddo enddo; enddo; enddo; enddo; enddo; enddo; enddo; enddo
end function math_rotate_forward3333 end function math_rotate_forward3333