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

use matmul instead of hand-written functions 

- performance is the same
- leaner code
- matmul works (was buggy a few years ago)
This commit is contained in:
Martin Diehl 2019-04-03 08:22:04 +02:00
parent 0cf2c7b9e6
commit 4604e65a42
17 changed files with 123 additions and 173 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
src/CPFEM.f90
View File

@ -259,7 +259,6 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
restartWrite restartWrite
use math, only: & use math, only: &
math_identity2nd, & math_identity2nd, &
math_mul33x33, &
math_det33, & math_det33, &
math_delta, & math_delta, &
math_sym3333to66, & math_sym3333to66, &
@ -557,7 +556,7 @@ subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt
endif endif
! translate from P to CS ! translate from P to CS
Kirchhoff = math_mul33x33(materialpoint_P(1:3,1:3,ip,elCP), transpose(materialpoint_F(1:3,1:3,ip,elCP))) Kirchhoff = matmul(materialpoint_P(1:3,1:3,ip,elCP), transpose(materialpoint_F(1:3,1:3,ip,elCP)))
J_inverse = 1.0_pReal / math_det33(materialpoint_F(1:3,1:3,ip,elCP)) J_inverse = 1.0_pReal / math_det33(materialpoint_F(1:3,1:3,ip,elCP))
CPFEM_cs(1:6,ip,elCP) = math_sym33to6(J_inverse * Kirchhoff,weighted=.false.) CPFEM_cs(1:6,ip,elCP) = math_sym33to6(J_inverse * Kirchhoff,weighted=.false.)

2
src/DAMASK_grid.f90
View File

@ -337,7 +337,7 @@ program DAMASK_spectral
endif endif
enddo; write(6,'(/)',advance='no') enddo; write(6,'(/)',advance='no')
enddo enddo
if (any(abs(math_mul33x33(newLoadCase%rotation, & if (any(abs(matmul(newLoadCase%rotation, &
transpose(newLoadCase%rotation))-math_I3) > & transpose(newLoadCase%rotation))-math_I3) > &
reshape(spread(tol_math_check,1,9),[ 3,3]))& reshape(spread(tol_math_check,1,9),[ 3,3]))&
.or. abs(math_det33(newLoadCase%rotation)) > & .or. abs(math_det33(newLoadCase%rotation)) > &

41
src/constitutive.f90
View File

@ -381,8 +381,6 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
S, Fi, ipc, ip, el) S, Fi, ipc, ip, el)
use prec, only: & use prec, only: &
pReal pReal
use math, only: &
math_mul33x33
use material, only: & use material, only: &
phasememberAt, & phasememberAt, &
phase_plasticity, & phase_plasticity, &
@ -439,7 +437,7 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
ho = material_homogenizationAt(el) ho = material_homogenizationAt(el)
tme = thermalMapping(ho)%p(ip,el) tme = thermalMapping(ho)%p(ip,el)
Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),S) Mp = matmul(matmul(transpose(Fi),Fi),S)
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
@ -483,9 +481,9 @@ subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
#else #else
do concurrent(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt) do concurrent(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt)
#endif #endif
dLp_dFi(i,j,1:3,1:3) = math_mul33x33(math_mul33x33(Fi,S),transpose(dLp_dMp(i,j,1:3,1:3))) + & dLp_dFi(i,j,1:3,1:3) = matmul(matmul(Fi,S),transpose(dLp_dMp(i,j,1:3,1:3))) + &
math_mul33x33(math_mul33x33(Fi,dLp_dMp(i,j,1:3,1:3)),S) matmul(matmul(Fi,dLp_dMp(i,j,1:3,1:3)),S)
dLp_dS(i,j,1:3,1:3) = math_mul33x33(math_mul33x33(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3)) ! ToDo: @PS: why not: dLp_dMp:(FiT Fi) dLp_dS(i,j,1:3,1:3) = matmul(matmul(transpose(Fi),Fi),dLp_dMp(i,j,1:3,1:3)) ! ToDo: @PS: why not: dLp_dMp:(FiT Fi)
#if defined(__INTEL_COMPILER) || defined(__PGI) #if defined(__INTEL_COMPILER) || defined(__PGI)
end forall end forall
#else #else
@ -506,8 +504,7 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
use math, only: & use math, only: &
math_I3, & math_I3, &
math_inv33, & math_inv33, &
math_det33, & math_det33
math_mul33x33
use material, only: & use material, only: &
phasememberAt, & phasememberAt, &
phase_plasticity, & phase_plasticity, &
@ -591,11 +588,11 @@ subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
FiInv = math_inv33(Fi) FiInv = math_inv33(Fi)
detFi = math_det33(Fi) detFi = math_det33(Fi)
Li = math_mul33x33(math_mul33x33(Fi,Li),FiInv)*detFi !< push forward to intermediate configuration Li = matmul(matmul(Fi,Li),FiInv)*detFi !< push forward to intermediate configuration
temp_33 = math_mul33x33(FiInv,Li) temp_33 = matmul(FiInv,Li)
do i = 1_pInt,3_pInt; do j = 1_pInt,3_pInt do i = 1_pInt,3_pInt; do j = 1_pInt,3_pInt
dLi_dS(1:3,1:3,i,j) = math_mul33x33(math_mul33x33(Fi,dLi_dS(1:3,1:3,i,j)),FiInv)*detFi dLi_dS(1:3,1:3,i,j) = matmul(matmul(Fi,dLi_dS(1:3,1:3,i,j)),FiInv)*detFi
dLi_dFi(1:3,1:3,i,j) = dLi_dFi(1:3,1:3,i,j) + Li*FiInv(j,i) dLi_dFi(1:3,1:3,i,j) = dLi_dFi(1:3,1:3,i,j) + Li*FiInv(j,i)
dLi_dFi(1:3,i,1:3,j) = dLi_dFi(1:3,i,1:3,j) + math_I3*temp_33(j,i) + Li*FiInv(j,i) dLi_dFi(1:3,i,1:3,j) = dLi_dFi(1:3,i,1:3,j) + math_I3*temp_33(j,i) + Li*FiInv(j,i)
end do; end do end do; end do
@ -689,7 +686,6 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
use prec, only: & use prec, only: &
pReal pReal
use math, only : & use math, only : &
math_mul33x33, &
math_mul3333xx33, & math_mul3333xx33, &
math_66toSym3333, & math_66toSym3333, &
math_I3 math_I3
@ -733,14 +729,14 @@ subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
end select degradationType end select degradationType
enddo DegradationLoop enddo DegradationLoop
E = 0.5_pReal*(math_mul33x33(transpose(Fe),Fe)-math_I3) !< Green-Lagrange strain in unloaded configuration E = 0.5_pReal*(matmul(transpose(Fe),Fe)-math_I3) !< Green-Lagrange strain in unloaded configuration
S = math_mul3333xx33(C,math_mul33x33(math_mul33x33(transpose(Fi),E),Fi)) !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration S = math_mul3333xx33(C,matmul(matmul(transpose(Fi),E),Fi)) !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration
dS_dFe = 0.0_pReal dS_dFe = 0.0_pReal
forall (i=1_pInt:3_pInt, j=1_pInt:3_pInt) forall (i=1_pInt:3_pInt, j=1_pInt:3_pInt)
dS_dFe(i,j,1:3,1:3) = & dS_dFe(i,j,1:3,1:3) = &
math_mul33x33(Fe,math_mul33x33(math_mul33x33(Fi,C(i,j,1:3,1:3)),transpose(Fi))) !< dS_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko matmul(Fe,matmul(matmul(Fi,C(i,j,1:3,1:3)),transpose(Fi))) !< dS_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko
dS_dFi(i,j,1:3,1:3) = 2.0_pReal*math_mul33x33(math_mul33x33(E,Fi),C(i,j,1:3,1:3)) !< dS_ij/dFi_kl = C_ijln * E_km * Fe_mn dS_dFi(i,j,1:3,1:3) = 2.0_pReal*matmul(matmul(E,Fi),C(i,j,1:3,1:3)) !< dS_ij/dFi_kl = C_ijln * E_km * Fe_mn
end forall end forall
end subroutine constitutive_hooke_SandItsTangents end subroutine constitutive_hooke_SandItsTangents
@ -756,9 +752,6 @@ subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, ipc, ip,
debug_level, & debug_level, &
debug_constitutive, & debug_constitutive, &
debug_levelBasic debug_levelBasic
use math, only: &
math_mul33x33, &
math_mul33x33
use mesh, only: & use mesh, only: &
theMesh theMesh
use material, only: & use material, only: &
@ -829,7 +822,7 @@ subroutine constitutive_collectDotState(S, FeArray, Fi, FpArray, subdt, ipc, ip,
ho = material_homogenizationAt(el) ho = material_homogenizationAt(el)
tme = thermalMapping(ho)%p(ip,el) tme = thermalMapping(ho)%p(ip,el)
Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),S) Mp = matmul(matmul(transpose(Fi),Fi),S)
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
@ -897,8 +890,6 @@ subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el)
debug_level, & debug_level, &
debug_constitutive, & debug_constitutive, &
debug_levelBasic debug_levelBasic
use math, only: &
math_mul33x33
use material, only: & use material, only: &
phasememberAt, & phasememberAt, &
phase_plasticityInstance, & phase_plasticityInstance, &
@ -931,7 +922,7 @@ subroutine constitutive_collectDeltaState(S, Fe, Fi, ipc, ip, el)
i, & i, &
instance, of instance, of
Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),S) Mp = matmul(matmul(transpose(Fi),Fi),S)
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
@ -966,8 +957,6 @@ end subroutine constitutive_collectDeltaState
function constitutive_postResults(S, Fi, ipc, ip, el) function constitutive_postResults(S, Fi, ipc, ip, el)
use prec, only: & use prec, only: &
pReal pReal
use math, only: &
math_mul33x33
use material, only: & use material, only: &
phasememberAt, & phasememberAt, &
phase_plasticityInstance, & phase_plasticityInstance, &
@ -1035,7 +1024,7 @@ function constitutive_postResults(S, Fi, ipc, ip, el)
constitutive_postResults = 0.0_pReal constitutive_postResults = 0.0_pReal
Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),S) Mp = matmul(matmul(transpose(Fi),Fi),S)
ho = material_homogenizationAt(el) ho = material_homogenizationAt(el)
tme = thermalMapping(ho)%p(ip,el) tme = thermalMapping(ho)%p(ip,el)

