Merge branch 'development' into 19-NewStylePhenopowerlaw

This commit is contained in:
Martin Diehl 2018-09-16 22:03:05 +02:00
commit 5f06a35900
2 changed files with 206 additions and 143 deletions

View File

@ -19,16 +19,16 @@ module constitutive
constitutive_init, &
constitutive_homogenizedC, &
constitutive_microstructure, &
constitutive_LpAndItsTangent, &
constitutive_LiAndItsTangent, &
constitutive_LpAndItsTangents, &
constitutive_LiAndItsTangents, &
constitutive_initialFi, &
constitutive_TandItsTangent, &
constitutive_SandItsTangents, &
constitutive_collectDotState, &
constitutive_collectDeltaState, &
constitutive_postResults
private :: &
constitutive_hooke_TandItsTangent
constitutive_hooke_SandItsTangents
contains
@ -346,6 +346,7 @@ end subroutine constitutive_init
!--------------------------------------------------------------------------------------------------
!> @brief returns the homogenize elasticity matrix
!> ToDo: homogenizedC66 would be more consistent
!--------------------------------------------------------------------------------------------------
function constitutive_homogenizedC(ipc,ip,el)
use prec, only: &
@ -430,7 +431,7 @@ end subroutine constitutive_microstructure
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient
!--------------------------------------------------------------------------------------------------
subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v, Fi, ipc, ip, el)
subroutine constitutive_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, S6, Fi, ipc, ip, el)
use prec, only: &
pReal
use math, only: &
@ -470,20 +471,21 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v
ip, & !< integration point
el !< element
real(pReal), intent(in), dimension(6) :: &
Tstar_v !< 2nd Piola-Kirchhoff stress
S6 !< 2nd Piola-Kirchhoff stress (vector notation)
real(pReal), intent(in), dimension(3,3) :: &
Fi !< intermediate deformation gradient
real(pReal), intent(out), dimension(3,3) :: &
Lp !< plastic velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
dLp_dTstar3333, & !< derivative of Lp with respect to Tstar (4th-order tensor)
dLp_dFi3333 !< derivative of Lp with respect to Fi (4th-order tensor)
real(pReal), dimension(6) :: &
Mstar_v !< Mandel stress work conjugate with Lp
dLp_dS, &
dLp_dFi !< derivative of Lp with respect to Fi
real(pReal), dimension(3,3,3,3) :: &
dLp_dMp !< derivative of Lp with respect to Mandel stress
real(pReal), dimension(9,9) :: &
dLp_dMstar !< derivative of Lp with respect to Mstar (4th-order tensor)
dLp_dMp99 !< derivative of Lp with respect to Mstar (matrix notation)
real(pReal), dimension(3,3) :: &
temp_33
Mp, & !< Mandel stress work conjugate with Lp
S !< 2nd Piola-Kirchhoff stress
integer(pInt) :: &
ho, & !< homogenization
tme !< thermal member position
@ -493,47 +495,57 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v
ho = material_homog(ip,el)
tme = thermalMapping(ho)%p(ip,el)
Mstar_v = math_Mandel33to6(math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(Tstar_v)))
S = math_Mandel6to33(S6)
Mp = math_mul33x33(math_mul33x33(transpose(Fi),Fi),S)
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_NONE_ID) plasticityType
Lp = 0.0_pReal
dLp_dMstar = 0.0_pReal
dLp_dMp = 0.0_pReal
case (PLASTICITY_ISOTROPIC_ID) plasticityType
call plastic_isotropic_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v,ipc,ip,el)
call plastic_isotropic_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp),ipc,ip,el)
dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
call plastic_phenopowerlaw_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v,ipc,ip,el)
call plastic_phenopowerlaw_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp),ipc,ip,el)
dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget
case (PLASTICITY_KINEHARDENING_ID) plasticityType
call plastic_kinehardening_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v,ipc,ip,el)
call plastic_kinehardening_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp),ipc,ip,el)
dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget
case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v, &
call plastic_nonlocal_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp), &
temperature(ho)%p(tme),ip,el)
dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget
case (PLASTICITY_DISLOTWIN_ID) plasticityType
call plastic_dislotwin_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v, &
call plastic_dislotwin_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp), &
temperature(ho)%p(tme),ipc,ip,el)
dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget
case (PLASTICITY_DISLOUCLA_ID) plasticityType
call plastic_disloucla_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v, &
call plastic_disloucla_LpAndItsTangent (Lp,dLp_dMp99, math_Mandel33to6(Mp), &
temperature(ho)%p(tme), ipc,ip,el)
dLp_dMp = math_Plain99to3333(dLp_dMp99) ! ToDo: We revert here the last statement in plastic_xx_LpAndItsTanget
end select plasticityType
dLp_dTstar3333 = math_Plain99to3333(dLp_dMstar)
temp_33 = math_mul33x33(Fi,math_Mandel6to33(Tstar_v))
forall(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt) &
dLp_dFi3333(i,j,1:3,1:3) = math_mul33x33(temp_33,transpose(dLp_dTstar3333(i,j,1:3,1:3))) + &
math_mul33x33(math_mul33x33(Fi,dLp_dTstar3333(i,j,1:3,1:3)),math_Mandel6to33(Tstar_v))
do concurrent(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt)
dLp_dFi(i,j,1:3,1:3) = math_mul33x33(math_mul33x33(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)
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)
enddo
temp_33 = math_mul33x33(transpose(Fi),Fi)
forall(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt) &
dLp_dTstar3333(i,j,1:3,1:3) = math_mul33x33(temp_33,dLp_dTstar3333(i,j,1:3,1:3))
end subroutine constitutive_LpAndItsTangent
end subroutine constitutive_LpAndItsTangents
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient
!--------------------------------------------------------------------------------------------------
subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar3333, dLi_dFi3333, Tstar_v, Fi, ipc, ip, el)
subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, S6, Fi, ipc, ip, el)
use prec, only: &
pReal
use math, only: &
@ -571,18 +583,18 @@ subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar3333, dLi_dFi3333, Tstar_v
ip, & !< integration point
el !< element
real(pReal), intent(in), dimension(6) :: &
Tstar_v !< 2nd Piola-Kirchhoff stress
S6 !< 2nd Piola-Kirchhoff stress (vector notation)
real(pReal), intent(in), dimension(3,3) :: &
Fi !< intermediate deformation gradient
real(pReal), intent(out), dimension(3,3) :: &
Li !< intermediate velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: &
dLi_dTstar3333, & !< derivative of Li with respect to Tstar (4th-order tensor)
dLi_dFi3333
dLi_dS, & !< derivative of Li with respect to S
dLi_dFi
real(pReal), dimension(3,3) :: &
my_Li !< intermediate velocity gradient
real(pReal), dimension(3,3,3,3) :: &
my_dLi_dTstar
my_dLi_dS
real(pReal), dimension(3,3) :: &
FiInv, &
temp_33
@ -594,51 +606,51 @@ subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar3333, dLi_dFi3333, Tstar_v
i, j
Li = 0.0_pReal
dLi_dTstar3333 = 0.0_pReal
dLi_dFi3333 = 0.0_pReal
dLi_dS = 0.0_pReal
dLi_dFi = 0.0_pReal
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_isotropic_ID) plasticityType
call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dTstar, Tstar_v, ipc, ip, el)
call plastic_isotropic_LiAndItsTangent(my_Li, my_dLi_dS, S6, ipc, ip, el)
case default plasticityType
my_Li = 0.0_pReal
my_dLi_dTstar = 0.0_pReal
my_dLi_dS = 0.0_pReal
end select plasticityType
Li = Li + my_Li
dLi_dTstar3333 = dLi_dTstar3333 + my_dLi_dTstar
dLi_dS = dLi_dS + my_dLi_dS
KinematicsLoop: do k = 1_pInt, phase_Nkinematics(material_phase(ipc,ip,el))
kinematicsType: select case (phase_kinematics(k,material_phase(ipc,ip,el)))
case (KINEMATICS_cleavage_opening_ID) kinematicsType
call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dTstar, Tstar_v, ipc, ip, el)
call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S6, ipc, ip, el)
case (KINEMATICS_slipplane_opening_ID) kinematicsType
call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dTstar, Tstar_v, ipc, ip, el)
call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S6, ipc, ip, el)
case (KINEMATICS_thermal_expansion_ID) kinematicsType
call kinematics_thermal_expansion_LiAndItsTangent(my_Li, my_dLi_dTstar, ipc, ip, el)
call kinematics_thermal_expansion_LiAndItsTangent(my_Li, my_dLi_dS, ipc, ip, el)
case (KINEMATICS_vacancy_strain_ID) kinematicsType
call kinematics_vacancy_strain_LiAndItsTangent(my_Li, my_dLi_dTstar, ipc, ip, el)
call kinematics_vacancy_strain_LiAndItsTangent(my_Li, my_dLi_dS, ipc, ip, el)
case (KINEMATICS_hydrogen_strain_ID) kinematicsType
call kinematics_hydrogen_strain_LiAndItsTangent(my_Li, my_dLi_dTstar, ipc, ip, el)
call kinematics_hydrogen_strain_LiAndItsTangent(my_Li, my_dLi_dS, ipc, ip, el)
case default kinematicsType
my_Li = 0.0_pReal
my_dLi_dTstar = 0.0_pReal
my_dLi_dS = 0.0_pReal
end select kinematicsType
Li = Li + my_Li
dLi_dTstar3333 = dLi_dTstar3333 + my_dLi_dTstar
dLi_dS = dLi_dS + my_dLi_dS
enddo KinematicsLoop
FiInv = math_inv33(Fi)
detFi = math_det33(Fi)
Li = math_mul33x33(math_mul33x33(Fi,Li),FiInv)*detFi !< push forward to intermediate configuration
temp_33 = math_mul33x33(FiInv,Li)
forall(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt)
dLi_dTstar3333(1:3,1:3,i,j) = math_mul33x33(math_mul33x33(Fi,dLi_dTstar3333(1:3,1:3,i,j)),FiInv)*detFi
dLi_dFi3333 (1:3,1:3,i,j) = dLi_dFi3333(1:3,1:3,i,j) + Li*FiInv(j,i)
dLi_dFi3333 (1:3,i,1:3,j) = dLi_dFi3333(1:3,i,1:3,j) + math_I3*temp_33(j,i) + Li*FiInv(j,i)
end forall
do concurrent(i = 1_pInt:3_pInt, 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_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)
end do
end subroutine constitutive_LiAndItsTangent
end subroutine constitutive_LiAndItsTangents
!--------------------------------------------------------------------------------------------------
@ -696,10 +708,10 @@ end function constitutive_initialFi
!--------------------------------------------------------------------------------------------------
!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to
!> the elastic deformation gradient depending on the selected elastic law (so far no case switch
!! because only Hooke is implemented
!> the elastic/intermediate deformation gradients depending on the selected elastic law
!! (so far no case switch because only Hooke is implemented)
!--------------------------------------------------------------------------------------------------
subroutine constitutive_TandItsTangent(T, dT_dFe, dT_dFi, Fe, Fi, ipc, ip, el)
subroutine constitutive_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip, el)
use prec, only: &
pReal
@ -712,30 +724,28 @@ subroutine constitutive_TandItsTangent(T, dT_dFe, dT_dFi, Fe, Fi, ipc, ip, el)
Fe, & !< elastic deformation gradient
Fi !< intermediate deformation gradient
real(pReal), intent(out), dimension(3,3) :: &
T !< 2nd Piola-Kirchhoff stress tensor
S !< 2nd Piola-Kirchhoff stress tensor
real(pReal), intent(out), dimension(3,3,3,3) :: &
dT_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient
dT_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient
dS_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient
dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient
call constitutive_hooke_TandItsTangent(T, dT_dFe, dT_dFi, Fe, Fi, ipc, ip, el)
call constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip, el)
end subroutine constitutive_TandItsTangent
end subroutine constitutive_SandItsTangents
!--------------------------------------------------------------------------------------------------
!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to
!> the elastic deformation gradient using hookes law
!> the elastic and intermeidate deformation gradients using Hookes law
!--------------------------------------------------------------------------------------------------
subroutine constitutive_hooke_TandItsTangent(T, dT_dFe, dT_dFi, Fe, Fi, ipc, ip, el)
subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip, el)
use prec, only: &
pReal
use math, only : &
math_mul3x3, &
math_mul33x33, &
math_mul3333xx33, &
math_Mandel66to3333, &
math_trace33, &
math_I3
use material, only: &
material_phase, &
@ -758,10 +768,10 @@ subroutine constitutive_hooke_TandItsTangent(T, dT_dFe, dT_dFi, Fe, Fi, ipc, ip,
Fe, & !< elastic deformation gradient
Fi !< intermediate deformation gradient
real(pReal), intent(out), dimension(3,3) :: &
T !< 2nd Piola-Kirchhoff stress tensor in lattice configuration
S !< 2nd Piola-Kirchhoff stress tensor in lattice configuration
real(pReal), intent(out), dimension(3,3,3,3) :: &
dT_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient
dT_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient
dS_dFe, & !< derivative of 2nd P-K stress with respect to elastic deformation gradient
dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient
real(pReal), dimension(3,3) :: E
real(pReal), dimension(3,3,3,3) :: C
integer(pInt) :: &
@ -784,22 +794,22 @@ subroutine constitutive_hooke_TandItsTangent(T, dT_dFe, dT_dFi, Fe, Fi, ipc, ip,
enddo DegradationLoop
E = 0.5_pReal*(math_mul33x33(transpose(Fe),Fe)-math_I3) !< Green-Lagrange strain in unloaded configuration
T = 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,math_mul33x33(math_mul33x33(transpose(Fi),E),Fi)) !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration
dT_dFe = 0.0_pReal
dS_dFe = 0.0_pReal
forall (i=1_pInt:3_pInt, j=1_pInt:3_pInt)
dT_dFe(i,j,1:3,1:3) = &
math_mul33x33(Fe,math_mul33x33(math_mul33x33(Fi,C(i,j,1:3,1:3)),transpose(Fi))) !< dT_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko
dT_dFi(i,j,1:3,1:3) = 2.0_pReal*math_mul33x33(math_mul33x33(E,Fi),C(i,j,1:3,1:3)) !< dT_ij/dFi_kl = C_ijln * E_km * Fe_mn
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
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
end forall
end subroutine constitutive_hooke_TandItsTangent
end subroutine constitutive_hooke_SandItsTangents
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfracArray,ipc, ip, el)
subroutine constitutive_collectDotState(S6, FeArray, Fi, FpArray, subdt, subfracArray,ipc, ip, el)
use prec, only: &
pReal, &
pLongInt
@ -807,6 +817,10 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra
debug_level, &
debug_constitutive, &
debug_levelBasic
use math, only: &
math_Mandel6to33, &
math_Mandel33to6, &
math_mul33x33
use mesh, only: &
mesh_NcpElems, &
mesh_maxNips
@ -863,8 +877,12 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra
real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
FeArray, & !< elastic deformation gradient
FpArray !< plastic deformation gradient
real(pReal), intent(in), dimension(3,3) :: &
Fi !< intermediate deformation gradient
real(pReal), intent(in), dimension(6) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel)
S6 !< 2nd Piola Kirchhoff stress (vector notation)
real(pReal), dimension(3,3) :: &
Mstar
integer(pInt) :: &
ho, & !< homogenization
tme, & !< thermal member position
@ -873,35 +891,50 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra
ho = material_homog( ip,el)
tme = thermalMapping(ho)%p(ip,el)
Mstar = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(S6))
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_ISOTROPIC_ID) plasticityType
call plastic_isotropic_dotState (Tstar_v,ipc,ip,el)
call plastic_isotropic_dotState (math_Mandel33to6(Mstar),ipc,ip,el)
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
call plastic_phenopowerlaw_dotState(Tstar_v,ipc,ip,el)
call plastic_phenopowerlaw_dotState(math_Mandel33to6(Mstar),ipc,ip,el)
case (PLASTICITY_KINEHARDENING_ID) plasticityType
call plastic_kinehardening_dotState(Tstar_v,ipc,ip,el)
call plastic_kinehardening_dotState(math_Mandel33to6(Mstar),ipc,ip,el)
case (PLASTICITY_DISLOTWIN_ID) plasticityType
call plastic_dislotwin_dotState (Tstar_v,temperature(ho)%p(tme), &
call plastic_dislotwin_dotState (math_Mandel33to6(Mstar),temperature(ho)%p(tme), &
ipc,ip,el)
case (PLASTICITY_DISLOUCLA_ID) plasticityType
call plastic_disloucla_dotState (Tstar_v,temperature(ho)%p(tme), &
call plastic_disloucla_dotState (math_Mandel33to6(Mstar),temperature(ho)%p(tme), &
ipc,ip,el)
case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_dotState (Tstar_v,FeArray,FpArray,temperature(ho)%p(tme), &
call plastic_nonlocal_dotState (math_Mandel33to6(Mstar),FeArray,FpArray,temperature(ho)%p(tme), &
subdt,subfracArray,ip,el)
end select plasticityType
SourceLoop: do s = 1_pInt, phase_Nsources(material_phase(ipc,ip,el))
sourceType: select case (phase_source(s,material_phase(ipc,ip,el)))
case (SOURCE_damage_anisoBrittle_ID) sourceType
call source_damage_anisoBrittle_dotState (Tstar_v, ipc, ip, el)
call source_damage_anisoBrittle_dotState (S6, ipc, ip, el) !< correct stress?
case (SOURCE_damage_isoDuctile_ID) sourceType
call source_damage_isoDuctile_dotState ( ipc, ip, el)
case (SOURCE_damage_anisoDuctile_ID) sourceType
call source_damage_anisoDuctile_dotState ( ipc, ip, el)
case (SOURCE_thermal_externalheat_ID) sourceType
call source_thermal_externalheat_dotState( ipc, ip, el)
end select sourceType
enddo SourceLoop
end subroutine constitutive_collectDotState
@ -910,7 +943,7 @@ end subroutine constitutive_collectDotState
!> @brief for constitutive models having an instantaneous change of state
!> will return false if delta state is not needed/supported by the constitutive model
!--------------------------------------------------------------------------------------------------
subroutine constitutive_collectDeltaState(Tstar_v, Fe, ipc, ip, el)
subroutine constitutive_collectDeltaState(S6, Fe, Fi, ipc, ip, el)
use prec, only: &
pReal, &
pLongInt
@ -918,6 +951,10 @@ subroutine constitutive_collectDeltaState(Tstar_v, Fe, ipc, ip, el)
debug_level, &
debug_constitutive, &
debug_levelBasic
use math, only: &
math_Mandel6to33, &
math_Mandel33to6, &
math_mul33x33
use material, only: &
phase_plasticity, &
phase_source, &
@ -945,29 +982,43 @@ subroutine constitutive_collectDeltaState(Tstar_v, Fe, ipc, ip, el)
ip, & !< integration point
el !< element
real(pReal), intent(in), dimension(6) :: &
Tstar_v !< 2nd Piola-Kirchhoff stress
S6 !< 2nd Piola Kirchhoff stress (vector notation)
real(pReal), intent(in), dimension(3,3) :: &
Fe !< elastic deformation gradient
Fe, & !< elastic deformation gradient
Fi !< intermediate deformation gradient
real(pReal), dimension(3,3) :: &
Mstar
integer(pInt) :: &
s !< counter in source loop
Mstar = math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(S6))
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_KINEHARDENING_ID) plasticityType
call plastic_kinehardening_deltaState(Tstar_v,ipc,ip,el)
call plastic_kinehardening_deltaState(math_Mandel33to6(Mstar),ipc,ip,el)
case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_deltaState(Tstar_v,ip,el)
call plastic_nonlocal_deltaState(math_Mandel33to6(Mstar),ip,el)
end select plasticityType
SourceLoop: do s = 1_pInt, phase_Nsources(material_phase(ipc,ip,el))
sourceLoop: do s = 1_pInt, phase_Nsources(material_phase(ipc,ip,el))
sourceType: select case (phase_source(s,material_phase(ipc,ip,el)))
case (SOURCE_damage_isoBrittle_ID) sourceType
call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(ipc,ip,el), Fe, &
ipc, ip, el)
case (SOURCE_vacancy_irradiation_ID) sourceType
call source_vacancy_irradiation_deltaState(ipc, ip, el)
case (SOURCE_vacancy_thermalfluc_ID) sourceType
call source_vacancy_thermalfluc_deltaState(ipc, ip, el)
end select sourceType
enddo SourceLoop
end subroutine constitutive_collectDeltaState
@ -976,7 +1027,7 @@ end subroutine constitutive_collectDeltaState
!--------------------------------------------------------------------------------------------------
!> @brief returns array of constitutive results
!--------------------------------------------------------------------------------------------------
function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el)
function constitutive_postResults(S6, FeArray, ipc, ip, el)
use prec, only: &
pReal
use mesh, only: &
@ -1036,7 +1087,7 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el)
real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
FeArray !< elastic deformation gradient
real(pReal), intent(in), dimension(6) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel)
S6 !< 2nd Piola Kirchhoff stress (vector notation)
integer(pInt) :: &
startPos, endPos
integer(pInt) :: &
@ -1054,22 +1105,22 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el)
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_ISOTROPIC_ID) plasticityType
constitutive_postResults(startPos:endPos) = plastic_isotropic_postResults(Tstar_v,ipc,ip,el)
constitutive_postResults(startPos:endPos) = plastic_isotropic_postResults(S6,ipc,ip,el)
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
constitutive_postResults(startPos:endPos) = &
plastic_phenopowerlaw_postResults(Tstar_v,ipc,ip,el)
plastic_phenopowerlaw_postResults(S6,ipc,ip,el)
case (PLASTICITY_KINEHARDENING_ID) plasticityType
constitutive_postResults(startPos:endPos) = &
plastic_kinehardening_postResults(Tstar_v,ipc,ip,el)
plastic_kinehardening_postResults(S6,ipc,ip,el)
case (PLASTICITY_DISLOTWIN_ID) plasticityType
constitutive_postResults(startPos:endPos) = &
plastic_dislotwin_postResults(Tstar_v,temperature(ho)%p(tme),ipc,ip,el)
plastic_dislotwin_postResults(S6,temperature(ho)%p(tme),ipc,ip,el)
case (PLASTICITY_DISLOUCLA_ID) plasticityType
constitutive_postResults(startPos:endPos) = &
plastic_disloucla_postResults(Tstar_v,temperature(ho)%p(tme),ipc,ip,el)
plastic_disloucla_postResults(S6,temperature(ho)%p(tme),ipc,ip,el)
case (PLASTICITY_NONLOCAL_ID) plasticityType
constitutive_postResults(startPos:endPos) = &
plastic_nonlocal_postResults (Tstar_v,FeArray,ip,el)
plastic_nonlocal_postResults (S6,FeArray,ip,el)
end select plasticityType
SourceLoop: do s = 1_pInt, phase_Nsources(material_phase(ipc,ip,el))

