polishing

This commit is contained in:
Martin Diehl 2020-12-23 08:12:56 +01:00
parent 7ee52afda2
commit 2947e7c444
2 changed files with 85 additions and 81 deletions

View File

@ -1250,97 +1250,97 @@ function crystallite_stressTangent(co,ip,el) result(dPdF)
pp = material_phaseAt(co,el) pp = material_phaseAt(co,el)
m = material_phaseMemberAt(co,ip,el) m = material_phaseMemberAt(co,ip,el)
call constitutive_hooke_SandItsTangents(devNull,dSdFe,dSdFi, & call constitutive_hooke_SandItsTangents(devNull,dSdFe,dSdFi, &
crystallite_Fe(1:3,1:3,co,ip,el), & crystallite_Fe(1:3,1:3,co,ip,el), &
constitutive_mech_Fi(pp)%data(1:3,1:3,m),co,ip,el) constitutive_mech_Fi(pp)%data(1:3,1:3,m),co,ip,el)
call constitutive_LiAndItsTangents(devNull,dLidS,dLidFi, & call constitutive_LiAndItsTangents(devNull,dLidS,dLidFi, &
crystallite_S (1:3,1:3,co,ip,el), & crystallite_S (1:3,1:3,co,ip,el), &
constitutive_mech_Fi(pp)%data(1:3,1:3,m), & constitutive_mech_Fi(pp)%data(1:3,1:3,m), &
co,ip,el) co,ip,el)
invFp = math_inv33(constitutive_mech_Fp(pp)%data(1:3,1:3,m)) invFp = math_inv33(constitutive_mech_Fp(pp)%data(1:3,1:3,m))
invFi = math_inv33(constitutive_mech_Fi(pp)%data(1:3,1:3,m)) invFi = math_inv33(constitutive_mech_Fi(pp)%data(1:3,1:3,m))
invSubFp0 = math_inv33(crystallite_subFp0(1:3,1:3,co,ip,el)) invSubFp0 = math_inv33(crystallite_subFp0(1:3,1:3,co,ip,el))
invSubFi0 = math_inv33(crystallite_subFi0(1:3,1:3,co,ip,el)) invSubFi0 = math_inv33(crystallite_subFi0(1:3,1:3,co,ip,el))
if (sum(abs(dLidS)) < tol_math_check) then if (sum(abs(dLidS)) < tol_math_check) then
dFidS = 0.0_pReal dFidS = 0.0_pReal
else else
lhs_3333 = 0.0_pReal; rhs_3333 = 0.0_pReal lhs_3333 = 0.0_pReal; rhs_3333 = 0.0_pReal
do o=1,3; do p=1,3 do o=1,3; do p=1,3
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(co,ip,el)*matmul(invSubFi0,dLidFi(1:3,1:3,o,p)) + crystallite_subdt(co,ip,el)*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) &
+ invFi*invFi(p,o) + invFi*invFi(p,o)
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(co,ip,el)*matmul(invSubFi0,dLidS(1:3,1:3,o,p)) - crystallite_subdt(co,ip,el)*matmul(invSubFi0,dLidS(1:3,1:3,o,p))
enddo; enddo enddo; enddo
call math_invert(temp_99,error,math_3333to99(lhs_3333)) call math_invert(temp_99,error,math_3333to99(lhs_3333))
if (error) then if (error) then
call IO_warning(warning_ID=600,el=el,ip=ip,g=co, & call IO_warning(warning_ID=600,el=el,ip=ip,g=co, &
ext_msg='inversion error in analytic tangent calculation') ext_msg='inversion error in analytic tangent calculation')
dFidS = 0.0_pReal dFidS = 0.0_pReal
else else
dFidS = math_mul3333xx3333(math_99to3333(temp_99),rhs_3333) dFidS = math_mul3333xx3333(math_99to3333(temp_99),rhs_3333)
endif endif
dLidS = math_mul3333xx3333(dLidFi,dFidS) + dLidS dLidS = math_mul3333xx3333(dLidFi,dFidS) + dLidS
endif endif
call constitutive_plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, & call constitutive_plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, &
crystallite_S (1:3,1:3,co,ip,el), & crystallite_S (1:3,1:3,co,ip,el), &
constitutive_mech_Fi(pp)%data(1:3,1:3,m),co,ip,el) constitutive_mech_Fi(pp)%data(1:3,1:3,m),co,ip,el)
dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculate dSdF ! calculate dSdF
temp_33_1 = transpose(matmul(invFp,invFi)) temp_33_1 = transpose(matmul(invFp,invFi))
temp_33_2 = matmul(crystallite_subF(1:3,1:3,co,ip,el),invSubFp0) temp_33_2 = matmul(crystallite_subF(1:3,1:3,co,ip,el),invSubFp0)
temp_33_3 = matmul(matmul(crystallite_subF(1:3,1:3,co,ip,el),invFp), invSubFi0) temp_33_3 = matmul(matmul(crystallite_subF(1:3,1:3,co,ip,el),invFp), invSubFi0)
do o=1,3; do p=1,3 do o=1,3; do p=1,3
rhs_3333(p,o,1:3,1:3) = matmul(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) = matmul(matmul(temp_33_2,dLpdS(1:3,1:3,p,o)), invFi) & temp_3333(1:3,1:3,p,o) = matmul(matmul(temp_33_2,dLpdS(1:3,1:3,p,o)), invFi) &
+ matmul(temp_33_3,dLidS(1:3,1:3,p,o)) + matmul(temp_33_3,dLidS(1:3,1:3,p,o))
enddo; enddo enddo; enddo
lhs_3333 = crystallite_subdt(co,ip,el)*math_mul3333xx3333(dSdFe,temp_3333) & lhs_3333 = crystallite_subdt(co,ip,el)*math_mul3333xx3333(dSdFe,temp_3333) &
+ math_mul3333xx3333(dSdFi,dFidS) + math_mul3333xx3333(dSdFi,dFidS)
call math_invert(temp_99,error,math_eye(9)+math_3333to99(lhs_3333)) call math_invert(temp_99,error,math_eye(9)+math_3333to99(lhs_3333))
if (error) then if (error) then
call IO_warning(warning_ID=600,el=el,ip=ip,g=co, & call IO_warning(warning_ID=600,el=el,ip=ip,g=co, &
ext_msg='inversion error in analytic tangent calculation') ext_msg='inversion error in analytic tangent calculation')
dSdF = rhs_3333 dSdF = rhs_3333
else else
dSdF = math_mul3333xx3333(math_99to3333(temp_99),rhs_3333) dSdF = math_mul3333xx3333(math_99to3333(temp_99),rhs_3333)
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculate dFpinvdF ! calculate dFpinvdF
temp_3333 = math_mul3333xx3333(dLpdS,dSdF) temp_3333 = math_mul3333xx3333(dLpdS,dSdF)
do o=1,3; do p=1,3 do o=1,3; do p=1,3
dFpinvdF(1:3,1:3,p,o) = -crystallite_subdt(co,ip,el) & dFpinvdF(1:3,1:3,p,o) = -crystallite_subdt(co,ip,el) &
* matmul(invSubFp0, matmul(temp_3333(1:3,1:3,p,o),invFi)) * matmul(invSubFp0, matmul(temp_3333(1:3,1:3,p,o),invFi))
enddo; enddo enddo; enddo
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! assemble dPdF ! assemble dPdF
temp_33_1 = matmul(crystallite_S(1:3,1:3,co,ip,el),transpose(invFp)) temp_33_1 = matmul(crystallite_S(1:3,1:3,co,ip,el),transpose(invFp))
temp_33_2 = matmul(invFp,temp_33_1) temp_33_2 = matmul(invFp,temp_33_1)
temp_33_3 = matmul(crystallite_subF(1:3,1:3,co,ip,el),invFp) temp_33_3 = matmul(crystallite_subF(1:3,1:3,co,ip,el),invFp)
temp_33_4 = matmul(temp_33_3,crystallite_S(1:3,1:3,co,ip,el)) temp_33_4 = matmul(temp_33_3,crystallite_S(1:3,1:3,co,ip,el))
dPdF = 0.0_pReal dPdF = 0.0_pReal
do p=1,3 do p=1,3
dPdF(p,1:3,p,1:3) = transpose(temp_33_2) dPdF(p,1:3,p,1:3) = transpose(temp_33_2)
enddo enddo
do o=1,3; do p=1,3 do o=1,3; do p=1,3
dPdF(1:3,1:3,p,o) = dPdF(1:3,1:3,p,o) & dPdF(1:3,1:3,p,o) = dPdF(1:3,1:3,p,o) &
+ matmul(matmul(crystallite_subF(1:3,1:3,co,ip,el), & + matmul(matmul(crystallite_subF(1:3,1:3,co,ip,el), &
dFpinvdF(1:3,1:3,p,o)),temp_33_1) & dFpinvdF(1:3,1:3,p,o)),temp_33_1) &
+ matmul(matmul(temp_33_3,dSdF(1:3,1:3,p,o)), & + matmul(matmul(temp_33_3,dSdF(1:3,1:3,p,o)), &
transpose(invFp)) & transpose(invFp)) &
+ matmul(temp_33_4,transpose(dFpinvdF(1:3,1:3,p,o))) + matmul(temp_33_4,transpose(dFpinvdF(1:3,1:3,p,o)))
enddo; enddo enddo; enddo
end function crystallite_stressTangent end function crystallite_stressTangent
@ -1385,14 +1385,16 @@ end subroutine crystallite_orientations
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function crystallite_push33ToRef(co,ip,el, tensor33) function crystallite_push33ToRef(co,ip,el, tensor33)
real(pReal), dimension(3,3) :: crystallite_push33ToRef
real(pReal), dimension(3,3), intent(in) :: tensor33 real(pReal), dimension(3,3), intent(in) :: tensor33
real(pReal), dimension(3,3) :: T real(pReal), dimension(3,3) :: T
integer, intent(in):: & integer, intent(in):: &
el, & el, &
ip, & ip, &
co co
real(pReal), dimension(3,3) :: crystallite_push33ToRef
T = matmul(material_orientation0(co,ip,el)%asMatrix(), & ! ToDo: initial orientation correct? T = matmul(material_orientation0(co,ip,el)%asMatrix(), & ! ToDo: initial orientation correct?
transpose(math_inv33(crystallite_subF(1:3,1:3,co,ip,el)))) transpose(math_inv33(crystallite_subF(1:3,1:3,co,ip,el))))
crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T)) crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T))
@ -1410,6 +1412,7 @@ subroutine integrateSourceState(co,ip,el)
el, & !< element index in element loop el, & !< element index in element loop
ip, & !< integration point index in ip loop ip, & !< integration point index in ip loop
co !< grain index in grain loop co !< grain index in grain loop
integer :: & integer :: &
NiterationState, & !< number of iterations in state loop NiterationState, & !< number of iterations in state loop
ph, & ph, &
@ -1426,6 +1429,7 @@ subroutine integrateSourceState(co,ip,el)
logical :: & logical :: &
broken broken
ph = material_phaseAt(co,el) ph = material_phaseAt(co,el)
c = material_phaseMemberAt(co,ip,el) c = material_phaseMemberAt(co,ip,el)