81
src/crystallite.f90
View File

@ -144,8 +144,7 @@ subroutine crystallite_init
use math, only: & use math, only: &
math_I3, & math_I3, &
math_EulerToR, & math_EulerToR, &
math_inv33, & math_inv33
math_mul33x33
use mesh, only: & use mesh, only: &
theMesh, & theMesh, &
mesh_element mesh_element
@ -353,7 +352,7 @@ subroutine crystallite_init
crystallite_Fi0(1:3,1:3,c,i,e) = constitutive_initialFi(c,i,e) crystallite_Fi0(1:3,1:3,c,i,e) = constitutive_initialFi(c,i,e)
crystallite_F0(1:3,1:3,c,i,e) = math_I3 crystallite_F0(1:3,1:3,c,i,e) = math_I3
crystallite_localPlasticity(c,i,e) = phase_localPlasticity(material_phase(c,i,e)) crystallite_localPlasticity(c,i,e) = phase_localPlasticity(material_phase(c,i,e))
crystallite_Fe(1:3,1:3,c,i,e) = math_inv33(math_mul33x33(crystallite_Fi0(1:3,1:3,c,i,e), & crystallite_Fe(1:3,1:3,c,i,e) = math_inv33(matmul(crystallite_Fi0(1:3,1:3,c,i,e), &
crystallite_Fp0(1:3,1:3,c,i,e))) ! assuming that euler angles are given in internal strain free configuration crystallite_Fp0(1:3,1:3,c,i,e))) ! assuming that euler angles are given in internal strain free configuration
crystallite_Fp(1:3,1:3,c,i,e) = crystallite_Fp0(1:3,1:3,c,i,e) crystallite_Fp(1:3,1:3,c,i,e) = crystallite_Fp0(1:3,1:3,c,i,e)
crystallite_Fi(1:3,1:3,c,i,e) = crystallite_Fi0(1:3,1:3,c,i,e) crystallite_Fi(1:3,1:3,c,i,e) = crystallite_Fi0(1:3,1:3,c,i,e)
@ -430,8 +429,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
IO_warning, & IO_warning, &
IO_error IO_error
use math, only: & use math, only: &
math_inv33, & math_inv33
math_mul33x33
use mesh, only: & use mesh, only: &
theMesh, & theMesh, &
mesh_element mesh_element
@ -602,7 +600,7 @@ function crystallite_stress(dummyArgumentToPreventInternalCompilerErrorWithGCC)
crystallite_subF(1:3,1:3,c,i,e) = crystallite_subF0(1:3,1:3,c,i,e) & crystallite_subF(1:3,1:3,c,i,e) = crystallite_subF0(1:3,1:3,c,i,e) &
+ crystallite_subStep(c,i,e) * (crystallite_partionedF (1:3,1:3,c,i,e) & + crystallite_subStep(c,i,e) * (crystallite_partionedF (1:3,1:3,c,i,e) &
- crystallite_partionedF0(1:3,1:3,c,i,e)) - crystallite_partionedF0(1:3,1:3,c,i,e))
crystallite_Fe(1:3,1:3,c,i,e) = math_mul33x33(math_mul33x33(crystallite_subF (1:3,1:3,c,i,e), & crystallite_Fe(1:3,1:3,c,i,e) = matmul(matmul(crystallite_subF (1:3,1:3,c,i,e), &
crystallite_invFp(1:3,1:3,c,i,e)), & crystallite_invFp(1:3,1:3,c,i,e)), &
crystallite_invFi(1:3,1:3,c,i,e)) crystallite_invFi(1:3,1:3,c,i,e))
crystallite_subdt(c,i,e) = crystallite_subStep(c,i,e) * crystallite_dt(c,i,e) crystallite_subdt(c,i,e) = crystallite_subStep(c,i,e) * crystallite_dt(c,i,e)
@ -691,7 +689,6 @@ subroutine crystallite_stressTangent()
use math, only: & use math, only: &
math_inv33, & math_inv33, &
math_identity2nd, & math_identity2nd, &
math_mul33x33, &
math_3333to99, & math_3333to99, &
math_99to3333, & math_99to3333, &
math_I3, & math_I3, &
@ -753,11 +750,11 @@ subroutine crystallite_stressTangent()
lhs_3333 = 0.0_pReal; rhs_3333 = 0.0_pReal lhs_3333 = 0.0_pReal; rhs_3333 = 0.0_pReal
do o=1_pInt,3_pInt; do p=1_pInt,3_pInt do o=1_pInt,3_pInt; do p=1_pInt,3_pInt
lhs_3333(1:3,1:3,o,p) = lhs_3333(1:3,1:3,o,p) & lhs_3333(1:3,1:3,o,p) = lhs_3333(1:3,1:3,o,p) &
+ crystallite_subdt(c,i,e)*math_mul33x33(invSubFi0,dLidFi(1:3,1:3,o,p)) + crystallite_subdt(c,i,e)*matmul(invSubFi0,dLidFi(1:3,1:3,o,p))
lhs_3333(1:3,o,1:3,p) = lhs_3333(1:3,o,1:3,p) & lhs_3333(1:3,o,1:3,p) = lhs_3333(1:3,o,1:3,p) &
+ crystallite_invFi(1:3,1:3,c,i,e)*crystallite_invFi(p,o,c,i,e) + crystallite_invFi(1:3,1:3,c,i,e)*crystallite_invFi(p,o,c,i,e)
rhs_3333(1:3,1:3,o,p) = rhs_3333(1:3,1:3,o,p) & rhs_3333(1:3,1:3,o,p) = rhs_3333(1:3,1:3,o,p) &
- crystallite_subdt(c,i,e)*math_mul33x33(invSubFi0,dLidS(1:3,1:3,o,p)) - crystallite_subdt(c,i,e)*matmul(invSubFi0,dLidS(1:3,1:3,o,p))
enddo;enddo enddo;enddo
call math_invert2(temp_99,error,math_3333to99(lhs_3333)) call math_invert2(temp_99,error,math_3333to99(lhs_3333))
if (error) then if (error) then
@ -777,19 +774,19 @@ subroutine crystallite_stressTangent()
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculate dSdF ! calculate dSdF
temp_33_1 = transpose(math_mul33x33(crystallite_invFp(1:3,1:3,c,i,e), & temp_33_1 = transpose(matmul(crystallite_invFp(1:3,1:3,c,i,e), &
crystallite_invFi(1:3,1:3,c,i,e))) crystallite_invFi(1:3,1:3,c,i,e)))
temp_33_2 = math_mul33x33( crystallite_subF (1:3,1:3,c,i,e), & temp_33_2 = matmul( crystallite_subF (1:3,1:3,c,i,e), &
math_inv33(crystallite_subFp0(1:3,1:3,c,i,e))) math_inv33(crystallite_subFp0(1:3,1:3,c,i,e)))
temp_33_3 = math_mul33x33(math_mul33x33(crystallite_subF (1:3,1:3,c,i,e), & temp_33_3 = matmul(matmul(crystallite_subF (1:3,1:3,c,i,e), &
crystallite_invFp (1:3,1:3,c,i,e)), & crystallite_invFp (1:3,1:3,c,i,e)), &
math_inv33(crystallite_subFi0(1:3,1:3,c,i,e))) math_inv33(crystallite_subFi0(1:3,1:3,c,i,e)))
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt)
rhs_3333(p,o,1:3,1:3) = math_mul33x33(dSdFe(p,o,1:3,1:3),temp_33_1) rhs_3333(p,o,1:3,1:3) = matmul(dSdFe(p,o,1:3,1:3),temp_33_1)
temp_3333(1:3,1:3,p,o) = math_mul33x33(math_mul33x33(temp_33_2,dLpdS(1:3,1:3,p,o)), & temp_3333(1:3,1:3,p,o) = matmul(matmul(temp_33_2,dLpdS(1:3,1:3,p,o)), &
crystallite_invFi(1:3,1:3,c,i,e)) & crystallite_invFi(1:3,1:3,c,i,e)) &
+ math_mul33x33(temp_33_3,dLidS(1:3,1:3,p,o)) + matmul(temp_33_3,dLidS(1:3,1:3,p,o))
end forall end forall
lhs_3333 = crystallite_subdt(c,i,e)*math_mul3333xx3333(dSdFe,temp_3333) & lhs_3333 = crystallite_subdt(c,i,e)*math_mul3333xx3333(dSdFe,temp_3333) &
+ math_mul3333xx3333(dSdFi,dFidS) + math_mul3333xx3333(dSdFi,dFidS)
@ -809,20 +806,20 @@ subroutine crystallite_stressTangent()
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt)
dFpinvdF(1:3,1:3,p,o) & dFpinvdF(1:3,1:3,p,o) &
= -crystallite_subdt(c,i,e) & = -crystallite_subdt(c,i,e) &
* math_mul33x33(math_inv33(crystallite_subFp0(1:3,1:3,c,i,e)), & * matmul(math_inv33(crystallite_subFp0(1:3,1:3,c,i,e)), &
math_mul33x33(temp_3333(1:3,1:3,p,o),crystallite_invFi(1:3,1:3,c,i,e))) matmul(temp_3333(1:3,1:3,p,o),crystallite_invFi(1:3,1:3,c,i,e)))
end forall end forall
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! assemble dPdF ! assemble dPdF
temp_33_1 = math_mul33x33(crystallite_invFp(1:3,1:3,c,i,e), & temp_33_1 = matmul(crystallite_invFp(1:3,1:3,c,i,e), &
math_mul33x33(crystallite_S(1:3,1:3,c,i,e), & matmul(crystallite_S(1:3,1:3,c,i,e), &
transpose(crystallite_invFp(1:3,1:3,c,i,e)))) transpose(crystallite_invFp(1:3,1:3,c,i,e))))
temp_33_2 = math_mul33x33(crystallite_S(1:3,1:3,c,i,e), & temp_33_2 = matmul(crystallite_S(1:3,1:3,c,i,e), &
transpose(crystallite_invFp(1:3,1:3,c,i,e))) transpose(crystallite_invFp(1:3,1:3,c,i,e)))
temp_33_3 = math_mul33x33(crystallite_subF(1:3,1:3,c,i,e), & temp_33_3 = matmul(crystallite_subF(1:3,1:3,c,i,e), &
crystallite_invFp(1:3,1:3,c,i,e)) crystallite_invFp(1:3,1:3,c,i,e))
temp_33_4 = math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,c,i,e), & temp_33_4 = matmul(matmul(crystallite_subF(1:3,1:3,c,i,e), &
crystallite_invFp(1:3,1:3,c,i,e)), & crystallite_invFp(1:3,1:3,c,i,e)), &
crystallite_S(1:3,1:3,c,i,e)) crystallite_S(1:3,1:3,c,i,e))
@ -832,9 +829,9 @@ subroutine crystallite_stressTangent()
enddo enddo
forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt) forall(p=1_pInt:3_pInt, o=1_pInt:3_pInt)
crystallite_dPdF(1:3,1:3,p,o,c,i,e) = crystallite_dPdF(1:3,1:3,p,o,c,i,e) + & crystallite_dPdF(1:3,1:3,p,o,c,i,e) = crystallite_dPdF(1:3,1:3,p,o,c,i,e) + &
math_mul33x33(math_mul33x33(crystallite_subF(1:3,1:3,c,i,e),dFpinvdF(1:3,1:3,p,o)),temp_33_2) + & matmul(matmul(crystallite_subF(1:3,1:3,c,i,e),dFpinvdF(1:3,1:3,p,o)),temp_33_2) + &
math_mul33x33(math_mul33x33(temp_33_3,dSdF(1:3,1:3,p,o)),transpose(crystallite_invFp(1:3,1:3,c,i,e))) + & matmul(matmul(temp_33_3,dSdF(1:3,1:3,p,o)),transpose(crystallite_invFp(1:3,1:3,c,i,e))) + &
math_mul33x33(temp_33_4,transpose(dFpinvdF(1:3,1:3,p,o))) matmul(temp_33_4,transpose(dFpinvdF(1:3,1:3,p,o)))
end forall end forall
enddo; enddo enddo; enddo
@ -895,7 +892,6 @@ end subroutine crystallite_orientations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function crystallite_push33ToRef(ipc,ip,el, tensor33) function crystallite_push33ToRef(ipc,ip,el, tensor33)
use math, only: & use math, only: &
math_mul33x33, &
math_inv33, & math_inv33, &
math_EulerToR math_EulerToR
use material, only: & use material, only: &
@ -910,9 +906,9 @@ function crystallite_push33ToRef(ipc,ip,el, tensor33)
ip, & ! integration point index ip, & ! integration point index
ipc ! grain index ipc ! grain index
T = math_mul33x33(math_EulerToR(material_EulerAngles(1:3,ipc,ip,el)), & T = matmul(math_EulerToR(material_EulerAngles(1:3,ipc,ip,el)), &
transpose(math_inv33(crystallite_subF(1:3,1:3,ipc,ip,el)))) transpose(math_inv33(crystallite_subF(1:3,1:3,ipc,ip,el))))
crystallite_push33ToRef = math_mul33x33(transpose(T),math_mul33x33(tensor33,T)) crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T))
end function crystallite_push33ToRef end function crystallite_push33ToRef
@ -924,7 +920,6 @@ function crystallite_postResults(ipc, ip, el)
use math, only: & use math, only: &
math_qToEuler, & math_qToEuler, &
math_qToEulerAxisAngle, & math_qToEulerAxisAngle, &
math_mul33x33, &
math_det33, & math_det33, &
math_I3, & math_I3, &
inDeg inDeg
@ -1093,11 +1088,7 @@ logical function integrateStress(&
use constitutive, only: constitutive_LpAndItsTangents, & use constitutive, only: constitutive_LpAndItsTangents, &
constitutive_LiAndItsTangents, & constitutive_LiAndItsTangents, &
constitutive_SandItsTangents constitutive_SandItsTangents
use math, only: math_mul33x33, & use math, only: math_mul33xx33, &
#ifdef __PGI
norm2, &
#endif
math_mul33xx33, &
math_mul3333xx3333, & math_mul3333xx3333, &
math_inv33, & math_inv33, &
math_det33, & math_det33, &
@ -1203,7 +1194,7 @@ logical function integrateStress(&
#endif #endif
return return
endif failedInversionFp endif failedInversionFp
A = math_mul33x33(Fg_new,invFp_current) ! intermediate tensor needed later to calculate dFe_dLp A = matmul(Fg_new,invFp_current) ! intermediate tensor needed later to calculate dFe_dLp
invFi_current = math_inv33(crystallite_subFi0(1:3,1:3,ipc,ip,el)) invFi_current = math_inv33(crystallite_subFi0(1:3,1:3,ipc,ip,el))
failedInversionFi: if (all(dEq0(invFi_current))) then failedInversionFi: if (all(dEq0(invFi_current))) then
@ -1235,7 +1226,7 @@ logical function integrateStress(&
return return
endif LiLoopLimit endif LiLoopLimit
invFi_new = math_mul33x33(invFi_current,math_I3 - dt*Liguess) invFi_new = matmul(invFi_current,math_I3 - dt*Liguess)
Fi_new = math_inv33(invFi_new) Fi_new = math_inv33(invFi_new)
detInvFi = math_det33(invFi_new) detInvFi = math_det33(invFi_new)
@ -1260,7 +1251,7 @@ logical function integrateStress(&
!* calculate (elastic) 2nd Piola--Kirchhoff stress tensor and its tangent from constitutive law !* calculate (elastic) 2nd Piola--Kirchhoff stress tensor and its tangent from constitutive law
B = math_I3 - dt*Lpguess B = math_I3 - dt*Lpguess
Fe = math_mul33x33(math_mul33x33(A,B), invFi_new) Fe = matmul(matmul(A,B), invFi_new)
call constitutive_SandItsTangents(S, dS_dFe, dS_dFi, & call constitutive_SandItsTangents(S, dS_dFe, dS_dFi, &
Fe, Fi_new, ipc, ip, el) ! call constitutive law to calculate 2nd Piola-Kirchhoff stress and its derivative in unloaded configuration Fe, Fi_new, ipc, ip, el) ! call constitutive law to calculate 2nd Piola-Kirchhoff stress and its derivative in unloaded configuration
@ -1406,13 +1397,13 @@ logical function integrateStress(&
!* calculate Jacobian for correction term !* calculate Jacobian for correction term
if (mod(jacoCounterLi, iJacoLpresiduum) == 0_pInt) then if (mod(jacoCounterLi, iJacoLpresiduum) == 0_pInt) then
temp_33 = math_mul33x33(math_mul33x33(A,B),invFi_current) temp_33 = matmul(matmul(A,B),invFi_current)
forall(o=1_pInt:3_pInt,p=1_pInt:3_pInt) forall(o=1_pInt:3_pInt,p=1_pInt:3_pInt)
dFe_dLi(1:3,o,1:3,p) = -dt*math_I3(o,p)*temp_33 ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j) dFe_dLi(1:3,o,1:3,p) = -dt*math_I3(o,p)*temp_33 ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j)
dFi_dLi(1:3,o,1:3,p) = -dt*math_I3(o,p)*invFi_current dFi_dLi(1:3,o,1:3,p) = -dt*math_I3(o,p)*invFi_current
end forall end forall
forall(o=1_pInt:3_pInt,p=1_pInt:3_pInt) & forall(o=1_pInt:3_pInt,p=1_pInt:3_pInt) &
dFi_dLi(1:3,1:3,o,p) = math_mul33x33(math_mul33x33(Fi_new,dFi_dLi(1:3,1:3,o,p)),Fi_new) dFi_dLi(1:3,1:3,o,p) = matmul(matmul(Fi_new,dFi_dLi(1:3,1:3,o,p)),Fi_new)
dRLi_dLi = math_identity2nd(9_pInt) & dRLi_dLi = math_identity2nd(9_pInt) &
- math_3333to99(math_mul3333xx3333(dLi_dS, math_mul3333xx3333(dS_dFe, dFe_dLi) + & - math_3333to99(math_mul3333xx3333(dLi_dS, math_mul3333xx3333(dS_dFe, dFe_dLi) + &
@ -1449,7 +1440,7 @@ logical function integrateStress(&
enddo LiLoop enddo LiLoop
!* calculate new plastic and elastic deformation gradient !* calculate new plastic and elastic deformation gradient
invFp_new = math_mul33x33(invFp_current,B) invFp_new = matmul(invFp_current,B)
invFp_new = invFp_new / math_det33(invFp_new)**(1.0_pReal/3.0_pReal) ! regularize invFp_new = invFp_new / math_det33(invFp_new)**(1.0_pReal/3.0_pReal) ! regularize
Fp_new = math_inv33(invFp_new) Fp_new = math_inv33(invFp_new)
failedInversionInvFp: if (all(dEq0(Fp_new))) then failedInversionInvFp: if (all(dEq0(Fp_new))) then
@ -1465,13 +1456,13 @@ logical function integrateStress(&
#endif #endif
return return
endif failedInversionInvFp endif failedInversionInvFp
Fe_new = math_mul33x33(math_mul33x33(Fg_new,invFp_new),invFi_new) Fe_new = matmul(matmul(Fg_new,invFp_new),invFi_new)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! stress integration was successful ! stress integration was successful
integrateStress = .true. integrateStress = .true.
crystallite_P (1:3,1:3,ipc,ip,el) = math_mul33x33(math_mul33x33(Fg_new,invFp_new), & crystallite_P (1:3,1:3,ipc,ip,el) = matmul(matmul(Fg_new,invFp_new), &
math_mul33x33(S,transpose(invFp_new))) matmul(S,transpose(invFp_new)))
crystallite_S (1:3,1:3,ipc,ip,el) = S crystallite_S (1:3,1:3,ipc,ip,el) = S
crystallite_Lp (1:3,1:3,ipc,ip,el) = Lpguess crystallite_Lp (1:3,1:3,ipc,ip,el) = Lpguess
crystallite_Li (1:3,1:3,ipc,ip,el) = Liguess crystallite_Li (1:3,1:3,ipc,ip,el) = Liguess
@ -1489,9 +1480,9 @@ logical function integrateStress(&
write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST integrateStress >> P / MPa', & write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST integrateStress >> P / MPa', &
transpose(crystallite_P(1:3,1:3,ipc,ip,el))*1.0e-6_pReal transpose(crystallite_P(1:3,1:3,ipc,ip,el))*1.0e-6_pReal
write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST integrateStress >> Cauchy / MPa', & write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST integrateStress >> Cauchy / MPa', &
math_mul33x33(crystallite_P(1:3,1:3,ipc,ip,el), transpose(Fg_new)) * 1.0e-6_pReal / math_det33(Fg_new) matmul(crystallite_P(1:3,1:3,ipc,ip,el), transpose(Fg_new)) * 1.0e-6_pReal / math_det33(Fg_new)
write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST integrateStress >> Fe Lp Fe^-1', & write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST integrateStress >> Fe Lp Fe^-1', &
transpose(math_mul33x33(Fe_new, math_mul33x33(crystallite_Lp(1:3,1:3,ipc,ip,el), math_inv33(Fe_new)))) transpose(matmul(Fe_new, matmul(crystallite_Lp(1:3,1:3,ipc,ip,el), math_inv33(Fe_new))))
write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST integrateStress >> Fp',transpose(crystallite_Fp(1:3,1:3,ipc,ip,el)) write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST integrateStress >> Fp',transpose(crystallite_Fp(1:3,1:3,ipc,ip,el))
write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST integrateStress >> Fi',transpose(crystallite_Fi(1:3,1:3,ipc,ip,el)) write(6,'(a,/,3(12x,3(f12.7,1x)/))') '<< CRYST integrateStress >> Fi',transpose(crystallite_Fi(1:3,1:3,ipc,ip,el))
endif endif

4
src/grid_damage_spectral.f90
View File

@ -286,8 +286,6 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
use mesh, only: & use mesh, only: &
grid, & grid, &
grid3 grid3
use math, only: &
math_mul33x3
use spectral_utilities, only: & use spectral_utilities, only: &
scalarField_real, & scalarField_real, &
vectorField_real, & vectorField_real, &
@ -328,7 +326,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
cell = 0 cell = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1 cell = cell + 1
vectorField_real(1:3,i,j,k) = math_mul33x3(damage_nonlocal_getDiffusion33(1,cell) - D_ref, & vectorField_real(1:3,i,j,k) = matmul(damage_nonlocal_getDiffusion33(1,cell) - D_ref, &
vectorField_real(1:3,i,j,k)) vectorField_real(1:3,i,j,k))
enddo; enddo; enddo enddo; enddo; enddo
call utilities_FFTvectorForward call utilities_FFTvectorForward

3
src/grid_mech_FEM.f90
View File

@ -316,7 +316,6 @@ end function grid_mech_FEM_solution
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mech_FEM_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC) subroutine grid_mech_FEM_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC)
use math, only: & use math, only: &
math_mul33x33 ,&
math_rotate_backward33 math_rotate_backward33
use numerics, only: & use numerics, only: &
worldrank worldrank
@ -402,7 +401,7 @@ subroutine grid_mech_FEM_forward(guess,timeinc,timeinc_old,loadCaseTime,deformat
! calculate rate for aim ! calculate rate for aim
if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F
F_aimDot = & F_aimDot = &
F_aimDot + deformation_BC%maskFloat * math_mul33x33(deformation_BC%values, F_aim_lastInc) F_aimDot + deformation_BC%maskFloat * matmul(deformation_BC%values, F_aim_lastInc)
elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed
F_aimDot = & F_aimDot = &
F_aimDot + deformation_BC%maskFloat * deformation_BC%values F_aimDot + deformation_BC%maskFloat * deformation_BC%values

3
src/grid_mech_spectral_basic.f90
View File

@ -285,7 +285,6 @@ end function grid_mech_spectral_basic_solution
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_basic_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC) subroutine grid_mech_spectral_basic_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC)
use math, only: & use math, only: &
math_mul33x33 ,&
math_rotate_backward33 math_rotate_backward33
use numerics, only: & use numerics, only: &
worldrank worldrank
@ -370,7 +369,7 @@ subroutine grid_mech_spectral_basic_forward(guess,timeinc,timeinc_old,loadCaseTi
! calculate rate for aim ! calculate rate for aim
if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F
F_aimDot = & F_aimDot = &
F_aimDot + deformation_BC%maskFloat * math_mul33x33(deformation_BC%values, F_aim_lastInc) F_aimDot + deformation_BC%maskFloat * matmul(deformation_BC%values, F_aim_lastInc)
elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed
F_aimDot = & F_aimDot = &
F_aimDot + deformation_BC%maskFloat * deformation_BC%values F_aimDot + deformation_BC%maskFloat * deformation_BC%values

14
src/grid_mech_spectral_polarisation.f90
View File

@ -308,7 +308,6 @@ end function grid_mech_spectral_polarisation_solution
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC) subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loadCaseTime,deformation_BC,stress_BC,rotation_BC)
use math, only: & use math, only: &
math_mul33x33, &
math_mul3333xx33, & math_mul3333xx33, &
math_rotate_backward33 math_rotate_backward33
use numerics, only: & use numerics, only: &
@ -402,7 +401,7 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa
! calculate rate for aim ! calculate rate for aim
if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F if (deformation_BC%myType=='l') then ! calculate F_aimDot from given L and current F
F_aimDot = & F_aimDot = &
F_aimDot + deformation_BC%maskFloat * math_mul33x33(deformation_BC%values, F_aim_lastInc) F_aimDot + deformation_BC%maskFloat * matmul(deformation_BC%values, F_aim_lastInc)
elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed elseif(deformation_BC%myType=='fdot') then ! F_aimDot is prescribed
F_aimDot = & F_aimDot = &
F_aimDot + deformation_BC%maskFloat * deformation_BC%values F_aimDot + deformation_BC%maskFloat * deformation_BC%values
@ -435,9 +434,9 @@ subroutine grid_mech_spectral_polarisation_forward(guess,timeinc,timeinc_old,loa
else else
do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1)
F_lambda33 = reshape(F_tau(1:9,i,j,k)-F(1:9,i,j,k),[3,3]) F_lambda33 = reshape(F_tau(1:9,i,j,k)-F(1:9,i,j,k),[3,3])
F_lambda33 = math_mul3333xx33(S_scale,math_mul33x33(F_lambda33, & F_lambda33 = math_mul3333xx33(S_scale,matmul(F_lambda33, &
math_mul3333xx33(C_scale,& math_mul3333xx33(C_scale,&
math_mul33x33(transpose(F_lambda33),& matmul(transpose(F_lambda33),&
F_lambda33)-math_I3))*0.5_pReal)& F_lambda33)-math_I3))*0.5_pReal)&
+ math_I3 + math_I3
F_tau(1:9,i,j,k) = reshape(F_lambda33,[9])+F(1:9,i,j,k) F_tau(1:9,i,j,k) = reshape(F_lambda33,[9])+F(1:9,i,j,k)
@ -528,8 +527,7 @@ subroutine formResidual(in, FandF_tau, &
math_rotate_forward33, & math_rotate_forward33, &
math_rotate_backward33, & math_rotate_backward33, &
math_mul3333xx33, & math_mul3333xx33, &
math_invSym3333, & math_invSym3333
math_mul33x33
use debug, only: & use debug, only: &
debug_level, & debug_level, &
debug_spectral, & debug_spectral, &
@ -605,7 +603,7 @@ subroutine formResidual(in, FandF_tau, &
do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1)
tensorField_real(1:3,1:3,i,j,k) = & tensorField_real(1:3,1:3,i,j,k) = &
polarBeta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -& polarBeta*math_mul3333xx33(C_scale,F(1:3,1:3,i,j,k) - math_I3) -&
polarAlpha*math_mul33x33(F(1:3,1:3,i,j,k), & polarAlpha*matmul(F(1:3,1:3,i,j,k), &
math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3)) math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))
enddo; enddo; enddo enddo; enddo; enddo
@ -644,7 +642,7 @@ subroutine formResidual(in, FandF_tau, &
e = e + 1 e = e + 1
residual_F(1:3,1:3,i,j,k) = & residual_F(1:3,1:3,i,j,k) = &
math_mul3333xx33(math_invSym3333(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,e) + C_scale), & math_mul3333xx33(math_invSym3333(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,e) + C_scale), &
residual_F(1:3,1:3,i,j,k) - math_mul33x33(F(1:3,1:3,i,j,k), & residual_F(1:3,1:3,i,j,k) - matmul(F(1:3,1:3,i,j,k), &
math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))) & math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))) &
+ residual_F_tau(1:3,1:3,i,j,k) + residual_F_tau(1:3,1:3,i,j,k)
enddo; enddo; enddo enddo; enddo; enddo

4
src/grid_thermal_spectral.f90
View File

@ -295,8 +295,6 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
use mesh, only: & use mesh, only: &
grid, & grid, &
grid3 grid3
use math, only: &
math_mul33x3
use spectral_utilities, only: & use spectral_utilities, only: &
scalarField_real, & scalarField_real, &
vectorField_real, & vectorField_real, &
@ -338,7 +336,7 @@ subroutine formResidual(in,x_scal,f_scal,dummy,ierr)
cell = 0 cell = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1 cell = cell + 1
vectorField_real(1:3,i,j,k) = math_mul33x3(thermal_conduction_getConductivity33(1,cell) - D_ref, & vectorField_real(1:3,i,j,k) = matmul(thermal_conduction_getConductivity33(1,cell) - D_ref, &
vectorField_real(1:3,i,j,k)) vectorField_real(1:3,i,j,k))
enddo; enddo; enddo enddo; enddo; enddo
call utilities_FFTvectorForward call utilities_FFTvectorForward

9
src/homogenization_RGC.f90
View File

@ -932,8 +932,7 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function surfaceCorrection(avgF,instance,of) function surfaceCorrection(avgF,instance,of)
use math, only: & use math, only: &
math_invert33, & math_invert33
math_mul33x33
implicit none implicit none
real(pReal), dimension(3) :: surfaceCorrection real(pReal), dimension(3) :: surfaceCorrection
@ -947,7 +946,7 @@ function homogenization_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
integer(pInt) :: i,j,iBase integer(pInt) :: i,j,iBase
logical :: error logical :: error
call math_invert33(math_mul33x33(transpose(avgF),avgF),invC,detF,error) call math_invert33(matmul(transpose(avgF),avgF),invC,detF,error)
surfaceCorrection = 0.0_pReal surfaceCorrection = 0.0_pReal
do iBase = 1_pInt,3_pInt do iBase = 1_pInt,3_pInt
@ -1139,8 +1138,6 @@ end function relaxationVector
!> @brief identify the normal of an interface !> @brief identify the normal of an interface
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function interfaceNormal(intFace,instance,of) pure function interfaceNormal(intFace,instance,of)
use math, only: &
math_mul33x3
implicit none implicit none
real(pReal), dimension (3) :: interfaceNormal real(pReal), dimension (3) :: interfaceNormal
@ -1156,7 +1153,7 @@ pure function interfaceNormal(intFace,instance,of)
nPos = abs(intFace(1)) ! identify the position of the interface in global state array nPos = abs(intFace(1)) ! identify the position of the interface in global state array
interfaceNormal(nPos) = real(intFace(1)/abs(intFace(1)),pReal) ! get the normal vector w.r.t. cluster axis interfaceNormal(nPos) = real(intFace(1)/abs(intFace(1)),pReal) ! get the normal vector w.r.t. cluster axis
interfaceNormal = math_mul33x3(dependentState(instance)%orientation(1:3,1:3,of),interfaceNormal) ! map the normal vector into sample coordinate system (basis) interfaceNormal = matmul(dependentState(instance)%orientation(1:3,1:3,of),interfaceNormal) ! map the normal vector into sample coordinate system (basis)
end function interfaceNormal end function interfaceNormal

15
src/lattice.f90
View File

@ -655,7 +655,6 @@ subroutine lattice_initializeStructure(myPhase,CoverA)
use prec, only: & use prec, only: &
tol_math_check tol_math_check
use math, only: & use math, only: &
math_mul33x33, &
math_sym3333to66, & math_sym3333to66, &
math_Voigt66to3333, & math_Voigt66to3333, &
math_cross math_cross
@ -1141,8 +1140,7 @@ function lattice_C66_trans(Ntrans,C_parent66,structure_target, &
math_axisAngleToR, & math_axisAngleToR, &
math_sym3333to66, & math_sym3333to66, &
math_66toSym3333, & math_66toSym3333, &
math_rotate_forward3333, & math_rotate_forward3333
math_mul33x33
implicit none implicit none
integer, dimension(:), intent(in) :: Ntrans !< number of active twin systems per family integer, dimension(:), intent(in) :: Ntrans !< number of active twin systems per family
@ -1210,7 +1208,6 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc
INRAD, & INRAD, &
math_outer, & math_outer, &
math_cross, & math_cross, &
math_mul33x3, &
math_axisAngleToR math_axisAngleToR
implicit none implicit none
integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family integer, dimension(:), intent(in) :: Nslip !< number of active slip systems per family
@ -1232,7 +1229,7 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc
do i = 1,sum(Nslip) do i = 1,sum(Nslip)
direction = coordinateSystem(1:3,1,i) direction = coordinateSystem(1:3,1,i)
normal = coordinateSystem(1:3,2,i) normal = coordinateSystem(1:3,2,i)
np = math_mul33x3(math_axisAngleToR(direction,60.0_pReal*INRAD), normal) np = matmul(math_axisAngleToR(direction,60.0_pReal*INRAD), normal)
if (size(nonSchmidCoefficients)>0) nonSchmidMatrix(1:3,1:3,i) = nonSchmidMatrix(1:3,1:3,i) & if (size(nonSchmidCoefficients)>0) nonSchmidMatrix(1:3,1:3,i) = nonSchmidMatrix(1:3,1:3,i) &
+ nonSchmidCoefficients(1) * math_outer(direction, np) + nonSchmidCoefficients(1) * math_outer(direction, np)
if (size(nonSchmidCoefficients)>1) nonSchmidMatrix(1:3,1:3,i) = nonSchmidMatrix(1:3,1:3,i) & if (size(nonSchmidCoefficients)>1) nonSchmidMatrix(1:3,1:3,i) = nonSchmidMatrix(1:3,1:3,i) &
@ -2401,8 +2398,6 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc)
use math, only: & use math, only: &
math_cross, & math_cross, &
math_outer, & math_outer, &
math_mul33x33, &
math_mul33x3, &
math_axisAngleToR, & math_axisAngleToR, &
INRAD, & INRAD, &
MATH_I3 MATH_I3
@ -2508,8 +2503,8 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc)
U = (a_bcc/a_fcc)*math_outer(x,x) & U = (a_bcc/a_fcc)*math_outer(x,x) &
+ (a_bcc/a_fcc)*math_outer(y,y) * sqrt(2.0_pReal) & + (a_bcc/a_fcc)*math_outer(y,y) * sqrt(2.0_pReal) &
+ (a_bcc/a_fcc)*math_outer(z,z) * sqrt(2.0_pReal) + (a_bcc/a_fcc)*math_outer(z,z) * sqrt(2.0_pReal)
Q(1:3,1:3,i) = math_mul33x33(R,B) Q(1:3,1:3,i) = matmul(R,B)
S(1:3,1:3,i) = math_mul33x33(R,U) - MATH_I3 S(1:3,1:3,i) = matmul(R,U) - MATH_I3
enddo enddo
elseif (cOverA > 0.0_pReal .and. dEq0(a_bcc)) then ! fcc -> hex transformation elseif (cOverA > 0.0_pReal .and. dEq0(a_bcc)) then ! fcc -> hex transformation
ss = MATH_I3 ss = MATH_I3
@ -2525,7 +2520,7 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc)
Q(1:3,1,i) = x Q(1:3,1,i) = x
Q(1:3,2,i) = y Q(1:3,2,i) = y
Q(1:3,3,i) = z Q(1:3,3,i) = z
S(1:3,1:3,i) = math_mul33x33(Q(1:3,1:3,i), math_mul33x33(math_mul33x33(sd,ss), transpose(Q(1:3,1:3,i)))) - MATH_I3 ! ToDo: This is of interest for the Schmid matrix only S(1:3,1:3,i) = matmul(Q(1:3,1:3,i), matmul(matmul(sd,ss), transpose(Q(1:3,1:3,i)))) - MATH_I3 ! ToDo: This is of interest for the Schmid matrix only
enddo enddo
else else
call IO_error(0) !ToDo: define error call IO_error(0) !ToDo: define error