View File

@ -570,9 +570,9 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
phase_Nsources, &
phaseAt, phasememberAt
use constitutive, only: &
constitutive_TandItsTangent, &
constitutive_LpAndItsTangent, &
constitutive_LiAndItsTangent
constitutive_SandItsTangents, &
constitutive_LpAndItsTangents, &
constitutive_LiAndItsTangents
implicit none
logical, intent(in) :: &
@ -1077,7 +1077,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
invFp = math_inv33(crystallite_partionedFp0(1:3,1:3,c,i,e))
Fe_guess = math_mul33x33(math_mul33x33(crystallite_partionedF(1:3,1:3,c,i,e), invFp), &
math_inv33(crystallite_partionedFi0(1:3,1:3,c,i,e)))
call constitutive_TandItsTangent(Tstar,dSdFe,dSdFi,Fe_guess,crystallite_partionedFi0(1:3,1:3,c,i,e),c,i,e)
call constitutive_SandItsTangents(Tstar,dSdFe,dSdFi,Fe_guess,crystallite_partionedFi0(1:3,1:3,c,i,e),c,i,e)
crystallite_P(1:3,1:3,c,i,e) = math_mul33x33(math_mul33x33(crystallite_partionedF(1:3,1:3,c,i,e), invFp), &
math_mul33x33(Tstar,transpose(invFp)))
endif
@ -1111,10 +1111,10 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
myNcomponents = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
do c = 1_pInt,myNcomponents
call constitutive_TandItsTangent(temp_33,dSdFe,dSdFi,crystallite_Fe(1:3,1:3,c,i,e), &
call constitutive_SandItsTangents(temp_33,dSdFe,dSdFi,crystallite_Fe(1:3,1:3,c,i,e), &
crystallite_Fi(1:3,1:3,c,i,e),c,i,e) ! call constitutive law to calculate elastic stress tangent
call constitutive_LiAndItsTangent(temp_33,dLidS,dLidFi,crystallite_Tstar_v(1:6,c,i,e), &
call constitutive_LiAndItsTangents(temp_33,dLidS,dLidFi,crystallite_Tstar_v(1:6,c,i,e), &
crystallite_Fi(1:3,1:3,c,i,e), &
c,i,e) ! call constitutive law to calculate Li tangent in lattice configuration
if (sum(abs(dLidS)) < tol_math_check) then
@ -1141,7 +1141,7 @@ subroutine crystallite_stressAndItsTangent(updateJaco)
dLidS = math_mul3333xx3333(dLidFi,dFidS) + dLidS
endif
call constitutive_LpAndItsTangent(temp_33,dLpdS,dLpdFi,crystallite_Tstar_v(1:6,c,i,e), &
call constitutive_LpAndItsTangents(temp_33,dLpdS,dLpdFi,crystallite_Tstar_v(1:6,c,i,e), &
crystallite_Fi(1:3,1:3,c,i,e),c,i,e) ! call constitutive law to calculate Lp tangent in lattice configuration
dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS
@ -1312,6 +1312,7 @@ subroutine crystallite_integrateStateRK4()
if (crystallite_todo(g,i,e)) &
call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), &
crystallite_Fe, &
crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_Fp, &
crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e)
enddo; enddo; enddo
@ -1454,6 +1455,7 @@ subroutine crystallite_integrateStateRK4()
if (crystallite_todo(g,i,e)) &
call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), &
crystallite_Fe, &
crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_Fp, &
timeStepFraction(n)*crystallite_subdt(g,i,e), & ! fraction of original timestep
crystallite_subFrac, g,i,e)
@ -1621,6 +1623,7 @@ subroutine crystallite_integrateStateRKCK45()
if (crystallite_todo(g,i,e)) &
call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), &
crystallite_Fe, &
crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_Fp, &
crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e)
enddo; enddo; enddo
@ -1773,6 +1776,7 @@ subroutine crystallite_integrateStateRKCK45()
if (crystallite_todo(g,i,e)) &
call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), &
crystallite_Fe, &
crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_Fp, &
C(stage)*crystallite_subdt(g,i,e), & ! fraction of original timestep
crystallite_subFrac, g,i,e)
@ -2098,6 +2102,7 @@ subroutine crystallite_integrateStateAdaptiveEuler()
if (crystallite_todo(g,i,e)) &
call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), &
crystallite_Fe, &
crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_Fp, &
crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e)
enddo; enddo; enddo
@ -2215,6 +2220,7 @@ subroutine crystallite_integrateStateAdaptiveEuler()
if (crystallite_todo(g,i,e)) &
call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), &
crystallite_Fe, &
crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_Fp, &
crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e)
enddo; enddo; enddo
@ -2410,6 +2416,7 @@ eIter = FEsolving_execElem(1:2)
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) &
call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), &
crystallite_Fe, &
crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_Fp, &
crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e)
enddo; enddo; enddo
@ -2669,6 +2676,7 @@ subroutine crystallite_integrateStateFPI()
if (crystallite_todo(g,i,e)) &
call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), &
crystallite_Fe, &
crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_Fp, &
crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e)
enddo; enddo; enddo
@ -2791,6 +2799,7 @@ subroutine crystallite_integrateStateFPI()
if (crystallite_todo(g,i,e) .and. .not. crystallite_converged(g,i,e)) &
call constitutive_collectDotState(crystallite_Tstar_v(1:6,g,i,e), &
crystallite_Fe, &
crystallite_Fi(1:3,1:3,g,i,e), &
crystallite_Fp, &
crystallite_subdt(g,i,e), crystallite_subFrac, g,i,e)
enddo; enddo; enddo
@ -3049,7 +3058,10 @@ logical function crystallite_stateJump(ipc,ip,el)
c = phasememberAt(ipc,ip,el)
p = phaseAt(ipc,ip,el)
call constitutive_collectDeltaState(crystallite_Tstar_v(1:6,ipc,ip,el), crystallite_Fe(1:3,1:3,ipc,ip,el), ipc,ip,el)
call constitutive_collectDeltaState(crystallite_Tstar_v(1:6,ipc,ip,el), &
crystallite_Fe(1:3,1:3,ipc,ip,el), &
crystallite_Fi(1:3,1:3,ipc,ip,el), &
ipc,ip,el)
myOffsetPlasticDeltaState = plasticState(p)%offsetDeltaState
mySizePlasticDeltaState = plasticState(p)%sizeDeltaState
@ -3156,9 +3168,9 @@ logical function crystallite_integrateStress(&
debug_levelExtensive, &
debug_levelSelective
use constitutive, only: constitutive_LpAndItsTangent, &
constitutive_LiAndItsTangent, &
constitutive_TandItsTangent
use constitutive, only: constitutive_LpAndItsTangents, &
constitutive_LiAndItsTangents, &
constitutive_SandItsTangents
use math, only: math_mul33x33, &
math_mul33xx33, &
math_mul3333xx3333, &
@ -3220,15 +3232,15 @@ logical function crystallite_integrateStress(&
real(pReal), dimension(9,9) :: dRLp_dLp, & ! partial derivative of residuum (Jacobian for NEwton-Raphson scheme)
dRLp_dLp2, & ! working copy of dRdLp
dRLi_dLi ! partial derivative of residuumI (Jacobian for NEwton-Raphson scheme)
real(pReal), dimension(3,3,3,3):: dT_dFe3333, & ! partial derivative of 2nd Piola-Kirchhoff stress
dT_dFi3333, &
dFe_dLp3333, & ! partial derivative of elastic deformation gradient
dFe_dLi3333, &
dFi_dLi3333, &
dLp_dFi3333, &
dLi_dFi3333, &
dLp_dT3333, &
dLi_dT3333
real(pReal), dimension(3,3,3,3):: dS_dFe, & ! partial derivative of 2nd Piola-Kirchhoff stress
dS_dFi, &
dFe_dLp, & ! partial derivative of elastic deformation gradient
dFe_dLi, &
dFi_dLi, &
dLp_dFi, &
dLi_dFi, &
dLp_dS, &
dLi_dS
real(pReal) detInvFi, & ! determinant of InvFi
steplengthLp, &
steplengthLi, &
@ -3349,7 +3361,7 @@ logical function crystallite_integrateStress(&
B = math_I3 - dt*Lpguess
Fe = math_mul33x33(math_mul33x33(A,B), invFi_new) ! current elastic deformation tensor
call constitutive_TandItsTangent(Tstar, dT_dFe3333, dT_dFi3333, &
call constitutive_SandItsTangents(Tstar, 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
Tstar_v = math_Mandel33to6(Tstar)
@ -3366,7 +3378,7 @@ logical function crystallite_integrateStress(&
write(6,'(a,/,6(e20.10,1x))') '<< CRYST >> Tstar', Tstar_v
endif
#endif
call constitutive_LpAndItsTangent(Lp_constitutive, dLp_dT3333, dLp_dFi3333, &
call constitutive_LpAndItsTangents(Lp_constitutive, dLp_dS, dLp_dFi, &
Tstar_v, Fi_new, ipc, ip, el)
#ifdef DEBUG
@ -3417,18 +3429,18 @@ logical function crystallite_integrateStress(&
!* calculate Jacobian for correction term
if (mod(jacoCounterLp, iJacoLpresiduum) == 0_pInt) then
dFe_dLp3333 = 0.0_pReal
dFe_dLp = 0.0_pReal
forall(o=1_pInt:3_pInt,p=1_pInt:3_pInt) &
dFe_dLp3333(o,1:3,p,1:3) = A(o,p)*transpose(invFi_new) ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j)
dFe_dLp3333 = - dt * dFe_dLp3333
dFe_dLp(o,1:3,p,1:3) = A(o,p)*transpose(invFi_new) ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j)
dFe_dLp = - dt * dFe_dLp
dRLp_dLp = math_identity2nd(9_pInt) &
- math_Plain3333to99(math_mul3333xx3333(math_mul3333xx3333(dLp_dT3333,dT_dFe3333),dFe_dLp3333))
- math_Plain3333to99(math_mul3333xx3333(math_mul3333xx3333(dLp_dS,dS_dFe),dFe_dLp))
#ifdef DEBUG
if (iand(debug_level(debug_crystallite), debug_levelExtensive) /= 0_pInt &
.and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,'(a,/,9(12x,9(e12.4,1x)/))') '<< CRYST >> dLp_dT', math_Plain3333to99(dLp_dT3333)
write(6,'(a,1x,e20.10)') '<< CRYST >> dLp_dT norm', norm2(math_Plain3333to99(dLp_dT3333))
write(6,'(a,/,9(12x,9(e12.4,1x)/))') '<< CRYST >> dLp_dS', math_Plain3333to99(dLp_dS)
write(6,'(a,1x,e20.10)') '<< CRYST >> dLp_dS norm', norm2(math_Plain3333to99(dLp_dS))
write(6,'(a,/,9(12x,9(e12.4,1x)/))') '<< CRYST >> dRLp_dLp', dRLp_dLp - math_identity2nd(9_pInt)
write(6,'(a,1x,e20.10)') '<< CRYST >> dRLp_dLp norm', norm2(dRLp_dLp - math_identity2nd(9_pInt))
endif
@ -3446,9 +3458,9 @@ logical function crystallite_integrateStress(&
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,*)
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dR_dLp',transpose(dRLp_dLp)
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dFe_dLp',transpose(math_Plain3333to99(dFe_dLp3333))
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dT_dFe_constitutive',transpose(math_Plain3333to99(dT_dFe3333))
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dLp_dT_constitutive',transpose(math_Plain3333to99(dLp_dT3333))
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dFe_dLp',transpose(math_Plain3333to99(dFe_dLp))
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dS_dFe_constitutive',transpose(math_Plain3333to99(dS_dFe))
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dLp_dS_constitutive',transpose(math_Plain3333to99(dLp_dS))
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> A',transpose(A)
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> B',transpose(B)
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Lp_constitutive',transpose(Lp_constitutive)
@ -3468,7 +3480,7 @@ logical function crystallite_integrateStress(&
!* calculate intermediate velocity gradient and its tangent from constitutive law
call constitutive_LiAndItsTangent(Li_constitutive, dLi_dT3333, dLi_dFi3333, &
call constitutive_LiAndItsTangents(Li_constitutive, dLi_dS, dLi_dFi, &
Tstar_v, Fi_new, ipc, ip, el)
#ifdef DEBUG
@ -3503,19 +3515,19 @@ logical function crystallite_integrateStress(&
if (mod(jacoCounterLi, iJacoLpresiduum) == 0_pInt) then
temp_33 = math_mul33x33(math_mul33x33(A,B),invFi_current)
dFe_dLi3333 = 0.0_pReal
dFi_dLi3333 = 0.0_pReal
dFe_dLi = 0.0_pReal
dFi_dLi = 0.0_pReal
forall(o=1_pInt:3_pInt,p=1_pInt:3_pInt)
dFe_dLi3333(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_dLi3333(1:3,o,1:3,p) = -dt*math_I3(o,p)*invFi_current
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
end forall
forall(o=1_pInt:3_pInt,p=1_pInt:3_pInt) &
dFi_dLi3333(1:3,1:3,o,p) = math_mul33x33(math_mul33x33(Fi_new,dFi_dLi3333(1:3,1:3,o,p)),Fi_new)
dFi_dLi(1:3,1:3,o,p) = math_mul33x33(math_mul33x33(Fi_new,dFi_dLi(1:3,1:3,o,p)),Fi_new)
dRLi_dLi = math_identity2nd(9_pInt) &
- math_Plain3333to99(math_mul3333xx3333(dLi_dT3333, math_mul3333xx3333(dT_dFe3333, dFe_dLi3333) + &
math_mul3333xx3333(dT_dFi3333, dFi_dLi3333))) &
- math_Plain3333to99(math_mul3333xx3333(dLi_dFi3333, dFi_dLi3333))
- math_Plain3333to99(math_mul3333xx3333(dLi_dS, math_mul3333xx3333(dS_dFe, dFe_dLi) + &
math_mul3333xx3333(dS_dFi, dFi_dLi))) &
- math_Plain3333to99(math_mul3333xx3333(dLi_dFi, dFi_dLi))
work = math_plain33to9(residuumLi)
call dgesv(9,1,dRLi_dLi,9,ipiv,work,9,ierr) ! solve dRLi/dLp * delta Li = -res for delta Li
if (ierr /= 0_pInt) then
@ -3528,9 +3540,9 @@ logical function crystallite_integrateStress(&
.or. .not. iand(debug_level(debug_crystallite), debug_levelSelective) /= 0_pInt)) then
write(6,*)
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dR_dLi',transpose(dRLi_dLi)
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dFe_dLi',transpose(math_Plain3333to99(dFe_dLi3333))
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dT_dFi_constitutive',transpose(math_Plain3333to99(dT_dFi3333))
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dLi_dT_constitutive',transpose(math_Plain3333to99(dLi_dT3333))
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dFe_dLi',transpose(math_Plain3333to99(dFe_dLi))
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dS_dFi_constitutive',transpose(math_Plain3333to99(dS_dFi))
write(6,'(a,/,9(12x,9(e15.3,1x)/))') '<< CRYST >> dLi_dS_constitutive',transpose(math_Plain3333to99(dLi_dS))
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Li_constitutive',transpose(Li_constitutive)
write(6,'(a,/,3(12x,3(e20.7,1x)/))') '<< CRYST >> Liguess',transpose(Liguess)
endif