View File

@ -1185,16 +1185,15 @@ end subroutine integrateStateRKCK45
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine integrateStateRK(co,ip,el,A,B,CC,DB) subroutine integrateStateRK(co,ip,el,A,B,CC,DB)
real(pReal), dimension(:,:), intent(in) :: A real(pReal), dimension(:,:), intent(in) :: A
real(pReal), dimension(:), intent(in) :: B, CC real(pReal), dimension(:), intent(in) :: B, CC
real(pReal), dimension(:), intent(in), optional :: DB real(pReal), dimension(:), intent(in), optional :: DB
integer, intent(in) :: & integer, intent(in) :: &
el, & !< element index in element loop el, & !< element index in element loop
ip, & !< integration point index in ip loop ip, & !< integration point index in ip loop
co !< grain index in grain loop co !< grain index in grain loop
integer :: &
integer :: &
stage, & ! stage index in integration stage loop stage, & ! stage index in integration stage loop
n, & n, &
ph, & ph, &
@ -1202,7 +1201,8 @@ subroutine integrateStateRK(co,ip,el,A,B,CC,DB)
sizeDotState sizeDotState
logical :: & logical :: &
broken broken
real(pReal), dimension(constitutive_plasticity_maxSizeDotState,size(B)) :: plastic_RKdotState real(pReal), dimension(constitutive_plasticity_maxSizeDotState,size(B)) :: plastic_RKdotState
ph = material_phaseAt(co,el) ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el) me = material_phaseMemberAt(co,ip,el)