44
src/math.f90
View File

@ -277,7 +277,7 @@ subroutine math_check
! +++ check rotation sense of q and R +++ ! +++ check rotation sense of q and R +++
v = halton([2,8,5]) ! random vector v = halton([2,8,5]) ! random vector
R = math_qToR(q) R = math_qToR(q)
if (any(abs(math_mul33x3(R,v) - math_qRot(q,v)) > tol_math_check)) then if (any(abs(matmul(R,v) - math_qRot(q,v)) > tol_math_check)) then
write (error_msg, '(a)' ) 'R(q)*v has different sense than q*v' write (error_msg, '(a)' ) 'R(q)*v has different sense than q*v'
call IO_error(401,ext_msg=error_msg) call IO_error(401,ext_msg=error_msg)
endif endif
@ -700,7 +700,7 @@ pure function math_exp33(A,n)
do i = 1, merge(n,5,present(n)) do i = 1, merge(n,5,present(n))
invFac = invFac/real(i,pReal) ! invfac = 1/i! invFac = invFac/real(i,pReal) ! invfac = 1/i!
B = math_mul33x33(B,A) B = matmul(B,A)
math_exp33 = math_exp33 + invFac*B ! exp = SUM (A^i)/i! math_exp33 = math_exp33 + invFac*B ! exp = SUM (A^i)/i!
enddo enddo
@ -1754,7 +1754,7 @@ real(pReal) pure function math_EulerMisorientation(EulerA,EulerB)
real(pReal), dimension(3), intent(in) :: EulerA,EulerB real(pReal), dimension(3), intent(in) :: EulerA,EulerB
real(pReal) :: cosTheta real(pReal) :: cosTheta
cosTheta = (math_trace33(math_mul33x33(math_EulerToR(EulerB), & cosTheta = (math_trace33(matmul(math_EulerToR(EulerB), &
transpose(math_EulerToR(EulerA)))) - 1.0_pReal) * 0.5_pReal transpose(math_EulerToR(EulerA)))) - 1.0_pReal) * 0.5_pReal
math_EulerMisorientation = acos(math_clip(cosTheta,-1.0_pReal,1.0_pReal)) math_EulerMisorientation = acos(math_clip(cosTheta,-1.0_pReal,1.0_pReal))
@ -1807,7 +1807,7 @@ function math_sampleGaussOri(center,FWHM)
angle = math_EulerMisorientation([0.0_pReal,0.0_pReal,0.0_pReal],math_RtoEuler(R)) angle = math_EulerMisorientation([0.0_pReal,0.0_pReal,0.0_pReal],math_RtoEuler(R))
if (rnd(4) <= exp(-4.0_pReal*log(2.0_pReal)*(angle/FWHM)**2_pReal)) exit ! rejection sampling (Gaussian) if (rnd(4) <= exp(-4.0_pReal*log(2.0_pReal)*(angle/FWHM)**2_pReal)) exit ! rejection sampling (Gaussian)
enddo GaussConvolution enddo GaussConvolution
math_sampleGaussOri = math_RtoEuler(math_mul33x33(R,math_EulerToR(center))) math_sampleGaussOri = math_RtoEuler(matmul(R,math_EulerToR(center)))
endif endif
end function math_sampleGaussOri end function math_sampleGaussOri
@ -1842,7 +1842,7 @@ function math_sampleFiberOri(alpha,beta,FWHM)
R = math_EulerAxisAngleToR(math_crossproduct(fInC,fInS),-acos(dot_product(fInC,fInS))) !< rotation to align fiber axis in crystal and sample system R = math_EulerAxisAngleToR(math_crossproduct(fInC,fInS),-acos(dot_product(fInC,fInS))) !< rotation to align fiber axis in crystal and sample system
rnd = halton([7,10,3]) rnd = halton([7,10,3])
R = math_mul33x33(R,math_EulerAxisAngleToR(fInS,rnd(1)*2.0_pReal*PI)) !< additional rotation (0..360deg) perpendicular to fiber axis R = matmul(R,math_EulerAxisAngleToR(fInS,rnd(1)*2.0_pReal*PI)) !< additional rotation (0..360deg) perpendicular to fiber axis
if (FWHM > 0.1_pReal*INRAD) then if (FWHM > 0.1_pReal*INRAD) then
reducedTo2D: do i=1,3 reducedTo2D: do i=1,3
@ -1863,7 +1863,7 @@ function math_sampleFiberOri(alpha,beta,FWHM)
u(j) = fInS(j) u(j) = fInS(j)
rejectionSampling: if (rnd(3) <= exp(-4.0_pReal*log(2.0_pReal)*(angle/FWHM)**2_pReal)) then rejectionSampling: if (rnd(3) <= exp(-4.0_pReal*log(2.0_pReal)*(angle/FWHM)**2_pReal)) then
R = math_mul33x33(R,math_EulerAxisAngleToR(math_crossproduct(u,fInS),angle)) ! tilt around direction of smallest component R = matmul(R,math_EulerAxisAngleToR(math_crossproduct(u,fInS),angle)) ! tilt around direction of smallest component
exit exit
endif rejectionSampling endif rejectionSampling
rnd = halton([7,10,3]) rnd = halton([7,10,3])
@ -2079,23 +2079,23 @@ pure function math_eigenvectorBasisSym33(m)
N(1:3,1:3,2) = m-values(2)*math_I3 N(1:3,1:3,2) = m-values(2)*math_I3
N(1:3,1:3,3) = m-values(3)*math_I3 N(1:3,1:3,3) = m-values(3)*math_I3
twoSimilarEigenvalues: if(abs(values(1)-values(2)) < TOL) then twoSimilarEigenvalues: if(abs(values(1)-values(2)) < TOL) then
EB(1:3,1:3,3)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,2))/ & EB(1:3,1:3,3)=matmul(N(1:3,1:3,1),N(1:3,1:3,2))/ &
((values(3)-values(1))*(values(3)-values(2))) ((values(3)-values(1))*(values(3)-values(2)))
EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,3) EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,3)
elseif(abs(values(2)-values(3)) < TOL) then twoSimilarEigenvalues elseif(abs(values(2)-values(3)) < TOL) then twoSimilarEigenvalues
EB(1:3,1:3,1)=math_mul33x33(N(1:3,1:3,2),N(1:3,1:3,3))/ & EB(1:3,1:3,1)=matmul(N(1:3,1:3,2),N(1:3,1:3,3))/ &
((values(1)-values(2))*(values(1)-values(3))) ((values(1)-values(2))*(values(1)-values(3)))
EB(1:3,1:3,2)=math_I3-EB(1:3,1:3,1) EB(1:3,1:3,2)=math_I3-EB(1:3,1:3,1)
elseif(abs(values(3)-values(1)) < TOL) then twoSimilarEigenvalues elseif(abs(values(3)-values(1)) < TOL) then twoSimilarEigenvalues
EB(1:3,1:3,2)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,3))/ & EB(1:3,1:3,2)=matmul(N(1:3,1:3,1),N(1:3,1:3,3))/ &
((values(2)-values(1))*(values(2)-values(3))) ((values(2)-values(1))*(values(2)-values(3)))
EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,2) EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,2)
else twoSimilarEigenvalues else twoSimilarEigenvalues
EB(1:3,1:3,1)=math_mul33x33(N(1:3,1:3,2),N(1:3,1:3,3))/ & EB(1:3,1:3,1)=matmul(N(1:3,1:3,2),N(1:3,1:3,3))/ &
((values(1)-values(2))*(values(1)-values(3))) ((values(1)-values(2))*(values(1)-values(3)))
EB(1:3,1:3,2)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,3))/ & EB(1:3,1:3,2)=matmul(N(1:3,1:3,1),N(1:3,1:3,3))/ &
((values(2)-values(1))*(values(2)-values(3))) ((values(2)-values(1))*(values(2)-values(3)))
EB(1:3,1:3,3)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,2))/ & EB(1:3,1:3,3)=matmul(N(1:3,1:3,1),N(1:3,1:3,2))/ &
((values(3)-values(1))*(values(3)-values(2))) ((values(3)-values(1))*(values(3)-values(2)))
endif twoSimilarEigenvalues endif twoSimilarEigenvalues
endif threeSimilarEigenvalues endif threeSimilarEigenvalues
@ -2144,23 +2144,23 @@ pure function math_eigenvectorBasisSym33_log(m)
N(1:3,1:3,2) = m-values(2)*math_I3 N(1:3,1:3,2) = m-values(2)*math_I3
N(1:3,1:3,3) = m-values(3)*math_I3 N(1:3,1:3,3) = m-values(3)*math_I3
twoSimilarEigenvalues: if(abs(values(1)-values(2)) < TOL) then twoSimilarEigenvalues: if(abs(values(1)-values(2)) < TOL) then
EB(1:3,1:3,3)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,2))/ & EB(1:3,1:3,3)=matmul(N(1:3,1:3,1),N(1:3,1:3,2))/ &
((values(3)-values(1))*(values(3)-values(2))) ((values(3)-values(1))*(values(3)-values(2)))
EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,3) EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,3)
elseif(abs(values(2)-values(3)) < TOL) then twoSimilarEigenvalues elseif(abs(values(2)-values(3)) < TOL) then twoSimilarEigenvalues
EB(1:3,1:3,1)=math_mul33x33(N(1:3,1:3,2),N(1:3,1:3,3))/ & EB(1:3,1:3,1)=matmul(N(1:3,1:3,2),N(1:3,1:3,3))/ &
((values(1)-values(2))*(values(1)-values(3))) ((values(1)-values(2))*(values(1)-values(3)))
EB(1:3,1:3,2)=math_I3-EB(1:3,1:3,1) EB(1:3,1:3,2)=math_I3-EB(1:3,1:3,1)
elseif(abs(values(3)-values(1)) < TOL) then twoSimilarEigenvalues elseif(abs(values(3)-values(1)) < TOL) then twoSimilarEigenvalues
EB(1:3,1:3,2)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,3))/ & EB(1:3,1:3,2)=matmul(N(1:3,1:3,1),N(1:3,1:3,3))/ &
((values(2)-values(1))*(values(2)-values(3))) ((values(2)-values(1))*(values(2)-values(3)))
EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,2) EB(1:3,1:3,1)=math_I3-EB(1:3,1:3,2)
else twoSimilarEigenvalues else twoSimilarEigenvalues
EB(1:3,1:3,1)=math_mul33x33(N(1:3,1:3,2),N(1:3,1:3,3))/ & EB(1:3,1:3,1)=matmul(N(1:3,1:3,2),N(1:3,1:3,3))/ &
((values(1)-values(2))*(values(1)-values(3))) ((values(1)-values(2))*(values(1)-values(3)))
EB(1:3,1:3,2)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,3))/ & EB(1:3,1:3,2)=matmul(N(1:3,1:3,1),N(1:3,1:3,3))/ &
((values(2)-values(1))*(values(2)-values(3))) ((values(2)-values(1))*(values(2)-values(3)))
EB(1:3,1:3,3)=math_mul33x33(N(1:3,1:3,1),N(1:3,1:3,2))/ & EB(1:3,1:3,3)=matmul(N(1:3,1:3,1),N(1:3,1:3,2))/ &
((values(3)-values(1))*(values(3)-values(2))) ((values(3)-values(1))*(values(3)-values(2)))
endif twoSimilarEigenvalues endif twoSimilarEigenvalues
endif threeSimilarEigenvalues endif threeSimilarEigenvalues
@ -2186,14 +2186,14 @@ function math_rotationalPart33(m)
real(pReal), dimension(3,3) :: math_rotationalPart33 real(pReal), dimension(3,3) :: math_rotationalPart33
real(pReal), dimension(3,3) :: U , Uinv real(pReal), dimension(3,3) :: U , Uinv
U = math_eigenvectorBasisSym33(math_mul33x33(transpose(m),m)) U = math_eigenvectorBasisSym33(matmul(transpose(m),m))
Uinv = math_inv33(U) Uinv = math_inv33(U)
inversionFailed: if (all(dEq0(Uinv))) then inversionFailed: if (all(dEq0(Uinv))) then
math_rotationalPart33 = math_I3 math_rotationalPart33 = math_I3
call IO_warning(650) call IO_warning(650)
else inversionFailed else inversionFailed
math_rotationalPart33 = math_mul33x33(m,Uinv) math_rotationalPart33 = matmul(m,Uinv)
endif inversionFailed endif inversionFailed
end function math_rotationalPart33 end function math_rotationalPart33
@ -2586,7 +2586,7 @@ pure function math_rotate_forward33(tensor,rot_tensor)
real(pReal), dimension(3,3) :: math_rotate_forward33 real(pReal), dimension(3,3) :: math_rotate_forward33
real(pReal), dimension(3,3), intent(in) :: tensor, rot_tensor real(pReal), dimension(3,3), intent(in) :: tensor, rot_tensor
math_rotate_forward33 = math_mul33x33(rot_tensor,math_mul33x33(tensor,transpose(rot_tensor))) math_rotate_forward33 = matmul(rot_tensor,matmul(tensor,transpose(rot_tensor)))
end function math_rotate_forward33 end function math_rotate_forward33
@ -2600,7 +2600,7 @@ pure function math_rotate_backward33(tensor,rot_tensor)
real(pReal), dimension(3,3) :: math_rotate_backward33 real(pReal), dimension(3,3) :: math_rotate_backward33
real(pReal), dimension(3,3), intent(in) :: tensor, rot_tensor real(pReal), dimension(3,3), intent(in) :: tensor, rot_tensor
math_rotate_backward33 = math_mul33x33(transpose(rot_tensor),math_mul33x33(tensor,rot_tensor)) math_rotate_backward33 = matmul(transpose(rot_tensor),matmul(tensor,rot_tensor))
end function math_rotate_backward33 end function math_rotate_backward33

