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

added new math function "math_deviatoric33" which returns the deviatoric part of a 3x3 tensor 

renamed some math functions, so that we have a universal naming scheme: for matrix multiplications use an "x" (e.g. math_mul33x3); don't use the "x" to describe the shape of the tensor that the function is applied to (e.g. math_invert33 instead of math_invert3x3)
This commit is contained in:
Christoph Kords 2012-01-26 13:50:00 +00:00
parent df931890e0
commit 1330576a01
9 changed files with 233 additions and 225 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
code/CPFEM.f90
View File

@ -259,8 +259,8 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP,
restartWrite
use math, only: math_identity2nd, &
math_mul33x33, &
math_det3x3, &
math_transpose3x3, &
math_det33, &
math_transpose33, &
math_I3, &
math_Mandel3333to66, &
math_Mandel66to3333, &
@ -493,8 +493,8 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP,
if (debug_verbosity > 0) then
!$OMP CRITICAL (write2out)
write(6,'(a,x,i8,x,i2)') '<< CPFEM >> OUTDATED at element ip',cp_en,IP
write(6,'(a,/,3(12(x),3(f10.6,x),/))') '<< CPFEM >> FFN1 old:',math_transpose3x3(materialpoint_F(1:3,1:3,IP,cp_en))
write(6,'(a,/,3(12(x),3(f10.6,x),/))') '<< CPFEM >> FFN1 now:',math_transpose3x3(ffn1)
write(6,'(a,/,3(12(x),3(f10.6,x),/))') '<< CPFEM >> FFN1 old:',math_transpose33(materialpoint_F(1:3,1:3,IP,cp_en))
write(6,'(a,/,3(12(x),3(f10.6,x),/))') '<< CPFEM >> FFN1 now:',math_transpose33(ffn1)
!$OMP END CRITICAL (write2out)
endif
outdatedFFN1 = .true.
@ -569,8 +569,8 @@ subroutine CPFEM_general(mode, coords, ffn, ffn1, Temperature, dt, element, IP,
CPFEM_dcsde(1:6,1:6,IP,cp_en) = CPFEM_odd_jacobian * math_identity2nd(6)
else
! translate from P to CS
Kirchhoff = math_mul33x33(materialpoint_P(1:3,1:3,IP, cp_en), math_transpose3x3(materialpoint_F(1:3,1:3,IP,cp_en)))
J_inverse = 1.0_pReal / math_det3x3(materialpoint_F(1:3,1:3,IP,cp_en))
Kirchhoff = math_mul33x33(materialpoint_P(1:3,1:3,IP, cp_en), math_transpose33(materialpoint_F(1:3,1:3,IP,cp_en)))
J_inverse = 1.0_pReal / math_det33(materialpoint_F(1:3,1:3,IP,cp_en))
CPFEM_cs(1:6,IP,cp_en) = math_Mandel33to6(J_inverse * Kirchhoff)
! translate from dP/dF to dCS/dE

50
code/DAMASK_spectral.f90
View File

@ -411,13 +411,13 @@ program DAMASK_spectral
else
print '(a)','deformation gradient rate:'
endif
print '(3(3(f12.6,x)/)$)', merge(math_transpose3x3(bc(loadcase)%deformation),&
print '(3(3(f12.6,x)/)$)', merge(math_transpose33(bc(loadcase)%deformation),&
reshape(spread(DAMASK_NaN,1,9),(/3,3/)),transpose(bc(loadcase)%maskDeformation))
print '(a,/,3(3(f12.6,x)/)$)','stress / GPa:',1e-9*merge(math_transpose3x3(bc(loadcase)%stress)&
print '(a,/,3(3(f12.6,x)/)$)','stress / GPa:',1e-9*merge(math_transpose33(bc(loadcase)%stress)&
,reshape(spread(DAMASK_NaN,1,9),(/3,3/))&
,transpose(bc(loadcase)%maskStress))
if (any(bc(loadcase)%rotation /= math_I3)) &
print '(a,3(3(f12.6,x)/)$)','rotation of loadframe:',math_transpose3x3(bc(loadcase)%rotation)
print '(a,3(3(f12.6,x)/)$)','rotation of loadframe:',math_transpose33(bc(loadcase)%rotation)
print '(a,f12.6)','temperature:',bc(loadcase)%temperature
print '(a,f12.6)','time: ',bc(loadcase)%time
print '(a,i5)' ,'increments: ',bc(loadcase)%incs
@ -428,9 +428,9 @@ program DAMASK_spectral
if (any(bc(loadcase)%maskStress .and. transpose(bc(loadcase)%maskStress) .and. &
reshape((/.false.,.true.,.true.,.true.,.false.,.true.,.true.,.true.,.false./),(/3,3/)))) &
errorID = 38_pInt ! no rotation is allowed by stress BC
if (any(abs(math_mul33x33(bc(loadcase)%rotation,math_transpose3x3(bc(loadcase)%rotation))&
if (any(abs(math_mul33x33(bc(loadcase)%rotation,math_transpose33(bc(loadcase)%rotation))&
-math_I3) > reshape(spread(rotation_tol,1,9),(/3,3/)))&
.or. abs(math_det3x3(bc(loadcase)%rotation)) > 1.0_pReal + rotation_tol)&
.or. abs(math_det33(bc(loadcase)%rotation)) > 1.0_pReal + rotation_tol)&
errorID = 46_pInt ! given rotation matrix contains strain
if (bc(loadcase)%time < 0.0_pReal) errorID = 34_pInt ! negative time increment
if (bc(loadcase)%incs < 1_pInt) errorID = 35_pInt ! non-positive incs count
@ -634,7 +634,7 @@ program DAMASK_spectral
do l = 1_pInt ,3_pInt; do m = 1_pInt,3_pInt
xiDyad(l,m) = xi(l,i,j,k)*xi(m,i,j,k)
enddo; enddo
temp33_Real = math_inv3x3(math_mul3333xx33(c0_reference, xiDyad))
temp33_Real = math_inv33(math_mul3333xx33(c0_reference, xiDyad))
else
xiDyad = 0.0_pReal
@ -704,10 +704,10 @@ program DAMASK_spectral
temp33_Real = defgrad(i,j,k,1:3,1:3)
if (bc(loadcase)%velGradApplied) & ! use velocity gradient to calculate new deformation gradient (if not guessing)
fDot = math_mul33x33(bc(loadcase)%deformation,&
math_rotate_forward3x3(defgradold(i,j,k,1:3,1:3),bc(loadcase)%rotation))
math_rotate_forward33(defgradold(i,j,k,1:3,1:3),bc(loadcase)%rotation))
defgrad(i,j,k,1:3,1:3) = defgrad(i,j,k,1:3,1:3) & ! decide if guessing along former trajectory or apply homogeneous addon
+ guessmode * (defgrad(i,j,k,1:3,1:3) - defgradold(i,j,k,1:3,1:3))& ! guessing...
+ math_rotate_backward3x3((1.0_pReal-guessmode) * mask_defgrad * fDot,&
+ math_rotate_backward33((1.0_pReal-guessmode) * mask_defgrad * fDot,&
bc(loadcase)%rotation) *timeinc ! apply the prescribed value where deformation is given if not guessing
defgradold(i,j,k,1:3,1:3) = temp33_Real
enddo; enddo; enddo
@ -725,7 +725,7 @@ program DAMASK_spectral
!--------------------------------------------------------------------------------------------------
! calculate reduced compliance
c_prev = math_rotate_forward3x3x3x3(c_current*wgt,bc(loadcase)%rotation) ! calculate stiffness from former inc
c_prev = math_rotate_forward3333(c_current*wgt,bc(loadcase)%rotation) ! calculate stiffness from former inc
if(size_reduced > 0_pInt) then ! calculate compliance in case stress BC is applied
c_prev99 = math_Plain3333to99(c_prev)
k = 0_pInt ! build reduced stiffness
@ -788,7 +788,7 @@ program DAMASK_spectral
defgrad_av_lab(m,n) = sum(defgrad(1:res(1),1:res(2),1:res(3),m,n)) * wgt
enddo; enddo
print '(a,/,3(3(f12.7,x)/)$)', 'deformation gradient:',&
math_transpose3x3(math_rotate_forward3x3(defgrad_av_lab,bc(loadcase)%rotation))
math_transpose33(math_rotate_forward33(defgrad_av_lab,bc(loadcase)%rotation))
print '(a)', ''
print '(a)', '... update stress field P(F) ................................'
@ -874,8 +874,8 @@ program DAMASK_spectral
!--------------------------------------------------------------------------------------------------
! calculating RMS divergence criterion in Fourier space
p_hat_avg = sqrt(maxval (math_eigenvalues3x3(math_mul33x33(real(tensorField_complex(1,1,1,1:3,1:3)),& ! L_2 norm of average stress (freq 0,0,0) in fourier space,
math_transpose3x3(real(tensorField_complex(1,1,1,1:3,1:3))))))) ! ignore imaginary part as it is always zero for real only input
p_hat_avg = sqrt(maxval (math_eigenvalues33(math_mul33x33(real(tensorField_complex(1,1,1,1:3,1:3)),& ! L_2 norm of average stress (freq 0,0,0) in fourier space,
math_transpose33(real(tensorField_complex(1,1,1,1:3,1:3))))))) ! ignore imaginary part as it is always zero for real only input
err_div_RMS = 0.0_pReal
do k = 1_pInt, res(3); do j = 1_pInt, res(2)
do i = 2_pInt, res1_red -1_pInt ! Has somewhere a conj. complex counterpart. Therefore count it twice.
@ -914,8 +914,8 @@ program DAMASK_spectral
err_real_div_RMS = err_real_div_RMS + sum(divergence_real(i,j,k,1:3)**2.0_pReal) ! avg of L_2 norm of div(stress) in real space
err_real_div_max = max(err_real_div_max, sqrt(sum(divergence_real(i,j,k,1:3)**2.0_pReal))) ! maximum of L two norm of div(stress) in real space
enddo; enddo; enddo
p_real_avg = sqrt(maxval (math_eigenvalues3x3(math_mul33x33(pstress_av_lab,& ! L_2 norm of average stress in real space,
math_transpose3x3(pstress_av_lab)))))
p_real_avg = sqrt(maxval (math_eigenvalues33(math_mul33x33(pstress_av_lab,& ! L_2 norm of average stress in real space,
math_transpose33(pstress_av_lab)))))
err_real_div_RMS = sqrt(wgt*err_real_div_RMS)*correctionFactor ! RMS in real space
err_real_div_max = err_real_div_max *correctionFactor
err_div_max = err_div_max *correctionFactor
@ -943,7 +943,7 @@ program DAMASK_spectral
do l = 1_pInt,3_pInt; do m = 1_pInt,3_pInt
xiDyad(l,m) = xi(l,i,j,k)*xi(m,i,j,k)
enddo; enddo
temp33_Real = math_inv3x3(math_mul3333xx33(c0_reference, xiDyad))
temp33_Real = math_inv33(math_mul3333xx33(c0_reference, xiDyad))
else
xiDyad = 0.0_pReal
temp33_Real = 0.0_pReal
@ -996,9 +996,9 @@ program DAMASK_spectral
temp33_Real = 0.0_pReal
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
maxCorrectionSym = max(maxCorrectionSym,&
maxval(math_symmetric3x3(tensorField_real(i,j,k,1:3,1:3))))
maxval(math_symmetric33(tensorField_real(i,j,k,1:3,1:3))))
maxCorrectionSkew = max(maxCorrectionSkew,&
maxval(math_skew3x3(tensorField_real(i,j,k,1:3,1:3))))
maxval(math_skew33(tensorField_real(i,j,k,1:3,1:3))))
temp33_Real = temp33_Real + tensorField_real(i,j,k,1:3,1:3)
enddo; enddo; enddo
print '(a,x,es10.4)' , 'max symmetrix correction of deformation:',&
@ -1006,8 +1006,8 @@ program DAMASK_spectral
print '(a,x,es10.4)' , 'max skew correction of deformation:',&
maxCorrectionSkew*wgt
print '(a,x,es10.4)' , 'max sym/skew of avg correction: ',&
maxval(math_symmetric3x3(temp33_real))/&
maxval(math_skew3x3(temp33_real))
maxval(math_symmetric33(temp33_real))/&
maxval(math_skew33(temp33_real))
endif
!--------------------------------------------------------------------------------------------------
@ -1028,8 +1028,8 @@ program DAMASK_spectral
!--------------------------------------------------------------------------------------------------
! stress BC handling
pstress_av = math_rotate_forward3x3(pstress_av_lab,bc(loadcase)%rotation)
print '(a,/,3(3(f12.7,x)/)$)', 'Piola-Kirchhoff stress / MPa: ',math_transpose3x3(pstress_av)/1.e6
pstress_av = math_rotate_forward33(pstress_av_lab,bc(loadcase)%rotation)
print '(a,/,3(3(f12.7,x)/)$)', 'Piola-Kirchhoff stress / MPa: ',math_transpose33(pstress_av)/1.e6
if(size_reduced > 0_pInt) then ! calculate stress BC if applied
err_stress = maxval(abs(mask_stress * (pstress_av - bc(loadcase)%stress))) ! maximum deviaton (tensor norm not applicable)
@ -1039,14 +1039,14 @@ program DAMASK_spectral
print '(a,es10.4,a,f6.2)', 'error stress = ',err_stress, ', ', err_stress/err_stress_tol
defgradAimCorr = - math_mul3333xx33(s_prev, ((pstress_av - bc(loadcase)%stress))) ! residual on given stress components
defgradAim = defgradAim + defgradAimCorr
print '(a,/,3(3(f12.7,x)/)$)', 'new deformation aim: ', math_transpose3x3(defgradAim)
print '(a,x,es10.4)' , 'with determinant: ', math_det3x3(defgradAim)
print '(a,/,3(3(f12.7,x)/)$)', 'new deformation aim: ', math_transpose33(defgradAim)
print '(a,x,es10.4)' , 'with determinant: ', math_det33(defgradAim)
else
err_stress_tol = 0.0_pReal
endif
! homogeneous correction towards avg deformation gradient -------------
defgradAim_lab = math_rotate_backward3x3(defgradAim,bc(loadcase)%rotation) ! boundary conditions from load frame into lab (Fourier) frame
defgradAim_lab = math_rotate_backward33(defgradAim,bc(loadcase)%rotation) ! boundary conditions from load frame into lab (Fourier) frame
do m = 1_pInt,3_pInt; do n = 1_pInt,3_pInt
defgrad(1:res(1),1:res(2),1:res(3),m,n) = &
defgrad(1:res(1),1:res(2),1:res(3),m,n) + (defgradAim_lab(m,n) - defgrad_av_lab(m,n)) ! anticipated target minus current state
@ -1058,7 +1058,7 @@ program DAMASK_spectral
defgradDetMax = -huge(1.0_pReal)
defgradDetMin = +huge(1.0_pReal)
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
defgradDet = math_det3x3(defgrad(i,j,k,1:3,1:3))
defgradDet = math_det33(defgrad(i,j,k,1:3,1:3))
defgradDetMax = max(defgradDetMax,defgradDet)
defgradDetMin = min(defgradDetMin,defgradDet)
enddo; enddo; enddo

10
code/constitutive_dislotwin.f90
View File

@ -932,7 +932,7 @@ subroutine constitutive_dislotwin_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,Temperat
!*********************************************************************
use prec, only: pReal,pInt,p_vec
use math, only: math_Plain3333to99, math_Mandel6to33, math_Mandel33to6, &
math_spectralDecompositionSym3x3, math_tensorproduct, math_symmetric3x3,math_mul33x3
math_spectralDecompositionSym33, math_tensorproduct, math_symmetric33,math_mul33x3
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance
use lattice, only: lattice_Sslip,lattice_Sslip_v,lattice_Stwin,lattice_Stwin_v,lattice_maxNslipFamily,lattice_maxNtwinFamily, &
@ -1040,12 +1040,12 @@ enddo
if(constitutive_dislotwin_sbVelocity(myInstance) /= 0.0_pReal) then
gdot_sb = 0.0_pReal
dgdot_dtausb = 0.0_pReal
call math_spectralDecompositionSym3x3(math_Mandel6to33(Tstar_v),eigValues,eigVectors, error)
call math_spectralDecompositionSym33(math_Mandel6to33(Tstar_v),eigValues,eigVectors, error)
do j = 1,6
sb_s = 0.5_pReal*sqrt(2.0_pReal)*math_mul33x3(eigVectors,sb_sComposition(1:3,j))
sb_m = 0.5_pReal*sqrt(2.0_pReal)*math_mul33x3(eigVectors,sb_mComposition(1:3,j))
sb_Smatrix = math_tensorproduct(sb_s,sb_m)
constitutive_dislotwin_sbSv(1:6,j,g,ip,el) = math_Mandel33to6(math_symmetric3x3(sb_Smatrix))
constitutive_dislotwin_sbSv(1:6,j,g,ip,el) = math_Mandel33to6(math_symmetric33(sb_Smatrix))
!* Calculation of Lp
!* Resolved shear stress on shear banding system
@ -1353,7 +1353,7 @@ function constitutive_dislotwin_postResults(Tstar_v,Temperature,dt,state,g,ip,el
!* - el : current element *
!*********************************************************************
use prec, only: pReal,pInt,p_vec
use math, only: pi,math_Mandel6to33, math_spectralDecompositionSym3x3
use math, only: pi,math_Mandel6to33, math_spectralDecompositionSym33
use mesh, only: mesh_NcpElems,mesh_maxNips
use material, only: homogenization_maxNgrains,material_phase,phase_constitutionInstance,phase_Noutput
use lattice, only: lattice_Sslip_v,lattice_Stwin_v,lattice_maxNslipFamily,lattice_maxNtwinFamily, &
@ -1387,7 +1387,7 @@ c = 0_pInt
constitutive_dislotwin_postResults = 0.0_pReal
!* Spectral decomposition of stress
call math_spectralDecompositionSym3x3(math_Mandel6to33(Tstar_v),eigValues,eigVectors, error)
call math_spectralDecompositionSym33(math_Mandel6to33(Tstar_v),eigValues,eigVectors, error)
do o = 1,phase_Noutput(material_phase(g,ip,el))
select case(constitutive_dislotwin_output(o,myInstance))

3
code/constitutive_j2.f90
View File

@ -326,8 +326,7 @@ pure subroutine constitutive_j2_LpAndItsTangent(Lp, dLp_dTstar_99, Tstar_dev_v,
p_vec
use math, only: math_mul6x6, &
math_Mandel6to33, &
math_Plain3333to99, &
math_spectral1
math_Plain3333to99
use lattice, only: lattice_Sslip, &
lattice_Sslip_v
use mesh, only: mesh_NcpElems, &

67
code/constitutive_nonlocal.f90
View File

@ -144,7 +144,7 @@ use prec, only: pInt, pReal
use math, only: math_Mandel3333to66, &
math_Voigt66to3333, &
math_mul3x3, &
math_transpose3x3
math_transpose33
use IO, only: IO_lc, &
IO_getTag, &
IO_isBlank, &
@ -720,9 +720,9 @@ do i = 1,maxNinstance
!*** rotation matrix from lattice configuration to slip system
constitutive_nonlocal_lattice2slip(1:3,1:3,s1,i) &
= math_transpose3x3( reshape((/ lattice_sd(1:3, constitutive_nonlocal_slipSystemLattice(s1,i), myStructure), &
-lattice_st(1:3, constitutive_nonlocal_slipSystemLattice(s1,i), myStructure), &
lattice_sn(1:3, constitutive_nonlocal_slipSystemLattice(s1,i), myStructure)/), (/3,3/)))
= math_transpose33( reshape((/ lattice_sd(1:3, constitutive_nonlocal_slipSystemLattice(s1,i), myStructure), &
-lattice_st(1:3, constitutive_nonlocal_slipSystemLattice(s1,i), myStructure), &
lattice_sn(1:3, constitutive_nonlocal_slipSystemLattice(s1,i), myStructure)/), (/3,3/)))
enddo
enddo
@ -886,9 +886,9 @@ use math, only: math_Mandel33to6, &
math_mul33x3, &
math_mul3x3, &
math_norm3, &
math_inv3x3, &
math_invert3x3, &
math_transpose3x3, &
math_inv33, &
math_invert33, &
math_transpose33, &
pi
use debug, only: debug_verbosity, &
debug_selectiveDebugger, &
@ -1034,8 +1034,8 @@ forall (s = 1:ns) &
tauBack = 0.0_pReal
if (.not. phase_localConstitution(phase)) then
call math_invert3x3(Fe, invFe, detFe, inversionError)
call math_invert3x3(Fp, invFp, detFp, inversionError)
call math_invert33(Fe, invFe, detFe, inversionError)
call math_invert33(Fp, invFp, detFp, inversionError)
ipCoords = mesh_ipCenterOfGravity(1:3,ip,el)
rhoExcess(1,1:ns) = rhoSgl(1:ns,1) - rhoSgl(1:ns,2)
rhoExcess(2,1:ns) = rhoSgl(1:ns,3) - rhoSgl(1:ns,4)
@ -1048,7 +1048,7 @@ if (.not. phase_localConstitution(phase)) then
do n = 1,FE_NipNeighbors(mesh_element(2,el))
neighboring_el = mesh_ipNeighborhood(1,n,ip,el)
neighboring_ip = mesh_ipNeighborhood(2,n,ip,el)
areaNormal_latticeConf(1:3,n) = detFp * math_mul33x3(math_transpose3x3(invFp), mesh_ipAreaNormal(1:3,n,ip,el)) ! calculate the normal of the interface in lattice configuration
areaNormal_latticeConf(1:3,n) = detFp * math_mul33x3(math_transpose33(invFp), mesh_ipAreaNormal(1:3,n,ip,el)) ! calculate the normal of the interface in lattice configuration
areaNormal_latticeConf(1:3,n) = areaNormal_latticeConf(1:3,n) / math_norm3(areaNormal_latticeConf(1:3,n)) ! normalize the surface normal to unit length
if (neighboring_el > 0 .and. neighboring_ip > 0) then
neighboring_phase = material_phase(g,neighboring_ip,neighboring_el)
@ -1143,7 +1143,7 @@ if (.not. phase_localConstitution(phase)) then
sampledPoint(1:3,1) = + gradientDistanceInter(1) * m(1:3,s,c)
sampledPoint(1:3,2) = - gradientDistanceInter(2) * m(1:3,s,c)
do side = 1,2
rhoExcessAtSampledPoint(side) = math_mul3x3(math_mul33x3(math_inv3x3(connections(1:3,1:3,side)), &
rhoExcessAtSampledPoint(side) = math_mul3x3(math_mul33x3(math_inv33(connections(1:3,1:3,side)), &
rhoExcessDifferences(1:3,side)), &
sampledPoint(1:3,side)) &
+ rhoExcess(c,s)
@ -1207,8 +1207,6 @@ subroutine constitutive_nonlocal_kinetics(v, tau, c, Temperature, state, g, ip,
use prec, only: pReal, &
pInt, &
p_vec
use math, only: math_mul6x6, &
math_Mandel6to33
use debug, only: debug_verbosity, &
debug_selectiveDebugger, &
debug_g, &
@ -1358,8 +1356,7 @@ use prec, only: pReal, &
pInt, &
p_vec
use math, only: math_Plain3333to99, &
math_mul6x6, &
math_Mandel6to33
math_mul6x6
use debug, only: debug_verbosity, &
debug_selectiveDebugger, &
debug_g, &
@ -1515,11 +1512,9 @@ use math, only: math_norm3, &
math_mul3x3, &
math_mul33x3, &
math_mul33x33, &
math_inv3x3, &
math_det3x3, &
math_Mandel6to33, &
math_QuaternionDisorientation, &
math_qRot, &
math_inv33, &
math_det33, &
math_transpose33, &
pi
use mesh, only: mesh_NcpElems, &
mesh_maxNips, &
@ -1825,9 +1820,9 @@ if (.not. phase_localConstitution(material_phase(g,ip,el))) then
forall (t = 1:4) &
neighboring_fluxdensity(1:ns,t) = state(g,neighboring_ip,neighboring_el)%p((t-1)*ns+1:t*ns) &
* state(g,neighboring_ip,neighboring_el)%p((12+t)*ns+1:(13+t)*ns)
normal_neighbor2me_defConf = math_det3x3(Favg) &
* math_mul33x3(math_inv3x3(transpose(Favg)), mesh_ipAreaNormal(1:3,neighboring_n,neighboring_ip,neighboring_el)) ! calculate the normal of the interface in (average) deformed configuration (now pointing from my neighbor to me!!!)
normal_neighbor2me = math_mul33x3(transpose(neighboring_Fe), normal_neighbor2me_defConf) / math_det3x3(neighboring_Fe) ! interface normal in the lattice configuration of my neighbor
normal_neighbor2me_defConf = math_det33(Favg) &
* math_mul33x3(math_inv33(transpose(Favg)), mesh_ipAreaNormal(1:3,neighboring_n,neighboring_ip,neighboring_el)) ! calculate the normal of the interface in (average) deformed configuration (now pointing from my neighbor to me!!!)
normal_neighbor2me = math_mul33x3(transpose(neighboring_Fe), normal_neighbor2me_defConf) / math_det33(neighboring_Fe) ! interface normal in the lattice configuration of my neighbor
area = mesh_ipArea(neighboring_n,neighboring_ip,neighboring_el) * math_norm3(normal_neighbor2me)
normal_neighbor2me = normal_neighbor2me / math_norm3(normal_neighbor2me) ! normalize the surface normal to unit length
do s = 1,ns
@ -1864,8 +1859,8 @@ if (.not. phase_localConstitution(material_phase(g,ip,el))) then
endif
if (considerLeavingFlux) then
normal_me2neighbor_defConf = math_det3x3(Favg) * math_mul33x3(math_inv3x3(transpose(Favg)), mesh_ipAreaNormal(1:3,n,ip,el)) ! calculate the normal of the interface in (average) deformed configuration (pointing from me to my neighbor!!!)
normal_me2neighbor = math_mul33x3(transpose(my_Fe), normal_me2neighbor_defConf) / math_det3x3(my_Fe) ! interface normal in my lattice configuration
normal_me2neighbor_defConf = math_det33(Favg) * math_mul33x3(math_inv33(math_transpose33(Favg)), mesh_ipAreaNormal(1:3,n,ip,el)) ! calculate the normal of the interface in (average) deformed configuration (pointing from me to my neighbor!!!)
normal_me2neighbor = math_mul33x3(math_transpose33(my_Fe), normal_me2neighbor_defConf) / math_det33(my_Fe) ! interface normal in my lattice configuration
area = mesh_ipArea(n,ip,el) * math_norm3(normal_me2neighbor)
normal_me2neighbor = normal_me2neighbor / math_norm3(normal_me2neighbor) ! normalize the surface normal to unit length
do s = 1,ns
@ -2187,8 +2182,8 @@ use prec, only: pReal, &
p_vec
use math, only: math_mul33x33, &
math_mul33x3, &
math_invert3x3, &
math_transpose3x3, &
math_invert33, &
math_transpose33, &
pi
use mesh, only: mesh_NcpElems, &
mesh_maxNips, &
@ -2296,7 +2291,7 @@ forall (t = 5:8) &
constitutive_nonlocal_dislocationstress = 0.0_pReal
if (.not. phase_localConstitution(phase)) then
call math_invert3x3(Fe(1:3,1:3,1,ip,el), invFe, detFe, inversionError)
call math_invert33(Fe(1:3,1:3,1,ip,el), invFe, detFe, inversionError)
! if (inversionError) then
! return
! endif
@ -2330,7 +2325,7 @@ ipLoop: do neighboring_ip = 1,FE_Nips(mesh_element(2,neighboring_el))
neighboring_instance = phase_constitutionInstance(neighboring_phase)
neighboring_latticeStruct = constitutive_nonlocal_structure(neighboring_instance)
neighboring_ns = constitutive_nonlocal_totalNslip(neighboring_instance)
call math_invert3x3(Fe(1:3,1:3,1,neighboring_ip,neighboring_el), neighboring_invFe, detFe, inversionError)
call math_invert33(Fe(1:3,1:3,1,neighboring_ip,neighboring_el), neighboring_invFe, detFe, inversionError)
! if (inversionError) then
! return
! endif
@ -2479,7 +2474,7 @@ ipLoop: do neighboring_ip = 1,FE_Nips(mesh_element(2,neighboring_el))
/ (4.0_pReal * pi * (1.0_pReal - nu)) &
* mesh_ipVolume(neighboring_ip,neighboring_el) / segmentLength ! reference volume is used here (according to the segment length calculation)
Tdislo_neighboringLattice = Tdislo_neighboringLattice &
+ math_mul33x33(math_transpose3x3(constitutive_nonlocal_lattice2slip(1:3,1:3,s,neighboring_instance)), &
+ math_mul33x33(math_transpose33(constitutive_nonlocal_lattice2slip(1:3,1:3,s,neighboring_instance)), &
math_mul33x33(sigma, constitutive_nonlocal_lattice2slip(1:3,1:3,s,neighboring_instance)))
enddo ! slip system loop
@ -2506,7 +2501,7 @@ ipLoop: do neighboring_ip = 1,FE_Nips(mesh_element(2,neighboring_el))
sigma(3,1) = sigma(1,3)
Tdislo_neighboringLattice = Tdislo_neighboringLattice &
+ math_mul33x33(math_transpose3x3(constitutive_nonlocal_lattice2slip(1:3,1:3,s,instance)), &
+ math_mul33x33(math_transpose33(constitutive_nonlocal_lattice2slip(1:3,1:3,s,instance)), &
math_mul33x33(sigma, constitutive_nonlocal_lattice2slip(1:3,1:3,s,instance)))
enddo ! slip system loop
@ -2525,7 +2520,7 @@ ipLoop: do neighboring_ip = 1,FE_Nips(mesh_element(2,neighboring_el))
constitutive_nonlocal_dislocationstress = constitutive_nonlocal_dislocationstress &
+ math_mul33x33(neighboringLattice2myLattice, &
math_mul33x33(Tdislo_neighboringLattice, &
math_transpose3x3(neighboringLattice2myLattice)))
math_transpose33(neighboringLattice2myLattice)))
enddo ipLoop
enddo ! element loop
@ -2543,15 +2538,9 @@ function constitutive_nonlocal_postResults(Tstar_v, Fe, Temperature, dt, state,
use prec, only: pReal, &
pInt, &
p_vec
use math, only: math_norm3, &
math_mul6x6, &
math_mul3x3, &
use math, only: math_mul6x6, &
math_mul33x3, &
math_mul33x33, &
math_inv3x3, &
math_det3x3, &
math_Mandel6to33, &
math_transpose3x3, &
pi
use mesh, only: mesh_NcpElems, &
mesh_maxNips, &

91
code/crystallite.f90
View File

@ -112,8 +112,8 @@ use numerics, only: subStepSizeCryst, &
stepIncreaseCryst
use math, only: math_I3, &
math_EulerToR, &
math_inv3x3, &
math_transpose3x3, &
math_inv33, &
math_transpose33, &
math_mul33xx33, &
math_mul33x33
use FEsolving, only: FEsolving_execElem, &
@ -322,7 +322,7 @@ close(file)
crystallite_F0(1:3,1:3,g,i,e) = math_I3
crystallite_localConstitution(g,i,e) = phase_localConstitution(material_phase(g,i,e))
!$OMP FLUSH(crystallite_Fp0)
crystallite_Fe(1:3,1:3,g,i,e) = math_transpose3x3(crystallite_Fp0(1:3,1:3,g,i,e))
crystallite_Fe(1:3,1:3,g,i,e) = math_transpose33(crystallite_Fp0(1:3,1:3,g,i,e))
enddo
enddo
enddo
@ -467,13 +467,12 @@ use debug, only: debug_verbosity, &
debug_g, &
debug_CrystalliteLoopDistribution
use IO, only: IO_warning
use math, only: math_inv3x3, &
math_transpose3x3, &
use math, only: math_inv33, &
math_transpose33, &
math_mul33x33, &
math_mul66x6, &
math_Mandel6to33, &
math_Mandel33to6, &
math_transpose3x3, &
math_I3
use FEsolving, only: FEsolving_execElem, &
FEsolving_execIP
@ -542,11 +541,11 @@ if (debug_verbosity > 4 .and. debug_e > 0 .and. debug_e <= mesh_NcpElems &
write (6,'(a,i8,x,i2,x,i3)') '<< CRYST >> crystallite start at el ip g ', debug_e, debug_i, debug_g
write (6,'(a,/,12(x),f14.9)') '<< CRYST >> Temp0', crystallite_partionedTemperature0(debug_g,debug_i,debug_e)
write (6,'(a,/,3(12(x),3(f14.9,x)/))') '<< CRYST >> F0 ', &
math_transpose3x3(crystallite_partionedF0(1:3,1:3,debug_g,debug_i,debug_e))
math_transpose33(crystallite_partionedF0(1:3,1:3,debug_g,debug_i,debug_e))
write (6,'(a,/,3(12(x),3(f14.9,x)/))') '<< CRYST >> Fp0', &
math_transpose3x3(crystallite_partionedFp0(1:3,1:3,debug_g,debug_i,debug_e))
math_transpose33(crystallite_partionedFp0(1:3,1:3,debug_g,debug_i,debug_e))
write (6,'(a,/,3(12(x),3(f14.9,x)/))') '<< CRYST >> Lp0', &
math_transpose3x3(crystallite_partionedLp0(1:3,1:3,debug_g,debug_i,debug_e))
math_transpose33(crystallite_partionedLp0(1:3,1:3,debug_g,debug_i,debug_e))
!$OMP END CRITICAL (write2out)
endif
@ -566,7 +565,7 @@ crystallite_subStep = 0.0_pReal
crystallite_subF0(1:3,1:3,g,i,e) = crystallite_partionedF0(1:3,1:3,g,i,e) ! ...def grad
crystallite_subTstar0_v(1:6,g,i,e) = crystallite_partionedTstar0_v(1:6,g,i,e) !...2nd PK stress
crystallite_subFe0(1:3,1:3,g,i,e) = math_mul33x33(crystallite_subF0(1:3,1:3,g,i,e), &
math_inv3x3(crystallite_subFp0(1:3,1:3,g,i,e))) ! only needed later on for stiffness calculation
math_inv33(crystallite_subFp0(1:3,1:3,g,i,e))) ! only needed later on for stiffness calculation
crystallite_subFrac(g,i,e) = 0.0_pReal
crystallite_subStep(g,i,e) = 1.0_pReal/subStepSizeCryst
@ -633,7 +632,7 @@ do while (any(crystallite_subStep(:,:,FEsolving_execELem(1):FEsolving_execElem(2
crystallite_subStep(g,i,e) = subStepSizeCryst * crystallite_subStep(g,i,e) ! cut step in half and restore...
crystallite_Temperature(g,i,e) = crystallite_subTemperature0(g,i,e) ! ...temperature
crystallite_Fp(1:3,1:3,g,i,e) = crystallite_subFp0(1:3,1:3,g,i,e) ! ...plastic def grad
crystallite_invFp(1:3,1:3,g,i,e) = math_inv3x3(crystallite_Fp(1:3,1:3,g,i,e))
crystallite_invFp(1:3,1:3,g,i,e) = math_inv33(crystallite_Fp(1:3,1:3,g,i,e))
crystallite_Lp(1:3,1:3,g,i,e) = crystallite_subLp0(1:3,1:3,g,i,e) ! ...plastic velocity grad
constitutive_state(g,i,e)%p = constitutive_subState0(g,i,e)%p ! ...microstructure
crystallite_Tstar_v(1:6,g,i,e) = crystallite_subTstar0_v(1:6,g,i,e) ! ...2nd PK stress
@ -698,7 +697,7 @@ enddo
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do g = 1,myNgrains
if (.not. crystallite_converged(g,i,e)) then ! respond fully elastically (might be not required due to becoming terminally ill anyway)
invFp = math_inv3x3(crystallite_partionedFp0(1:3,1:3,g,i,e))
invFp = math_inv33(crystallite_partionedFp0(1:3,1:3,g,i,e))
Fe_guess = math_mul33x33(crystallite_partionedF(1:3,1:3,g,i,e), invFp)
Tstar = math_Mandel6to33( math_mul66x6( 0.5_pReal*constitutive_homogenizedC(g,i,e), &
math_Mandel33to6( math_mul33x33(transpose(Fe_guess),Fe_guess) - math_I3 ) ) )
@ -709,9 +708,9 @@ enddo
.and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
write (6,'(a,i8,x,i2,x,i3)') '<< CRYST >> central solution of cryst_StressAndTangent at el ip g ',e,i,g
write (6,*)
write (6,'(a,/,3(12(x),3(f12.4,x)/))') '<< CRYST >> P / MPa', math_transpose3x3(crystallite_P(1:3,1:3,g,i,e)) / 1e6
write (6,'(a,/,3(12(x),3(f14.9,x)/))') '<< CRYST >> Fp', math_transpose3x3(crystallite_Fp(1:3,1:3,g,i,e))
write (6,'(a,/,3(12(x),3(f14.9,x)/))') '<< CRYST >> Lp', math_transpose3x3(crystallite_Lp(1:3,1:3,g,i,e))
write (6,'(a,/,3(12(x),3(f12.4,x)/))') '<< CRYST >> P / MPa', math_transpose33(crystallite_P(1:3,1:3,g,i,e)) / 1e6
write (6,'(a,/,3(12(x),3(f14.9,x)/))') '<< CRYST >> Fp', math_transpose33(crystallite_Fp(1:3,1:3,g,i,e))
write (6,'(a,/,3(12(x),3(f14.9,x)/))') '<< CRYST >> Lp', math_transpose33(crystallite_Lp(1:3,1:3,g,i,e))
write (6,*)
endif
#endif
@ -2605,11 +2604,11 @@ use math, only: math_mul33x33, &
math_mul33xx33, &
math_mul66x6, &
math_mul99x99, &
math_transpose3x3, &
math_inv3x3, &
math_invert3x3, &
math_transpose33, &
math_inv33, &
math_invert33, &
math_invert, &
math_det3x3, &
math_det33, &
math_norm33, &
math_I3, &
math_identity2nd, &
@ -2707,14 +2706,14 @@ Lpguess = crystallite_Lp(1:3,1:3,g,i,e) ! ... and tak
!* inversion of Fp_current...
invFp_current = math_inv3x3(Fp_current)
invFp_current = math_inv33(Fp_current)
if (all(invFp_current == 0.0_pReal)) then ! ... failed?
#ifndef _OPENMP
if (debug_verbosity > 4) then
write(6,'(a,i8,x,i2,x,i3)') '<< CRYST >> integrateStress failed on invFp_current inversion at el ip g ',e,i,g
if (debug_verbosity > 5 .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
write(6,*)
write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> invFp_new',math_transpose3x3(invFp_new(1:3,1:3))
write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> invFp_new',math_transpose33(invFp_new(1:3,1:3))
endif
endif
#endif
@ -2754,7 +2753,7 @@ LpLoop: do
endif
B = math_I3 - dt*Lpguess
BT = math_transpose3x3(B)
BT = math_transpose33(B)
AB = math_mul33x33(A,B)
BTA = math_mul33x33(BT,A)
@ -2787,8 +2786,8 @@ LpLoop: do
.and. numerics_integrationMode == 1_pInt) then
write(6,'(a,i3)') '<< CRYST >> iteration ', NiterationStress
write(6,*)
write(6,'(a,/,3(12(x),3(e20.7,x)/))') '<< CRYST >> Lp_constitutive', math_transpose3x3(Lp_constitutive)
write(6,'(a,/,3(12(x),3(e20.7,x)/))') '<< CRYST >> Lpguess', math_transpose3x3(Lpguess)
write(6,'(a,/,3(12(x),3(e20.7,x)/))') '<< CRYST >> Lp_constitutive', math_transpose33(Lp_constitutive)
write(6,'(a,/,3(12(x),3(e20.7,x)/))') '<< CRYST >> Lpguess', math_transpose33(Lpguess)
endif
#endif
@ -2903,10 +2902,10 @@ LpLoop: do
write(6,'(a,/,9(12(x),9(e15.3,x)/))') '<< CRYST >> dR_dLp',transpose(dR_dLp)
write(6,'(a,/,9(12(x),9(e15.3,x)/))') '<< CRYST >> dT_dLp',transpose(dT_dLp)
write(6,'(a,/,9(12(x),9(e15.3,x)/))') '<< CRYST >> dLp_dT_constitutive',transpose(dLp_dT_constitutive)
write(6,'(a,/,3(12(x),3(e20.7,x)/))') '<< CRYST >> AB',math_transpose3x3(AB)
write(6,'(a,/,3(12(x),3(e20.7,x)/))') '<< CRYST >> BTA',math_transpose3x3(BTA)
write(6,'(a,/,3(12(x),3(e20.7,x)/))') '<< CRYST >> Lp_constitutive',math_transpose3x3(Lp_constitutive)
write(6,'(a,/,3(12(x),3(e20.7,x)/))') '<< CRYST >> Lpguess',math_transpose3x3(Lpguess)
write(6,'(a,/,3(12(x),3(e20.7,x)/))') '<< CRYST >> AB',math_transpose33(AB)
write(6,'(a,/,3(12(x),3(e20.7,x)/))') '<< CRYST >> BTA',math_transpose33(BTA)
write(6,'(a,/,3(12(x),3(e20.7,x)/))') '<< CRYST >> Lp_constitutive',math_transpose33(Lp_constitutive)
write(6,'(a,/,3(12(x),3(e20.7,x)/))') '<< CRYST >> Lpguess',math_transpose33(Lpguess)
endif
endif
#endif
@ -2936,8 +2935,8 @@ enddo LpLoop
!* calculate new plastic and elastic deformation gradient
invFp_new = math_mul33x33(invFp_current,B)
invFp_new = invFp_new/math_det3x3(invFp_new)**(1.0_pReal/3.0_pReal) ! regularize by det
call math_invert3x3(invFp_new,Fp_new,det,error)
invFp_new = invFp_new/math_det33(invFp_new)**(1.0_pReal/3.0_pReal) ! regularize by det
call math_invert33(invFp_new,Fp_new,det,error)
if (error) then
#ifndef _OPENMP
if (debug_verbosity > 4) then
@ -2945,7 +2944,7 @@ if (error) then
' ; iteration ', NiterationStress
if (debug_verbosity > 5 .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger)) then
write(6,*)
write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> invFp_new',math_transpose3x3(invFp_new)
write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> invFp_new',math_transpose33(invFp_new)
endif
endif
#endif
@ -2957,7 +2956,7 @@ Fe_new = math_mul33x33(Fg_new,invFp_new) ! calc resu
!* add volumetric component to 2nd Piola-Kirchhoff stress and calculate 1st Piola-Kirchhoff stress
forall (n=1:3) Tstar_v(n) = Tstar_v(n) + p_hydro
crystallite_P(1:3,1:3,g,i,e) = math_mul33x33(Fe_new, math_mul33x33(math_Mandel6to33(Tstar_v), math_transpose3x3(invFp_new)))
crystallite_P(1:3,1:3,g,i,e) = math_mul33x33(Fe_new, math_mul33x33(math_Mandel6to33(Tstar_v), math_transpose33(invFp_new)))
!* store local values in global variables
@ -2975,12 +2974,12 @@ crystallite_integrateStress = .true.
#ifndef _OPENMP
if (debug_verbosity > 5 .and. ((e == debug_e .and. i == debug_i .and. g == debug_g) .or. .not. debug_selectiveDebugger) &
.and. numerics_integrationMode == 1_pInt) then
write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> P / MPa',math_transpose3x3(crystallite_P(1:3,1:3,g,i,e))/1e6
write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> P / MPa',math_transpose33(crystallite_P(1:3,1:3,g,i,e))/1e6
write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> Cauchy / MPa', &
math_mul33x33(crystallite_P(1:3,1:3,g,i,e), math_transpose3x3(Fg_new)) / 1e6 / math_det3x3(Fg_new)
math_mul33x33(crystallite_P(1:3,1:3,g,i,e), math_transpose33(Fg_new)) / 1e6 / math_det33(Fg_new)
write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> Fe Lp Fe^-1', &
math_transpose3x3(math_mul33x33(Fe_new, math_mul33x33(crystallite_Lp(1:3,1:3,g,i,e), math_inv3x3(Fe_new)))) ! transpose to get correct print out order
write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> Fp',math_transpose3x3(crystallite_Fp(1:3,1:3,g,i,e))
math_transpose33(math_mul33x33(Fe_new, math_mul33x33(crystallite_Lp(1:3,1:3,g,i,e), math_inv33(Fe_new)))) ! transpose to get correct print out order
write(6,'(a,/,3(12(x),3(f12.7,x)/))') '<< CRYST >> Fp',math_transpose33(crystallite_Fp(1:3,1:3,g,i,e))
endif
#endif
@ -3141,8 +3140,8 @@ function crystallite_postResults(&
use math, only: math_QuaternionToEuler, &
math_QuaternionToAxisAngle, &
math_mul33x33, &
math_transpose3x3, &
math_det3x3, &
math_transpose33, &
math_det33, &
math_I3, &
inDeg, &
math_Mandel6to33
@ -3191,7 +3190,7 @@ function crystallite_postResults(&
crystallite_postResults(c+1) = material_texture(g,i,e) ! textureID of grain
case ('volume')
mySize = 1_pInt
detF = math_det3x3(crystallite_partionedF(1:3,1:3,g,i,e)) ! V_current = det(F) * V_reference
detF = math_det33(crystallite_partionedF(1:3,1:3,g,i,e)) ! V_current = det(F) * V_reference
crystallite_postResults(c+1) = detF * mesh_ipVolume(i,e) / homogenization_Ngrains(mesh_element(3,e)) ! grain volume (not fraction but absolute)
case ('orientation')
mySize = 4_pInt
@ -3209,28 +3208,28 @@ function crystallite_postResults(&
case ('defgrad','f')
mySize = 9_pInt
crystallite_postResults(c+1:c+mySize) = reshape(math_transpose3x3(crystallite_partionedF(1:3,1:3,g,i,e)),(/mySize/))
crystallite_postResults(c+1:c+mySize) = reshape(math_transpose33(crystallite_partionedF(1:3,1:3,g,i,e)),(/mySize/))
case ('e')
mySize = 9_pInt
crystallite_postResults(c+1:c+mySize) = 0.5_pReal * reshape((math_mul33x33( &
math_transpose3x3(crystallite_partionedF(1:3,1:3,g,i,e)), &
math_transpose33(crystallite_partionedF(1:3,1:3,g,i,e)), &
crystallite_partionedF(1:3,1:3,g,i,e)) - math_I3),(/mySize/))
case ('fe')
mySize = 9_pInt
crystallite_postResults(c+1:c+mySize) = reshape(math_transpose3x3(crystallite_Fe(1:3,1:3,g,i,e)),(/mySize/))
crystallite_postResults(c+1:c+mySize) = reshape(math_transpose33(crystallite_Fe(1:3,1:3,g,i,e)),(/mySize/))
case ('ee')
Ee = 0.5_pReal * (math_mul33x33(math_transpose3x3(crystallite_Fe(1:3,1:3,g,i,e)), crystallite_Fe(1:3,1:3,g,i,e)) - math_I3)
Ee = 0.5_pReal * (math_mul33x33(math_transpose33(crystallite_Fe(1:3,1:3,g,i,e)), crystallite_Fe(1:3,1:3,g,i,e)) - math_I3)
mySize = 9_pInt
crystallite_postResults(c+1:c+mySize) = reshape(Ee,(/mySize/))
case ('fp')
mySize = 9_pInt
crystallite_postResults(c+1:c+mySize) = reshape(math_transpose3x3(crystallite_Fp(1:3,1:3,g,i,e)),(/mySize/))
crystallite_postResults(c+1:c+mySize) = reshape(math_transpose33(crystallite_Fp(1:3,1:3,g,i,e)),(/mySize/))
case ('lp')
mySize = 9_pInt
crystallite_postResults(c+1:c+mySize) = reshape(math_transpose3x3(crystallite_Lp(1:3,1:3,g,i,e)),(/mySize/))
crystallite_postResults(c+1:c+mySize) = reshape(math_transpose33(crystallite_Lp(1:3,1:3,g,i,e)),(/mySize/))
case ('p','firstpiola','1stpiola')
mySize = 9_pInt
crystallite_postResults(c+1:c+mySize) = reshape(math_transpose3x3(crystallite_P(1:3,1:3,g,i,e)),(/mySize/))
crystallite_postResults(c+1:c+mySize) = reshape(math_transpose33(crystallite_P(1:3,1:3,g,i,e)),(/mySize/))
case ('s','tstar','secondpiola','2ndpiola')
mySize = 9_pInt
crystallite_postResults(c+1:c+mySize) = reshape(math_Mandel6to33(crystallite_Tstar_v(1:6,g,i,e)),(/mySize/))

14
code/homogenization_RGC.f90
View File

@ -922,7 +922,7 @@ subroutine homogenization_RGC_stressPenalty(&
use prec, only: pReal,pInt,p_vec
use mesh, only: mesh_element
use constitutive, only: constitutive_homogenizedC
use math, only: math_civita,math_invert3x3
use math, only: math_civita,math_invert33
use material, only: homogenization_maxNgrains,homogenization_Ngrains
use numerics, only: xSmoo_RGC
implicit none
@ -1054,7 +1054,7 @@ subroutine homogenization_RGC_volumePenalty(&
)
use prec, only: pReal,pInt,p_vec
use mesh, only: mesh_element
use math, only: math_det3x3,math_inv3x3
use math, only: math_det33,math_inv33
use material, only: homogenization_maxNgrains,homogenization_Ngrains
use numerics, only: maxVolDiscr_RGC,volDiscrMod_RGC,volDiscrPow_RGC
@ -1072,9 +1072,9 @@ subroutine homogenization_RGC_volumePenalty(&
nGrain = homogenization_Ngrains(mesh_element(3,el))
!* Compute the volumes of grains and of cluster
vDiscrep = math_det3x3(fAvg) ! compute the volume of the cluster
vDiscrep = math_det33(fAvg) ! compute the volume of the cluster
do iGrain = 1,nGrain
gVol(iGrain) = math_det3x3(fDef(:,:,iGrain)) ! compute the volume of individual grains
gVol(iGrain) = math_det33(fDef(:,:,iGrain)) ! compute the volume of individual grains
vDiscrep = vDiscrep - gVol(iGrain)/dble(nGrain) ! calculate the difference/dicrepancy between
! the volume of the cluster and the the total volume of grains
enddo
@ -1084,7 +1084,7 @@ subroutine homogenization_RGC_volumePenalty(&
do iGrain = 1,nGrain
vPen(:,:,iGrain) = -1.0_pReal/dble(nGrain)*volDiscrMod_RGC*volDiscrPow_RGC/maxVolDiscr_RGC* &
sign((abs(vDiscrep)/maxVolDiscr_RGC)**(volDiscrPow_RGC - 1.0),vDiscrep)* &
gVol(iGrain)*transpose(math_inv3x3(fDef(:,:,iGrain)))
gVol(iGrain)*transpose(math_inv33(fDef(:,:,iGrain)))
!* Debugging the stress-like penalty of volume discrepancy
! if (ip == 1 .and. el == 1) then
@ -1108,7 +1108,7 @@ function homogenization_RGC_surfaceCorrection(&
)
use prec, only: pReal,pInt,p_vec
use math, only: math_invert3x3,math_mul33x33
use math, only: math_invert33,math_mul33x33
implicit none
!* Definition of variables
@ -1126,7 +1126,7 @@ function homogenization_RGC_surfaceCorrection(&
avgC = 0.0_pReal
avgC = math_mul33x33(transpose(avgF),avgF)
invC = 0.0_pReal
call math_invert3x3(avgC,invC,detF,error)
call math_invert33(avgC,invC,detF,error)
homogenization_RGC_surfaceCorrection = 0.0_pReal
do iBase = 1,3
intFace = (/iBase,1_pInt,1_pInt,1_pInt/)

4
code/lattice.f90
View File

@ -927,14 +927,14 @@ function lattice_initializeStructure(struct,CoverA)
lattice_st(1:3,i,myStructure) = st(1:3,i)
lattice_sn(1:3,i,myStructure) = sn(1:3,i)
lattice_Sslip(1:3,1:3,i,myStructure) = math_tensorproduct(sd(1:3,i),sn(1:3,i))
lattice_Sslip_v(1:6,i,myStructure) = math_Mandel33to6(math_symmetric3x3(lattice_Sslip(1:3,1:3,i,myStructure)))
lattice_Sslip_v(1:6,i,myStructure) = math_Mandel33to6(math_symmetric33(lattice_Sslip(1:3,1:3,i,myStructure)))
enddo
do i = 1,myNtwin ! store twin system vectors and Schmid plus rotation matrix for my structure
lattice_td(1:3,i,myStructure) = td(1:3,i)
lattice_tt(1:3,i,myStructure) = tt(1:3,i)
lattice_tn(1:3,i,myStructure) = tn(1:3,i)
lattice_Stwin(1:3,1:3,i,myStructure) = math_tensorproduct(td(1:3,i),tn(1:3,i))
lattice_Stwin_v(1:6,i,myStructure) = math_Mandel33to6(math_symmetric3x3(lattice_Stwin(1:3,1:3,i,myStructure)))
lattice_Stwin_v(1:6,i,myStructure) = math_Mandel33to6(math_symmetric33(lattice_Stwin(1:3,1:3,i,myStructure)))
lattice_Qtwin(1:3,1:3,i,myStructure) = math_AxisAngleToR(tn(1:3,i),180.0_pReal*inRad)
lattice_shearTwin(i,myStructure) = ts(i)
enddo

207
code/math.f90
View File

@ -495,7 +495,7 @@ real(pReal), dimension(4,36), parameter :: math_symOperations = &
!**************************************************************************
! matrix multiplication 33x33 = 3x3
! matrix multiplication 33x33 = 33
!**************************************************************************
pure function math_mul33x33(A,B)
@ -512,7 +512,7 @@ real(pReal), dimension(4,36), parameter :: math_symOperations = &
!**************************************************************************
! matrix multiplication 66x66 = 6x6
! matrix multiplication 66x66 = 66
!**************************************************************************
pure function math_mul66x66(A,B)
@ -530,7 +530,7 @@ real(pReal), dimension(4,36), parameter :: math_symOperations = &
!**************************************************************************
! matrix multiplication 99x99 = 9x9
! matrix multiplication 99x99 = 99
!**************************************************************************
pure function math_mul99x99(A,B)
@ -735,25 +735,25 @@ real(pReal), dimension(4,36), parameter :: math_symOperations = &
!**************************************************************************
! transposition of a 3x3 matrix
! transposition of a 33 matrix
!**************************************************************************
pure function math_transpose3x3(A)
pure function math_transpose33(A)
implicit none
real(pReal),dimension(3,3),intent(in) :: A
real(pReal),dimension(3,3) :: math_transpose3x3
real(pReal),dimension(3,3) :: math_transpose33
integer(pInt) :: i,j
forall(i=1_pInt:3_pInt, j=1_pInt:3_pInt) math_transpose3x3(i,j) = A(j,i)
forall(i=1_pInt:3_pInt, j=1_pInt:3_pInt) math_transpose33(i,j) = A(j,i)
endfunction math_transpose3x3
endfunction math_transpose33
!**************************************************************************
! Cramer inversion of 3x3 matrix (function)
! Cramer inversion of 33 matrix (function)
!**************************************************************************
pure function math_inv3x3(A)
pure function math_inv33(A)
! direct Cramer inversion of matrix A.
! returns all zeroes if not possible, i.e. if det close to zero
@ -762,37 +762,37 @@ pure function math_transpose3x3(A)
real(pReal),dimension(3,3),intent(in) :: A
real(pReal) :: DetA
real(pReal),dimension(3,3) :: math_inv3x3
real(pReal),dimension(3,3) :: math_inv33
math_inv3x3 = 0.0_pReal
math_inv33 = 0.0_pReal
DetA = A(1,1) * (A(2,2) * A(3,3) - A(2,3) * A(3,2))&
- A(1,2) * (A(2,1) * A(3,3) - A(2,3) * A(3,1))&
+ A(1,3) * (A(2,1) * A(3,2) - A(2,2) * A(3,1))
if (abs(DetA) > tiny(abs(DetA))) then
math_inv3x3(1,1) = ( A(2,2) * A(3,3) - A(2,3) * A(3,2)) / DetA
math_inv3x3(2,1) = (-A(2,1) * A(3,3) + A(2,3) * A(3,1)) / DetA
math_inv3x3(3,1) = ( A(2,1) * A(3,2) - A(2,2) * A(3,1)) / DetA
math_inv33(1,1) = ( A(2,2) * A(3,3) - A(2,3) * A(3,2)) / DetA
math_inv33(2,1) = (-A(2,1) * A(3,3) + A(2,3) * A(3,1)) / DetA
math_inv33(3,1) = ( A(2,1) * A(3,2) - A(2,2) * A(3,1)) / DetA
math_inv3x3(1,2) = (-A(1,2) * A(3,3) + A(1,3) * A(3,2)) / DetA
math_inv3x3(2,2) = ( A(1,1) * A(3,3) - A(1,3) * A(3,1)) / DetA
math_inv3x3(3,2) = (-A(1,1) * A(3,2) + A(1,2) * A(3,1)) / DetA
math_inv33(1,2) = (-A(1,2) * A(3,3) + A(1,3) * A(3,2)) / DetA
math_inv33(2,2) = ( A(1,1) * A(3,3) - A(1,3) * A(3,1)) / DetA
math_inv33(3,2) = (-A(1,1) * A(3,2) + A(1,2) * A(3,1)) / DetA
math_inv3x3(1,3) = ( A(1,2) * A(2,3) - A(1,3) * A(2,2)) / DetA
math_inv3x3(2,3) = (-A(1,1) * A(2,3) + A(1,3) * A(2,1)) / DetA
math_inv3x3(3,3) = ( A(1,1) * A(2,2) - A(1,2) * A(2,1)) / DetA
math_inv33(1,3) = ( A(1,2) * A(2,3) - A(1,3) * A(2,2)) / DetA
math_inv33(2,3) = (-A(1,1) * A(2,3) + A(1,3) * A(2,1)) / DetA
math_inv33(3,3) = ( A(1,1) * A(2,2) - A(1,2) * A(2,1)) / DetA
endif
endfunction math_inv3x3
endfunction math_inv33
!**************************************************************************
! Cramer inversion of 3x3 matrix (subroutine)
! Cramer inversion of 33 matrix (subroutine)
!**************************************************************************
PURE SUBROUTINE math_invert3x3(A, InvA, DetA, error)
PURE SUBROUTINE math_invert33(A, InvA, DetA, error)
! Bestimmung der Determinanten und Inversen einer 3x3-Matrix
! Bestimmung der Determinanten und Inversen einer 33-Matrix
! A = Matrix A
! InvA = Inverse of A
! DetA = Determinant of A
@ -827,7 +827,7 @@ pure function math_transpose3x3(A)
error = .false.
endif
ENDSUBROUTINE math_invert3x3
ENDSUBROUTINE math_invert33
!**************************************************************************
@ -1018,51 +1018,72 @@ pure function math_transpose3x3(A)
!********************************************************************
! symmetrize a 3x3 matrix
! symmetrize a 33 matrix
!********************************************************************
function math_symmetric3x3(m)
function math_symmetric33(m)
implicit none
real(pReal), dimension(3,3) :: math_symmetric3x3
real(pReal), dimension(3,3) :: math_symmetric33
real(pReal), dimension(3,3), intent(in) :: m
integer(pInt) :: i,j
forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt) math_symmetric3x3(i,j) = 0.5_pReal * (m(i,j) + m(j,i))
forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt) math_symmetric33(i,j) = 0.5_pReal * (m(i,j) + m(j,i))
endfunction math_symmetric3x3
endfunction math_symmetric33
!********************************************************************
! symmetrize a 6x6 matrix
! symmetrize a 66 matrix
!********************************************************************
pure function math_symmetric6x6(m)
pure function math_symmetric66(m)
implicit none
integer(pInt) :: i,j
real(pReal), dimension(6,6), intent(in) :: m
real(pReal), dimension(6,6) :: math_symmetric6x6
real(pReal), dimension(6,6) :: math_symmetric66
forall (i=1_pInt:6_pInt,j=1_pInt:6_pInt) math_symmetric6x6(i,j) = 0.5_pReal * (m(i,j) + m(j,i))
forall (i=1_pInt:6_pInt,j=1_pInt:6_pInt) math_symmetric66(i,j) = 0.5_pReal * (m(i,j) + m(j,i))
endfunction math_symmetric66
endfunction math_symmetric6x6
!********************************************************************
! skew part of a 3x3 matrix
! skew part of a 33 matrix
!********************************************************************
function math_skew3x3(m)
function math_skew33(m)
implicit none
real(pReal), dimension(3,3) :: math_skew3x3
real(pReal), dimension(3,3) :: math_skew33
real(pReal), dimension(3,3), intent(in) :: m
integer(pInt) :: i,j
forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt) math_skew3x3(i,j) = m(i,j) - 0.5_pReal * (m(i,j) + m(j,i))
forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt) math_skew33(i,j) = m(i,j) - 0.5_pReal * (m(i,j) + m(j,i))
endfunction math_skew33
endfunction math_skew3x3
!********************************************************************
! deviatoric part of a 33 matrix
!********************************************************************
function math_deviatoric33(m)
implicit none
real(pReal), dimension(3,3) :: math_deviatoric33
real(pReal), dimension(3,3), intent(in) :: m
integer(pInt) :: i
real(pReal) :: hydrostatic
hydrostatic = (m(1,1) + m(2,2) + m(3,3)) / 3.0_pReal
math_deviatoric33 = m
forall (i=1_pInt:3_pInt) math_deviatoric33(i,i) = m(i,i) - hydrostatic
endfunction math_deviatoric33
!********************************************************************
! equivalent scalar quantity of a full strain tensor
!********************************************************************
@ -1119,24 +1140,24 @@ pure function math_transpose3x3(A)
!********************************************************************
! determinant of a 3x3 matrix
! determinant of a 33 matrix
!********************************************************************
pure function math_det3x3(m)
pure function math_det33(m)
implicit none
real(pReal), dimension(3,3), intent(in) :: m
real(pReal) :: math_det3x3
real(pReal) :: math_det33
math_det3x3 = m(1,1)*(m(2,2)*m(3,3)-m(2,3)*m(3,2)) &
math_det33 = m(1,1)*(m(2,2)*m(3,3)-m(2,3)*m(3,2)) &
-m(1,2)*(m(2,1)*m(3,3)-m(2,3)*m(3,1)) &
+m(1,3)*(m(2,1)*m(3,2)-m(2,2)*m(3,1))
endfunction math_det3x3
endfunction math_det33
!********************************************************************
! norm of a 3x3 matrix
! norm of a 33 matrix
!********************************************************************
pure function math_norm33(m)
@ -1151,7 +1172,7 @@ pure function math_transpose3x3(A)
!********************************************************************
! euclidic norm of a 3x1 vector
! euclidic norm of a 3 vector
!********************************************************************
pure function math_norm3(v)
@ -1166,7 +1187,7 @@ pure function math_transpose3x3(A)
!********************************************************************
! convert 3x3 matrix into vector 9x1
! convert 33 matrix into vector 9
!********************************************************************
pure function math_Plain33to9(m33)
@ -1182,7 +1203,7 @@ pure function math_transpose3x3(A)
!********************************************************************
! convert Plain 9x1 back to 3x3 matrix
! convert Plain 9 back to 33 matrix
!********************************************************************
pure function math_Plain9to33(v9)
@ -1198,7 +1219,7 @@ pure function math_transpose3x3(A)
!********************************************************************
! convert symmetric 3x3 matrix into Mandel vector 6x1
! convert symmetric 33 matrix into Mandel vector 6
!********************************************************************
pure function math_Mandel33to6(m33)
@ -1214,7 +1235,7 @@ pure function math_transpose3x3(A)
!********************************************************************
! convert Mandel 6x1 back to symmetric 3x3 matrix
! convert Mandel 6 back to symmetric 33 matrix
!********************************************************************
pure function math_Mandel6to33(v6)
@ -1233,7 +1254,7 @@ pure function math_transpose3x3(A)
!********************************************************************
! convert 3x3x3x3 tensor into plain matrix 9x9
! convert 3333 tensor into plain matrix 99
!********************************************************************
pure function math_Plain3333to99(m3333)
@ -1249,7 +1270,7 @@ pure function math_transpose3x3(A)
endfunction math_Plain3333to99
!********************************************************************
! plain matrix 9x9 into 3x3x3x3 tensor
! plain matrix 99 into 3333 tensor
!********************************************************************
pure function math_Plain99to3333(m99)
@ -1266,7 +1287,7 @@ pure function math_transpose3x3(A)
!********************************************************************
! convert Mandel matrix 6x6 into Plain matrix 6x6
! convert Mandel matrix 66 into Plain matrix 66
!********************************************************************
pure function math_Mandel66toPlain66(m66)
@ -1284,7 +1305,7 @@ pure function math_transpose3x3(A)
!********************************************************************
! convert Plain matrix 6x6 into Mandel matrix 6x6
! convert Plain matrix 66 into Mandel matrix 66
!********************************************************************
pure function math_Plain66toMandel66(m66)
@ -1302,7 +1323,7 @@ pure function math_transpose3x3(A)
!********************************************************************
! convert symmetric 3x3x3x3 tensor into Mandel matrix 6x6
! convert symmetric 3333 tensor into Mandel matrix 66
!********************************************************************
pure function math_Mandel3333to66(m3333)
@ -1319,7 +1340,7 @@ pure function math_transpose3x3(A)
!********************************************************************
! convert Mandel matrix 6x6 back to symmetric 3x3x3x3 tensor
! convert Mandel matrix 66 back to symmetric 3333 tensor
!********************************************************************
pure function math_Mandel66to3333(m66)
@ -1340,7 +1361,7 @@ pure function math_transpose3x3(A)
!********************************************************************
! convert Voigt matrix 6x6 back to symmetric 3x3x3x3 tensor
! convert Voigt matrix 66 back to symmetric 3333 tensor
!********************************************************************
pure function math_Voigt66to3333(m66)
@ -2015,7 +2036,7 @@ endfunction math_sampleGaussVar
!****************************************************************
subroutine math_spectralDecompositionSym3x3(M,values,vectors,error)
subroutine math_spectralDecompositionSym33(M,values,vectors,error)
!****************************************************************
implicit none
@ -2048,11 +2069,11 @@ end subroutine
real(pReal) :: EW1, EW2, EW3, det
error = .false.
ce = math_mul33x33(math_transpose3x3(FE),FE)
ce = math_mul33x33(math_transpose33(FE),FE)
CALL math_spectral1(CE,EW1,EW2,EW3,EB1,EB2,EB3)
U=sqrt(EW1)*EB1+sqrt(EW2)*EB2+sqrt(EW3)*EB3
call math_invert3x3(U,UI,det,error)
call math_invert33(U,UI,det,error)
if (.not. error) R = math_mul33x33(FE,UI)
ENDSUBROUTINE math_pDecomposition
@ -2157,14 +2178,14 @@ end subroutine
!**********************************************************************
function math_eigenvalues3x3(M)
function math_eigenvalues33(M)
!**** Eigenvalues of symmetric 3X3 matrix M
implicit none
real(pReal), intent(in), dimension(3,3) :: M
real(pReal), dimension(3,3) :: EB1 = 0.0_pReal, EB2 = 0.0_pReal, EB3 = 0.0_pReal
real(pReal), dimension(3) :: math_eigenvalues3x3
real(pReal), dimension(3) :: math_eigenvalues33
real(pReal) :: HI1M, HI2M, HI3M, R, S, T, P, Q, RHO, PHI, Y1, Y2, Y3, arg
real(pReal), parameter :: TOL=1.e-14_pReal
@ -2177,9 +2198,9 @@ end subroutine
if((abs(P) < TOL) .and. (abs(Q) < TOL)) THEN
! three equivalent eigenvalues
math_eigenvalues3x3(1) = HI1M/3.0_pReal
math_eigenvalues3x3(2)=math_eigenvalues3x3(1)
math_eigenvalues3x3(3)=math_eigenvalues3x3(1)
math_eigenvalues33(1) = HI1M/3.0_pReal
math_eigenvalues33(2)=math_eigenvalues33(1)
math_eigenvalues33(3)=math_eigenvalues33(1)
! this is not really correct, but this way U is calculated
! correctly in PDECOMPOSITION (correct is EB?=I)
EB1(1,1)=1.0_pReal
@ -2194,11 +2215,11 @@ end subroutine
Y1=2*RHO**(1.0_pReal/3.0_pReal)*cos(PHI/3.0_pReal)
Y2=2*RHO**(1.0_pReal/3.0_pReal)*cos(PHI/3.0_pReal+2.0_pReal/3.0_pReal*PI)
Y3=2*RHO**(1.0_pReal/3.0_pReal)*cos(PHI/3.0_pReal+4.0_pReal/3.0_pReal*PI)
math_eigenvalues3x3(1) = Y1-R/3.0_pReal
math_eigenvalues3x3(2) = Y2-R/3.0_pReal
math_eigenvalues3x3(3) = Y3-R/3.0_pReal
math_eigenvalues33(1) = Y1-R/3.0_pReal
math_eigenvalues33(2) = Y2-R/3.0_pReal
math_eigenvalues33(3) = Y3-R/3.0_pReal
endif
endfunction math_eigenvalues3x3
endfunction math_eigenvalues33
!**********************************************************************
@ -2214,7 +2235,7 @@ end subroutine
HI1M=M(1,1)+M(2,2)+M(3,3)
HI2M=HI1M**2.0_pReal/2.0_pReal- (M(1,1)**2.0_pReal+M(2,2)**2.0_pReal+M(3,3)**2.0_pReal)&
/2.0_pReal-M(1,2)*M(2,1)-M(1,3)*M(3,1)-M(2,3)*M(3,2)
HI3M=math_det3x3(M)
HI3M=math_det33(M)
! QUESTION: is 3rd equiv det(M) ?? if yes, use function math_det !agreed on YES
ENDSUBROUTINE math_hi
@ -2810,65 +2831,65 @@ end subroutine
m(1:3,2) = v2-v3
m(1:3,3) = v3-v4
math_volTetrahedron = math_det3x3(m)/6.0_pReal
math_volTetrahedron = math_det33(m)/6.0_pReal
endfunction math_volTetrahedron
!**************************************************************************
! rotate 3x3 tensor forward
! rotate 33 tensor forward
!**************************************************************************
pure function math_rotate_forward3x3(tensor,rot_tensor)
pure function math_rotate_forward33(tensor,rot_tensor)
implicit none
real(pReal), dimension(3,3) :: math_rotate_forward3x3
real(pReal), dimension(3,3) :: math_rotate_forward33
real(pReal), dimension(3,3), intent(in) :: tensor, rot_tensor
math_rotate_forward3x3 = math_mul33x33(rot_tensor,&
math_mul33x33(tensor,math_transpose3x3(rot_tensor)))
math_rotate_forward33 = math_mul33x33(rot_tensor,&
math_mul33x33(tensor,math_transpose33(rot_tensor)))
endfunction math_rotate_forward3x3
endfunction math_rotate_forward33
!**************************************************************************
! rotate 3x3 tensor backward
! rotate 33 tensor backward
!**************************************************************************
pure function math_rotate_backward3x3(tensor,rot_tensor)
pure function math_rotate_backward33(tensor,rot_tensor)
implicit none
real(pReal), dimension(3,3) :: math_rotate_backward3x3
real(pReal), dimension(3,3) :: math_rotate_backward33
real(pReal), dimension(3,3), intent(in) :: tensor, rot_tensor
math_rotate_backward3x3 = math_mul33x33(math_transpose3x3(rot_tensor),&
math_rotate_backward33 = math_mul33x33(math_transpose33(rot_tensor),&
math_mul33x33(tensor,rot_tensor))
endfunction math_rotate_backward3x3
endfunction math_rotate_backward33
!**************************************************************************
! rotate 3x3x3x3 tensor
! rotate 3333 tensor
! C'_ijkl=g_im*g_jn*g_ko*g_lp*C_mnop
!**************************************************************************
pure function math_rotate_forward3x3x3x3(tensor,rot_tensor)
pure function math_rotate_forward3333(tensor,rot_tensor)
implicit none
real(pReal), dimension(3,3,3,3) :: math_rotate_forward3x3x3x3
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), intent(in) :: tensor
integer(pInt) :: i,j,k,l,m,n,o,p
math_rotate_forward3x3x3x3= 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 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_forward3x3x3x3(i,j,k,l) = tensor(i,j,k,l)+rot_tensor(m,i)*rot_tensor(n,j)*&
math_rotate_forward3333(i,j,k,l) = tensor(i,j,k,l)+rot_tensor(m,i)*rot_tensor(n,j)*&
rot_tensor(o,k)*rot_tensor(p,l)*tensor(m,n,o,p)
enddo; enddo; enddo; enddo; enddo; enddo; enddo; enddo
endfunction math_rotate_forward3x3x3x3
endfunction math_rotate_forward3333
!+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
@ -2950,7 +2971,7 @@ end subroutine
+ abs(math_volTetrahedron(coords(7,1:3),coords(1,1:3),coords(3,1:3),coords(2,1:3))) &
+ abs(math_volTetrahedron(coords(7,1:3),coords(5,1:3),coords(2,1:3),coords(6,1:3))) &
+ abs(math_volTetrahedron(coords(7,1:3),coords(5,1:3),coords(2,1:3),coords(1,1:3)))
volume_mismatch(i,j,k) = volume_mismatch(i,j,k)/math_det3x3(defgrad(i,j,k,1:3,1:3))
volume_mismatch(i,j,k) = volume_mismatch(i,j,k)/math_det33(defgrad(i,j,k,1:3,1:3))
enddo; enddo; enddo
volume_mismatch = volume_mismatch/vol_initial
@ -3630,7 +3651,7 @@ subroutine logstrain_spat(res,defgrad,logstrain_field)
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
call math_pDecomposition(defgrad(i,j,k,1:3,1:3),temp33_Real2,temp33_Real,errmatinv) !store R in temp33_Real
temp33_Real2 = math_inv3x3(temp33_Real)
temp33_Real2 = math_inv33(temp33_Real)
temp33_Real = math_mul33x33(defgrad(i,j,k,1:3,1:3),temp33_Real2) ! v = F o inv(R), store in temp33_Real2
call math_spectral1(temp33_Real,eigenvalue(1), eigenvalue(2), eigenvalue(3),&
eigenvectorbasis(1,1:3,1:3),eigenvectorbasis(2,1:3,1:3),eigenvectorbasis(3,1:3,1:3))
@ -3690,7 +3711,7 @@ subroutine calculate_cauchy(res,defgrad,p_stress,c_stress)
c_stress = 0.0_pInt
do k = 1_pInt, res(3); do j = 1_pInt, res(2); do i = 1_pInt, res(1)
jacobi = math_det3x3(defgrad(i,j,k,1:3,1:3))
jacobi = math_det33(defgrad(i,j,k,1:3,1:3))
c_stress(i,j,k,1:3,1:3) = matmul(p_stress(i,j,k,1:3,1:3),transpose(defgrad(i,j,k,1:3,1:3)))/jacobi
enddo; enddo; enddo