diff --git a/code/CPFEM.f90 b/code/CPFEM.f90 index 684112263..3db788dee 100644 --- a/code/CPFEM.f90 +++ b/code/CPFEM.f90 @@ -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 diff --git a/code/DAMASK_spectral.f90 b/code/DAMASK_spectral.f90 index 5eb9f4c07..f8771df77 100644 --- a/code/DAMASK_spectral.f90 +++ b/code/DAMASK_spectral.f90 @@ -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 diff --git a/code/constitutive_dislotwin.f90 b/code/constitutive_dislotwin.f90 index 2eb1a11e1..4a88346ce 100644 --- a/code/constitutive_dislotwin.f90 +++ b/code/constitutive_dislotwin.f90 @@ -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)) diff --git a/code/constitutive_j2.f90 b/code/constitutive_j2.f90 index 44dc9963f..63d801079 100644 --- a/code/constitutive_j2.f90 +++ b/code/constitutive_j2.f90 @@ -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, & diff --git a/code/constitutive_nonlocal.f90 b/code/constitutive_nonlocal.f90 index d56406c20..2de72e4b9 100644 --- a/code/constitutive_nonlocal.f90 +++ b/code/constitutive_nonlocal.f90 @@ -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, & diff --git a/code/crystallite.f90 b/code/crystallite.f90 index 2cf8a958b..749a5a9a7 100644 --- a/code/crystallite.f90 +++ b/code/crystallite.f90 @@ -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/)) diff --git a/code/homogenization_RGC.f90 b/code/homogenization_RGC.f90 index d03f18a73..374a34bcb 100644 --- a/code/homogenization_RGC.f90 +++ b/code/homogenization_RGC.f90 @@ -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/) diff --git a/code/lattice.f90 b/code/lattice.f90 index f3009f070..88861a87c 100644 --- a/code/lattice.f90 +++ b/code/lattice.f90 @@ -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 diff --git a/code/math.f90 b/code/math.f90 index e5d93d76a..27038e120 100644 --- a/code/math.f90 +++ b/code/math.f90 @@ -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