4
src/mesh_grid.f90
View File

@ -561,8 +561,6 @@ function mesh_nodesAroundCentres(gDim,Favg,centres) result(nodes)
debug_mesh, & debug_mesh, &
debug_level, & debug_level, &
debug_levelBasic debug_levelBasic
use math, only: &
math_mul33x3
implicit none implicit none
real(pReal), intent(in), dimension(:,:,:,:) :: & real(pReal), intent(in), dimension(:,:,:,:) :: &
@ -624,7 +622,7 @@ function mesh_nodesAroundCentres(gDim,Favg,centres) result(nodes)
lookup = me-diag+shift*iRes lookup = me-diag+shift*iRes
wrappedCentres(1:3,i+1_pInt, j+1_pInt, k+1_pInt) = & wrappedCentres(1:3,i+1_pInt, j+1_pInt, k+1_pInt) = &
centres(1:3,lookup(1)+1_pInt,lookup(2)+1_pInt,lookup(3)+1_pInt) & centres(1:3,lookup(1)+1_pInt,lookup(2)+1_pInt,lookup(3)+1_pInt) &
- math_mul33x3(Favg, real(shift,pReal)*gDim) - matmul(Favg, real(shift,pReal)*gDim)
endif endif
enddo; enddo; enddo enddo; enddo; enddo

