unified syntax

_v for vector representation of tensor conflicts with 3333 suffix for
4-th order tensors.
General idea: Mark symmetric second order tensors in vector notation with
'6' and fourth order tensors in second order matrix notation with '99'.

Append nothing for 'natural' representation, i.e. F and NOT F33, Lp_dTstar and NOT
Lp_dTstar3333
This commit is contained in:
Martin Diehl 2018-08-25 13:55:39 +02:00
parent 2ad13a08e8
commit d585deee7e
2 changed files with 47 additions and 46 deletions

@ -1 +1 @@
Subproject commit dfd67ea44ba88ee1e0a33266a3986c64137908cf Subproject commit 81fd7109fea8456b8eecaaef0eec041edcce7792

View File

@ -345,6 +345,7 @@ end subroutine constitutive_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns the homogenize elasticity matrix !> @brief returns the homogenize elasticity matrix
!> ToDo: homogenizedC66 would be more consistent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function constitutive_homogenizedC(ipc,ip,el) function constitutive_homogenizedC(ipc,ip,el)
use prec, only: & use prec, only: &
@ -429,7 +430,7 @@ end subroutine constitutive_microstructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient !> @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_LpAndItsTangent(Lp, dLp_dTstar, dLp_dFi, Tstar6, Fi, ipc, ip, el)
use prec, only: & use prec, only: &
pReal pReal
use math, only: & use math, only: &
@ -469,18 +470,18 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v
ip, & !< integration point ip, & !< integration point
el !< element el !< element
real(pReal), intent(in), dimension(6) :: & real(pReal), intent(in), dimension(6) :: &
Tstar_v !< 2nd Piola-Kirchhoff stress Tstar6 !< 2nd Piola-Kirchhoff stress
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
Fi !< intermediate deformation gradient Fi !< intermediate deformation gradient
real(pReal), intent(out), dimension(3,3) :: & real(pReal), intent(out), dimension(3,3) :: &
Lp !< plastic velocity gradient Lp !< plastic velocity gradient
real(pReal), intent(out), dimension(3,3,3,3) :: & real(pReal), intent(out), dimension(3,3,3,3) :: &
dLp_dTstar3333, & !< derivative of Lp with respect to Tstar (4th-order tensor) dLp_dTstar, & !< derivative of Lp with respect to Tstar (4th-order tensor)
dLp_dFi3333 !< derivative of Lp with respect to Fi (4th-order tensor) dLp_dFi !< derivative of Lp with respect to Fi (4th-order tensor)
real(pReal), dimension(6) :: & real(pReal), dimension(6) :: &
Mstar_v !< Mandel stress work conjugate with Lp Mstar6 !< Mandel stress work conjugate with Lp
real(pReal), dimension(9,9) :: & real(pReal), dimension(9,9) :: &
dLp_dMstar !< derivative of Lp with respect to Mstar (4th-order tensor) dLp_dMstar99 !< derivative of Lp with respect to Mstar (4th-order tensor)
real(pReal), dimension(3,3) :: & real(pReal), dimension(3,3) :: &
temp_33 temp_33
integer(pInt) :: & integer(pInt) :: &
@ -492,39 +493,39 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar3333, dLp_dFi3333, Tstar_v
ho = material_homog(ip,el) ho = material_homog(ip,el)
tme = thermalMapping(ho)%p(ip,el) tme = thermalMapping(ho)%p(ip,el)
Mstar_v = math_Mandel33to6(math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(Tstar_v))) Mstar6 = math_Mandel33to6(math_mul33x33(math_mul33x33(transpose(Fi),Fi),math_Mandel6to33(Tstar6)))
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_NONE_ID) plasticityType case (PLASTICITY_NONE_ID) plasticityType
Lp = 0.0_pReal Lp = 0.0_pReal
dLp_dMstar = 0.0_pReal dLp_dMstar99 = 0.0_pReal
case (PLASTICITY_ISOTROPIC_ID) plasticityType case (PLASTICITY_ISOTROPIC_ID) plasticityType
call plastic_isotropic_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v,ipc,ip,el) call plastic_isotropic_LpAndItsTangent (Lp,dLp_dMstar99,Mstar6,ipc,ip,el)
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
call plastic_phenopowerlaw_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v,ipc,ip,el) call plastic_phenopowerlaw_LpAndItsTangent (Lp,dLp_dMstar99,Mstar6,ipc,ip,el)
case (PLASTICITY_KINEHARDENING_ID) plasticityType case (PLASTICITY_KINEHARDENING_ID) plasticityType
call plastic_kinehardening_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v,ipc,ip,el) call plastic_kinehardening_LpAndItsTangent (Lp,dLp_dMstar99,Mstar6,ipc,ip,el)
case (PLASTICITY_NONLOCAL_ID) plasticityType case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v, & call plastic_nonlocal_LpAndItsTangent (Lp,dLp_dMstar99,Mstar6, &
temperature(ho)%p(tme),ip,el) temperature(ho)%p(tme),ip,el)
case (PLASTICITY_DISLOTWIN_ID) plasticityType case (PLASTICITY_DISLOTWIN_ID) plasticityType
call plastic_dislotwin_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v, & call plastic_dislotwin_LpAndItsTangent (Lp,dLp_dMstar99,Mstar6, &
temperature(ho)%p(tme),ipc,ip,el) temperature(ho)%p(tme),ipc,ip,el)
case (PLASTICITY_DISLOUCLA_ID) plasticityType case (PLASTICITY_DISLOUCLA_ID) plasticityType
call plastic_disloucla_LpAndItsTangent (Lp,dLp_dMstar,Mstar_v, & call plastic_disloucla_LpAndItsTangent (Lp,dLp_dMstar99,Mstar6, &
temperature(ho)%p(tme), ipc,ip,el) temperature(ho)%p(tme), ipc,ip,el)
end select plasticityType end select plasticityType
dLp_dTstar3333 = math_Plain99to3333(dLp_dMstar) dLp_dTstar = math_Plain99to3333(dLp_dMstar99)
temp_33 = math_mul33x33(Fi,math_Mandel6to33(Tstar_v)) temp_33 = math_mul33x33(Fi,math_Mandel6to33(Tstar6))
forall(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt) & 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))) + & dLp_dFi(i,j,1:3,1:3) = math_mul33x33(temp_33,transpose(dLp_dTstar(i,j,1:3,1:3))) + &
math_mul33x33(math_mul33x33(Fi,dLp_dTstar3333(i,j,1:3,1:3)),math_Mandel6to33(Tstar_v)) math_mul33x33(math_mul33x33(Fi,dLp_dTstar(i,j,1:3,1:3)),math_Mandel6to33(Tstar6))
temp_33 = math_mul33x33(transpose(Fi),Fi) temp_33 = math_mul33x33(transpose(Fi),Fi)
forall(i = 1_pInt:3_pInt, j = 1_pInt:3_pInt) & 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)) dLp_dTstar(i,j,1:3,1:3) = math_mul33x33(temp_33,dLp_dTstar(i,j,1:3,1:3))
end subroutine constitutive_LpAndItsTangent end subroutine constitutive_LpAndItsTangent
@ -532,7 +533,7 @@ end subroutine constitutive_LpAndItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient !> @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_LiAndItsTangent(Li, dLi_dTstar3333, dLi_dFi3333, Tstar6, Fi, ipc, ip, el)
use prec, only: & use prec, only: &
pReal pReal
use math, only: & use math, only: &
@ -570,7 +571,7 @@ subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar3333, dLi_dFi3333, Tstar_v
ip, & !< integration point ip, & !< integration point
el !< element el !< element
real(pReal), intent(in), dimension(6) :: & real(pReal), intent(in), dimension(6) :: &
Tstar_v !< 2nd Piola-Kirchhoff stress Tstar6 !< 2nd Piola-Kirchhoff stress
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
Fi !< intermediate deformation gradient Fi !< intermediate deformation gradient
real(pReal), intent(out), dimension(3,3) :: & real(pReal), intent(out), dimension(3,3) :: &
@ -598,7 +599,7 @@ subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar3333, dLi_dFi3333, Tstar_v
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_isotropic_ID) plasticityType 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_dTstar, Tstar6, ipc, ip, el)
case default plasticityType case default plasticityType
my_Li = 0.0_pReal my_Li = 0.0_pReal
my_dLi_dTstar = 0.0_pReal my_dLi_dTstar = 0.0_pReal
@ -610,9 +611,9 @@ subroutine constitutive_LiAndItsTangent(Li, dLi_dTstar3333, dLi_dFi3333, Tstar_v
KinematicsLoop: do k = 1_pInt, phase_Nkinematics(material_phase(ipc,ip,el)) KinematicsLoop: do k = 1_pInt, phase_Nkinematics(material_phase(ipc,ip,el))
kinematicsType: select case (phase_kinematics(k,material_phase(ipc,ip,el))) kinematicsType: select case (phase_kinematics(k,material_phase(ipc,ip,el)))
case (KINEMATICS_cleavage_opening_ID) kinematicsType 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_dTstar, Tstar6, ipc, ip, el)
case (KINEMATICS_slipplane_opening_ID) kinematicsType 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_dTstar, Tstar6, ipc, ip, el)
case (KINEMATICS_thermal_expansion_ID) kinematicsType 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_dTstar, ipc, ip, el)
case (KINEMATICS_vacancy_strain_ID) kinematicsType case (KINEMATICS_vacancy_strain_ID) kinematicsType
@ -798,7 +799,7 @@ end subroutine constitutive_hooke_TandItsTangent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure !> @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(Tstar6, FeArray, FpArray, subdt, subfracArray,ipc, ip, el)
use prec, only: & use prec, only: &
pReal, & pReal, &
pLongInt pLongInt
@ -863,7 +864,7 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra
FeArray, & !< elastic deformation gradient FeArray, & !< elastic deformation gradient
FpArray !< plastic deformation gradient FpArray !< plastic deformation gradient
real(pReal), intent(in), dimension(6) :: & real(pReal), intent(in), dimension(6) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel) Tstar6 !< 2nd Piola Kirchhoff stress tensor (Mandel)
integer(pInt) :: & integer(pInt) :: &
ho, & !< homogenization ho, & !< homogenization
tme, & !< thermal member position tme, & !< thermal member position
@ -874,26 +875,26 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, subdt, subfra
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_ISOTROPIC_ID) plasticityType case (PLASTICITY_ISOTROPIC_ID) plasticityType
call plastic_isotropic_dotState (Tstar_v,ipc,ip,el) call plastic_isotropic_dotState (Tstar6,ipc,ip,el)
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
call plastic_phenopowerlaw_dotState(Tstar_v,ipc,ip,el) call plastic_phenopowerlaw_dotState(Tstar6,ipc,ip,el)
case (PLASTICITY_KINEHARDENING_ID) plasticityType case (PLASTICITY_KINEHARDENING_ID) plasticityType
call plastic_kinehardening_dotState(Tstar_v,ipc,ip,el) call plastic_kinehardening_dotState(Tstar6,ipc,ip,el)
case (PLASTICITY_DISLOTWIN_ID) plasticityType case (PLASTICITY_DISLOTWIN_ID) plasticityType
call plastic_dislotwin_dotState (Tstar_v,temperature(ho)%p(tme), & call plastic_dislotwin_dotState (Tstar6,temperature(ho)%p(tme), &
ipc,ip,el) ipc,ip,el)
case (PLASTICITY_DISLOUCLA_ID) plasticityType case (PLASTICITY_DISLOUCLA_ID) plasticityType
call plastic_disloucla_dotState (Tstar_v,temperature(ho)%p(tme), & call plastic_disloucla_dotState (Tstar6,temperature(ho)%p(tme), &
ipc,ip,el) ipc,ip,el)
case (PLASTICITY_NONLOCAL_ID) plasticityType case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_dotState (Tstar_v,FeArray,FpArray,temperature(ho)%p(tme), & call plastic_nonlocal_dotState (Tstar6,FeArray,FpArray,temperature(ho)%p(tme), &
subdt,subfracArray,ip,el) subdt,subfracArray,ip,el)
end select plasticityType 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))) sourceType: select case (phase_source(s,material_phase(ipc,ip,el)))
case (SOURCE_damage_anisoBrittle_ID) sourceType case (SOURCE_damage_anisoBrittle_ID) sourceType
call source_damage_anisoBrittle_dotState (Tstar_v, ipc, ip, el) call source_damage_anisoBrittle_dotState (Tstar6, ipc, ip, el)
case (SOURCE_damage_isoDuctile_ID) sourceType case (SOURCE_damage_isoDuctile_ID) sourceType
call source_damage_isoDuctile_dotState ( ipc, ip, el) call source_damage_isoDuctile_dotState ( ipc, ip, el)
case (SOURCE_damage_anisoDuctile_ID) sourceType case (SOURCE_damage_anisoDuctile_ID) sourceType
@ -909,7 +910,7 @@ end subroutine constitutive_collectDotState
!> @brief for constitutive models having an instantaneous change of state !> @brief for constitutive models having an instantaneous change of state
!> will return false if delta state is not needed/supported by the constitutive model !> 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(Tstar6, Fe, ipc, ip, el)
use prec, only: & use prec, only: &
pReal, & pReal, &
pLongInt pLongInt
@ -944,7 +945,7 @@ subroutine constitutive_collectDeltaState(Tstar_v, Fe, ipc, ip, el)
ip, & !< integration point ip, & !< integration point
el !< element el !< element
real(pReal), intent(in), dimension(6) :: & real(pReal), intent(in), dimension(6) :: &
Tstar_v !< 2nd Piola-Kirchhoff stress Tstar6 !< 2nd Piola-Kirchhoff stress
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
Fe !< elastic deformation gradient Fe !< elastic deformation gradient
integer(pInt) :: & integer(pInt) :: &
@ -952,9 +953,9 @@ subroutine constitutive_collectDeltaState(Tstar_v, Fe, ipc, ip, el)
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_KINEHARDENING_ID) plasticityType case (PLASTICITY_KINEHARDENING_ID) plasticityType
call plastic_kinehardening_deltaState(Tstar_v,ipc,ip,el) call plastic_kinehardening_deltaState(Tstar6,ipc,ip,el)
case (PLASTICITY_NONLOCAL_ID) plasticityType case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_deltaState(Tstar_v,ip,el) call plastic_nonlocal_deltaState(Tstar6,ip,el)
end select plasticityType 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))
@ -975,7 +976,7 @@ end subroutine constitutive_collectDeltaState
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief returns array of constitutive results !> @brief returns array of constitutive results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el) function constitutive_postResults(Tstar6, FeArray, ipc, ip, el)
use prec, only: & use prec, only: &
pReal pReal
use mesh, only: & use mesh, only: &
@ -1035,7 +1036,7 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el)
real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: & real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems) :: &
FeArray !< elastic deformation gradient FeArray !< elastic deformation gradient
real(pReal), intent(in), dimension(6) :: & real(pReal), intent(in), dimension(6) :: &
Tstar_v !< 2nd Piola Kirchhoff stress tensor (Mandel) Tstar6 !< 2nd Piola Kirchhoff stress tensor (Mandel)
integer(pInt) :: & integer(pInt) :: &
startPos, endPos startPos, endPos
integer(pInt) :: & integer(pInt) :: &
@ -1053,22 +1054,22 @@ function constitutive_postResults(Tstar_v, FeArray, ipc, ip, el)
plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el))) plasticityType: select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_ISOTROPIC_ID) plasticityType case (PLASTICITY_ISOTROPIC_ID) plasticityType
constitutive_postResults(startPos:endPos) = plastic_isotropic_postResults(Tstar_v,ipc,ip,el) constitutive_postResults(startPos:endPos) = plastic_isotropic_postResults(Tstar6,ipc,ip,el)
case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType case (PLASTICITY_PHENOPOWERLAW_ID) plasticityType
constitutive_postResults(startPos:endPos) = & constitutive_postResults(startPos:endPos) = &
plastic_phenopowerlaw_postResults(Tstar_v,ipc,ip,el) plastic_phenopowerlaw_postResults(Tstar6,ipc,ip,el)
case (PLASTICITY_KINEHARDENING_ID) plasticityType case (PLASTICITY_KINEHARDENING_ID) plasticityType
constitutive_postResults(startPos:endPos) = & constitutive_postResults(startPos:endPos) = &
plastic_kinehardening_postResults(Tstar_v,ipc,ip,el) plastic_kinehardening_postResults(Tstar6,ipc,ip,el)
case (PLASTICITY_DISLOTWIN_ID) plasticityType case (PLASTICITY_DISLOTWIN_ID) plasticityType
constitutive_postResults(startPos:endPos) = & constitutive_postResults(startPos:endPos) = &
plastic_dislotwin_postResults(Tstar_v,temperature(ho)%p(tme),ipc,ip,el) plastic_dislotwin_postResults(Tstar6,temperature(ho)%p(tme),ipc,ip,el)
case (PLASTICITY_DISLOUCLA_ID) plasticityType case (PLASTICITY_DISLOUCLA_ID) plasticityType
constitutive_postResults(startPos:endPos) = & constitutive_postResults(startPos:endPos) = &
plastic_disloucla_postResults(Tstar_v,temperature(ho)%p(tme),ipc,ip,el) plastic_disloucla_postResults(Tstar6,temperature(ho)%p(tme),ipc,ip,el)
case (PLASTICITY_NONLOCAL_ID) plasticityType case (PLASTICITY_NONLOCAL_ID) plasticityType
constitutive_postResults(startPos:endPos) = & constitutive_postResults(startPos:endPos) = &
plastic_nonlocal_postResults (Tstar_v,FeArray,ip,el) plastic_nonlocal_postResults (Tstar6,FeArray,ip,el)
end select plasticityType 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))