7
src/plastic_dislotwin.f90
View File

@ -651,8 +651,7 @@ subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of)
math_eigenValuesVectorsSym, & math_eigenValuesVectorsSym, &
math_outer, & math_outer, &
math_symmetric33, & math_symmetric33, &
math_mul33xx33, & math_mul33xx33
math_mul33x3
implicit none implicit none
real(pReal), dimension(3,3), intent(out) :: Lp real(pReal), dimension(3,3), intent(out) :: Lp
@ -723,8 +722,8 @@ subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,instance,of)
call math_eigenValuesVectorsSym(Mp,eigValues,eigVectors,error) call math_eigenValuesVectorsSym(Mp,eigValues,eigVectors,error)
do i = 1,6 do i = 1,6
P_sb = 0.5_pReal * math_outer(math_mul33x3(eigVectors,sb_sComposition(1:3,i)),& P_sb = 0.5_pReal * math_outer(matmul(eigVectors,sb_sComposition(1:3,i)),&
math_mul33x3(eigVectors,sb_mComposition(1:3,i))) matmul(eigVectors,sb_mComposition(1:3,i)))
tau = math_mul33xx33(Mp,P_sb) tau = math_mul33xx33(Mp,P_sb)
significantShearBandStress: if (abs(tau) > tol_math_check) then significantShearBandStress: if (abs(tau) > tol_math_check) then

48
src/plastic_nonlocal.f90
View File

@ -838,8 +838,7 @@ subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el)
IO_error IO_error
use math, only: & use math, only: &
PI, & PI, &
math_mul33x3, & math_inner, &
math_mul3x3, &
math_inv33 math_inv33
#ifdef DEBUG #ifdef DEBUG
use debug, only: & use debug, only: &
@ -1004,10 +1003,10 @@ subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el)
neighbor_rhoTotal(2,:,n) = sum(abs(rho_neighbor(:,scr)),2) neighbor_rhoTotal(2,:,n) = sum(abs(rho_neighbor(:,scr)),2)
connection_latticeConf(1:3,n) = & connection_latticeConf(1:3,n) = &
math_mul33x3(invFe, mesh_ipCoordinates(1:3,neighbor_ip,neighbor_el) & matmul(invFe, mesh_ipCoordinates(1:3,neighbor_ip,neighbor_el) &
- mesh_ipCoordinates(1:3,ip,el)) - mesh_ipCoordinates(1:3,ip,el))
normal_latticeConf = math_mul33x3(transpose(invFp), mesh_ipAreaNormal(1:3,n,ip,el)) normal_latticeConf = matmul(transpose(invFp), mesh_ipAreaNormal(1:3,n,ip,el))
if (math_mul3x3(normal_latticeConf,connection_latticeConf(1:3,n)) < 0.0_pReal) & ! neighboring connection points in opposite direction to face normal: must be periodic image if (math_inner(normal_latticeConf,connection_latticeConf(1:3,n)) < 0.0_pReal) & ! neighboring connection points in opposite direction to face normal: must be periodic image
connection_latticeConf(1:3,n) = normal_latticeConf * mesh_ipVolume(ip,el)/mesh_ipArea(n,ip,el) ! instead take the surface normal scaled with the diameter of the cell connection_latticeConf(1:3,n) = normal_latticeConf * mesh_ipVolume(ip,el)/mesh_ipArea(n,ip,el) ! instead take the surface normal scaled with the diameter of the cell
else else
! local neighbor or different lattice structure or different constitution instance -> use central values instead ! local neighbor or different lattice structure or different constitution instance -> use central values instead
@ -1047,7 +1046,7 @@ subroutine plastic_nonlocal_dependentState(Fe, Fp, ip, el)
invConnections = math_inv33(connections) invConnections = math_inv33(connections)
if (all(dEq0(invConnections))) call IO_error(-1,ext_msg='back stress calculation: inversion error') if (all(dEq0(invConnections))) call IO_error(-1,ext_msg='back stress calculation: inversion error')
rhoExcessGradient(c) = math_mul3x3(m(1:3,s,c), math_mul33x3(invConnections,rhoExcessDifferences)) rhoExcessGradient(c) = math_inner(m(1:3,s,c), matmul(invConnections,rhoExcessDifferences))
enddo enddo
! ... plus gradient from deads ... ! ... plus gradient from deads ...
@ -1528,13 +1527,8 @@ subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
debug_e debug_e
#endif #endif
use math, only: & use math, only: &
#ifdef __PGI math_inner, &
norm2, &
#endif
math_mul3x3, &
math_mul33x3, &
math_mul33xx33, & math_mul33xx33, &
math_mul33x33, &
math_inv33, & math_inv33, &
math_det33, & math_det33, &
PI PI
@ -1756,7 +1750,7 @@ subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
m(1:3,1:ns,4) = prm%slip_transverse m(1:3,1:ns,4) = prm%slip_transverse
my_Fe = Fe(1:3,1:3,1,ip,el) my_Fe = Fe(1:3,1:3,1,ip,el)
my_F = math_mul33x33(my_Fe, Fp(1:3,1:3,1,ip,el)) my_F = matmul(my_Fe, Fp(1:3,1:3,1,ip,el))
neighbors: do n = 1,theMesh%elem%nIPneighbors neighbors: do n = 1,theMesh%elem%nIPneighbors
@ -1774,7 +1768,7 @@ subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient
neighbor_instance = phase_plasticityInstance(material_phase(1,neighbor_ip,neighbor_el)) neighbor_instance = phase_plasticityInstance(material_phase(1,neighbor_ip,neighbor_el))
neighbor_Fe = Fe(1:3,1:3,1,neighbor_ip,neighbor_el) neighbor_Fe = Fe(1:3,1:3,1,neighbor_ip,neighbor_el)
neighbor_F = math_mul33x33(neighbor_Fe, Fp(1:3,1:3,1,neighbor_ip,neighbor_el)) neighbor_F = matmul(neighbor_Fe, Fp(1:3,1:3,1,neighbor_ip,neighbor_el))
Favg = 0.5_pReal * (my_F + neighbor_F) Favg = 0.5_pReal * (my_F + neighbor_F)
else ! if no neighbor, take my value as average else ! if no neighbor, take my value as average
Favg = my_F Favg = my_F
@ -1809,9 +1803,9 @@ subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
where (neighbor_rhoSgl * mesh_ipVolume(neighbor_ip,neighbor_el) ** 0.667_pReal < prm%significantN & where (neighbor_rhoSgl * mesh_ipVolume(neighbor_ip,neighbor_el) ** 0.667_pReal < prm%significantN &
.or. neighbor_rhoSgl < prm%significantRho) & .or. neighbor_rhoSgl < prm%significantRho) &
neighbor_rhoSgl = 0.0_pReal neighbor_rhoSgl = 0.0_pReal
normal_neighbor2me_defConf = math_det33(Favg) * math_mul33x3(math_inv33(transpose(Favg)), & normal_neighbor2me_defConf = math_det33(Favg) * matmul(math_inv33(transpose(Favg)), &
mesh_ipAreaNormal(1:3,neighbor_n,neighbor_ip,neighbor_el)) ! calculate the normal of the interface in (average) deformed configuration (now pointing from my neighbor to me!!!) mesh_ipAreaNormal(1:3,neighbor_n,neighbor_ip,neighbor_el)) ! calculate the normal of the interface in (average) deformed configuration (now pointing from my neighbor to me!!!)
normal_neighbor2me = math_mul33x3(transpose(neighbor_Fe), normal_neighbor2me_defConf) & normal_neighbor2me = matmul(transpose(neighbor_Fe), normal_neighbor2me_defConf) &
/ math_det33(neighbor_Fe) ! interface normal in the lattice configuration of my neighbor / math_det33(neighbor_Fe) ! interface normal in the lattice configuration of my neighbor
area = mesh_ipArea(neighbor_n,neighbor_ip,neighbor_el) * norm2(normal_neighbor2me) area = mesh_ipArea(neighbor_n,neighbor_ip,neighbor_el) * norm2(normal_neighbor2me)
normal_neighbor2me = normal_neighbor2me / norm2(normal_neighbor2me) ! normalize the surface normal to unit length normal_neighbor2me = normal_neighbor2me / norm2(normal_neighbor2me) ! normalize the surface normal to unit length
@ -1819,10 +1813,10 @@ subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
do t = 1,4 do t = 1,4
c = (t + 1) / 2 c = (t + 1) / 2
topp = t + mod(t,2) - mod(t+1,2) topp = t + mod(t,2) - mod(t+1,2)
if (neighbor_v(s,t) * math_mul3x3(m(1:3,s,t), normal_neighbor2me) > 0.0_pReal & ! flux from my neighbor to me == entering flux for me if (neighbor_v(s,t) * math_inner(m(1:3,s,t), normal_neighbor2me) > 0.0_pReal & ! flux from my neighbor to me == entering flux for me
.and. v(s,t) * neighbor_v(s,t) >= 0.0_pReal ) then ! ... only if no sign change in flux density .and. v(s,t) * neighbor_v(s,t) >= 0.0_pReal ) then ! ... only if no sign change in flux density
lineLength = neighbor_rhoSgl(s,t) * neighbor_v(s,t) & lineLength = neighbor_rhoSgl(s,t) * neighbor_v(s,t) &
* math_mul3x3(m(1:3,s,t), normal_neighbor2me) * area ! positive line length that wants to enter through this interface * math_inner(m(1:3,s,t), normal_neighbor2me) * area ! positive line length that wants to enter through this interface
where (compatibility(c,1:ns,s,n,ip,el) > 0.0_pReal) & ! positive compatibility... where (compatibility(c,1:ns,s,n,ip,el) > 0.0_pReal) & ! positive compatibility...
rhoDotFlux(1:ns,t) = rhoDotFlux(1:ns,t) & rhoDotFlux(1:ns,t) = rhoDotFlux(1:ns,t) &
+ lineLength / mesh_ipVolume(ip,el) & ! ... transferring to equally signed mobile dislocation type + lineLength / mesh_ipVolume(ip,el) & ! ... transferring to equally signed mobile dislocation type
@ -1856,23 +1850,23 @@ subroutine plastic_nonlocal_dotState(Mp, Fe, Fp, Temperature, &
my_v = v my_v = v
normal_me2neighbor_defConf = math_det33(Favg) & normal_me2neighbor_defConf = math_det33(Favg) &
* math_mul33x3(math_inv33(transpose(Favg)), & * matmul(math_inv33(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!!!) 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) & normal_me2neighbor = matmul(transpose(my_Fe), normal_me2neighbor_defConf) &
/ math_det33(my_Fe) ! interface normal in my lattice configuration / math_det33(my_Fe) ! interface normal in my lattice configuration
area = mesh_ipArea(n,ip,el) * norm2(normal_me2neighbor) area = mesh_ipArea(n,ip,el) * norm2(normal_me2neighbor)
normal_me2neighbor = normal_me2neighbor / norm2(normal_me2neighbor) ! normalize the surface normal to unit length normal_me2neighbor = normal_me2neighbor / norm2(normal_me2neighbor) ! normalize the surface normal to unit length
do s = 1,ns do s = 1,ns
do t = 1,4 do t = 1,4
c = (t + 1) / 2 c = (t + 1) / 2
if (my_v(s,t) * math_mul3x3(m(1:3,s,t), normal_me2neighbor) > 0.0_pReal ) then ! flux from me to my neighbor == leaving flux for me (might also be a pure flux from my mobile density to dead density if interface not at all transmissive) if (my_v(s,t) * math_inner(m(1:3,s,t), normal_me2neighbor) > 0.0_pReal ) then ! flux from me to my neighbor == leaving flux for me (might also be a pure flux from my mobile density to dead density if interface not at all transmissive)
if (my_v(s,t) * neighbor_v(s,t) >= 0.0_pReal) then ! no sign change in flux density if (my_v(s,t) * neighbor_v(s,t) >= 0.0_pReal) then ! no sign change in flux density
transmissivity = sum(compatibility(c,1:ns,s,n,ip,el)**2.0_pReal) ! overall transmissivity from this slip system to my neighbor transmissivity = sum(compatibility(c,1:ns,s,n,ip,el)**2.0_pReal) ! overall transmissivity from this slip system to my neighbor
else ! sign change in flux density means sign change in stress which does not allow for dislocations to arive at the neighbor else ! sign change in flux density means sign change in stress which does not allow for dislocations to arive at the neighbor
transmissivity = 0.0_pReal transmissivity = 0.0_pReal
endif endif
lineLength = my_rhoSgl(s,t) * my_v(s,t) & lineLength = my_rhoSgl(s,t) * my_v(s,t) &
* math_mul3x3(m(1:3,s,t), normal_me2neighbor) * area ! positive line length of mobiles that wants to leave through this interface * math_inner(m(1:3,s,t), normal_me2neighbor) * area ! positive line length of mobiles that wants to leave through this interface
rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / mesh_ipVolume(ip,el) ! subtract dislocation flux from current type rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / mesh_ipVolume(ip,el) ! subtract dislocation flux from current type
rhoDotFlux(s,t+4) = rhoDotFlux(s,t+4) & rhoDotFlux(s,t+4) = rhoDotFlux(s,t+4) &
+ lineLength / mesh_ipVolume(ip,el) * (1.0_pReal - transmissivity) & + lineLength / mesh_ipVolume(ip,el) * (1.0_pReal - transmissivity) &
@ -2017,7 +2011,7 @@ end subroutine plastic_nonlocal_dotState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine plastic_nonlocal_updateCompatibility(orientation,i,e) subroutine plastic_nonlocal_updateCompatibility(orientation,i,e)
use math, only: & use math, only: &
math_mul3x3, & math_inner, &
math_qRot math_qRot
use rotations, only: & use rotations, only: &
rotation rotation
@ -2134,13 +2128,13 @@ subroutine plastic_nonlocal_updateCompatibility(orientation,i,e)
absoluteMisorientation = rot%asQuaternion() absoluteMisorientation = rot%asQuaternion()
mySlipSystems: do s1 = 1,ns mySlipSystems: do s1 = 1,ns
neighborSlipSystems: do s2 = 1,ns neighborSlipSystems: do s2 = 1,ns
my_compatibility(1,s2,s1,n) = math_mul3x3(prm%slip_normal(1:3,s1), & my_compatibility(1,s2,s1,n) = math_inner(prm%slip_normal(1:3,s1), &
math_qRot(absoluteMisorientation, prm%slip_normal(1:3,s2))) & math_qRot(absoluteMisorientation, prm%slip_normal(1:3,s2))) &
* abs(math_mul3x3(prm%slip_direction(1:3,s1), & * abs(math_inner(prm%slip_direction(1:3,s1), &
math_qRot(absoluteMisorientation, prm%slip_direction(1:3,s2)))) math_qRot(absoluteMisorientation, prm%slip_direction(1:3,s2))))
my_compatibility(2,s2,s1,n) = abs(math_mul3x3(prm%slip_normal(1:3,s1), & my_compatibility(2,s2,s1,n) = abs(math_inner(prm%slip_normal(1:3,s1), &
math_qRot(absoluteMisorientation, prm%slip_normal(1:3,s2)))) & math_qRot(absoluteMisorientation, prm%slip_normal(1:3,s2)))) &
* abs(math_mul3x3(prm%slip_direction(1:3,s1), & * abs(math_inner(prm%slip_direction(1:3,s1), &
math_qRot(absoluteMisorientation, prm%slip_direction(1:3,s2)))) math_qRot(absoluteMisorientation, prm%slip_direction(1:3,s2))))
enddo neighborSlipSystems enddo neighborSlipSystems

7
src/source_damage_isoBrittle.f90
View File

@ -174,8 +174,6 @@ subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el)
sourceState sourceState
use math, only : & use math, only : &
math_sym33to6, & math_sym33to6, &
math_mul33x33, &
math_mul66x6, &
math_I3 math_I3
implicit none implicit none
@ -200,9 +198,10 @@ subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el)
sourceOffset = source_damage_isoBrittle_offset(phase) sourceOffset = source_damage_isoBrittle_offset(phase)
strain = 0.5_pReal*math_sym33to6(math_mul33x33(transpose(Fe),Fe)-math_I3) strain = 0.5_pReal*math_sym33to6(matmul(transpose(Fe),Fe)-math_I3)
strainenergy = 2.0_pReal*sum(strain*math_mul66x6(C,strain))/param(instance)%critStrainEnergy strainenergy = 2.0_pReal*sum(strain*matmul(C,strain))/param(instance)%critStrainEnergy
! ToDo: check strainenergy = 2.0_pReal*dot_product(strain,matmul(C,strain))/param(instance)%critStrainEnergy
if (strainenergy > sourceState(phase)%p(sourceOffset)%subState0(1,constituent)) then if (strainenergy > sourceState(phase)%p(sourceOffset)%subState0(1,constituent)) then
sourceState(phase)%p(sourceOffset)%deltaState(1,constituent) = & sourceState(phase)%p(sourceOffset)%deltaState(1,constituent) = &

7
src/spectral_utilities.f90
View File

@ -610,7 +610,6 @@ end subroutine utilities_fourierGammaConvolution
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine utilities_fourierGreenConvolution(D_ref, mobility_ref, deltaT) subroutine utilities_fourierGreenConvolution(D_ref, mobility_ref, deltaT)
use math, only: & use math, only: &
math_mul33x3, &
PI PI
use mesh, only: & use mesh, only: &
grid, & grid, &
@ -1158,8 +1157,6 @@ subroutine utilities_updateIPcoords(F)
cNeq cNeq
use IO, only: & use IO, only: &
IO_error IO_error
use math, only: &
math_mul33x3
use mesh, only: & use mesh, only: &
grid, & grid, &
grid3, & grid3, &
@ -1200,12 +1197,12 @@ subroutine utilities_updateIPcoords(F)
if (grid3Offset == 0) offset_coords = vectorField_real(1:3,1,1,1) if (grid3Offset == 0) offset_coords = vectorField_real(1:3,1,1,1)
call MPI_Bcast(offset_coords,3,MPI_DOUBLE,0,PETSC_COMM_WORLD,ierr) call MPI_Bcast(offset_coords,3,MPI_DOUBLE,0,PETSC_COMM_WORLD,ierr)
if(ierr /=0) call IO_error(894, ext_msg='update_IPcoords') if(ierr /=0) call IO_error(894, ext_msg='update_IPcoords')
offset_coords = math_mul33x3(Favg,step/2.0_pReal) - offset_coords offset_coords = matmul(Favg,step/2.0_pReal) - offset_coords
m = 1 m = 1
do k = 1,grid3; do j = 1,grid(2); do i = 1,grid(1) do k = 1,grid3; do j = 1,grid(2); do i = 1,grid(1)
mesh_ipCoordinates(1:3,1,m) = vectorField_real(1:3,i,j,k) & mesh_ipCoordinates(1:3,1,m) = vectorField_real(1:3,i,j,k) &
+ offset_coords & + offset_coords &
+ math_mul33x3(Favg,step*real([i,j,k+grid3Offset]-1,pReal)) + matmul(Favg,step*real([i,j,k+grid3Offset]-1,pReal))
m = m+1 m = m+1
enddo; enddo; enddo enddo; enddo; enddo