new names

This commit is contained in:
Martin Diehl 2021-04-06 11:38:44 +02:00
parent 2b798589ce
commit 869976c7a0
22 changed files with 174 additions and 177 deletions

View File

@ -259,25 +259,25 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
NiterationMPstate, & NiterationMPstate, &
ip, & !< integration point number ip, & !< integration point number
el, & !< element number el, & !< element number
myNgrains, co, ce, ho, me, ph myNgrains, co, ce, ho, en, ph
logical :: & logical :: &
converged converged
logical, dimension(2) :: & logical, dimension(2) :: &
doneAndHappy doneAndHappy
!$OMP PARALLEL !$OMP PARALLEL
!$OMP DO PRIVATE(ce,me,ho,myNgrains,NiterationMPstate,converged,doneAndHappy) !$OMP DO PRIVATE(ce,en,ho,myNgrains,NiterationMPstate,converged,doneAndHappy)
do el = FEsolving_execElem(1),FEsolving_execElem(2) do el = FEsolving_execElem(1),FEsolving_execElem(2)
ho = material_homogenizationAt(el) ho = material_homogenizationAt(el)
myNgrains = homogenization_Nconstituents(ho) myNgrains = homogenization_Nconstituents(ho)
do ip = FEsolving_execIP(1),FEsolving_execIP(2) do ip = FEsolving_execIP(1),FEsolving_execIP(2)
ce = (el-1)*discretization_nIPs + ip ce = (el-1)*discretization_nIPs + ip
me = material_homogenizationMemberAt2(ce) en = material_homogenizationEntry(ce)
call phase_restore(ce,.false.) ! wrong name (is more a forward function) call phase_restore(ce,.false.) ! wrong name (is more a forward function)
if(homogState(ho)%sizeState > 0) homogState(ho)%state(:,me) = homogState(ho)%state0(:,me) if(homogState(ho)%sizeState > 0) homogState(ho)%state(:,en) = homogState(ho)%state0(:,en)
if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%state(:,me) = damageState_h(ho)%state0(:,me) if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%state(:,en) = damageState_h(ho)%state0(:,en)
call damage_partition(ce) call damage_partition(ce)
doneAndHappy = [.false.,.true.] doneAndHappy = [.false.,.true.]
@ -505,12 +505,12 @@ function damage_nonlocal_getDiffusion(ce)
ho, & ho, &
co co
ho = material_homogenizationAt2(ce) ho = material_homogenizationID(ce)
damage_nonlocal_getDiffusion = 0.0_pReal damage_nonlocal_getDiffusion = 0.0_pReal
do co = 1, homogenization_Nconstituents(ho) do co = 1, homogenization_Nconstituents(ho)
damage_nonlocal_getDiffusion = damage_nonlocal_getDiffusion + & damage_nonlocal_getDiffusion = damage_nonlocal_getDiffusion + &
crystallite_push33ToRef(co,ce,lattice_D(1:3,1:3,material_phaseAt2(co,ce))) crystallite_push33ToRef(co,ce,lattice_D(1:3,1:3,material_phaseID(co,ce)))
enddo enddo
damage_nonlocal_getDiffusion = & damage_nonlocal_getDiffusion = &

View File

@ -42,7 +42,7 @@ module subroutine damage_init()
allocate(current(configHomogenizations%length)) allocate(current(configHomogenizations%length))
do ho = 1, configHomogenizations%length do ho = 1, configHomogenizations%length
allocate(current(ho)%phi(count(material_homogenizationAt2==ho)), source=1.0_pReal) allocate(current(ho)%phi(count(material_homogenizationID==ho)), source=1.0_pReal)
configHomogenization => configHomogenizations%get(ho) configHomogenization => configHomogenizations%get(ho)
associate(prm => param(ho)) associate(prm => param(ho))
if (configHomogenization%contains('damage')) then if (configHomogenization%contains('damage')) then
@ -72,9 +72,9 @@ module subroutine damage_partition(ce)
integer :: co integer :: co
if(damageState_h(material_homogenizationAt2(ce))%sizeState < 1) return if(damageState_h(material_homogenizationID(ce))%sizeState < 1) return
phi = damagestate_h(material_homogenizationAt2(ce))%state(1,material_homogenizationMemberAt2(ce)) phi = damagestate_h(material_homogenizationID(ce))%state(1,material_homogenizationEntry(ce))
do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) do co = 1, homogenization_Nconstituents(material_homogenizationID(ce))
call phase_damage_set_phi(phi,co,ce) call phase_damage_set_phi(phi,co,ce)
enddo enddo
@ -94,11 +94,11 @@ module function damage_nonlocal_getMobility(ce) result(M)
M = 0.0_pReal M = 0.0_pReal
do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) do co = 1, homogenization_Nconstituents(material_homogenizationID(ce))
M = M + lattice_M(material_phaseAt2(co,ce)) M = M + lattice_M(material_phaseID(co,ce))
enddo enddo
M = M/real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal) M = M/real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
end function damage_nonlocal_getMobility end function damage_nonlocal_getMobility
@ -118,8 +118,8 @@ module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, p
dPhiDot_dPhi = 0.0_pReal dPhiDot_dPhi = 0.0_pReal
call phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ce) call phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ce)
phiDot = phiDot/real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal) phiDot = phiDot/real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal) dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
end subroutine damage_nonlocal_getSourceAndItsTangent end subroutine damage_nonlocal_getSourceAndItsTangent
@ -134,11 +134,11 @@ module subroutine damage_nonlocal_putNonLocalDamage(phi,ce)
phi phi
integer :: & integer :: &
ho, & ho, &
me en
ho = material_homogenizationAt2(ce) ho = material_homogenizationID(ce)
me = material_homogenizationMemberAt2(ce) en = material_homogenizationEntry(ce)
damagestate_h(ho)%state(1,me) = phi damagestate_h(ho)%state(1,en) = phi
end subroutine damage_nonlocal_putNonLocalDamage end subroutine damage_nonlocal_putNonLocalDamage

View File

@ -110,10 +110,10 @@ module subroutine mechanical_partition(subF,ce)
ce ce
integer :: co integer :: co
real(pReal), dimension (3,3,homogenization_Nconstituents(material_homogenizationAt2(ce))) :: Fs real(pReal), dimension (3,3,homogenization_Nconstituents(material_homogenizationID(ce))) :: Fs
chosenHomogenization: select case(homogenization_type(material_homogenizationAt2(ce))) chosenHomogenization: select case(homogenization_type(material_homogenizationID(ce)))
case (HOMOGENIZATION_NONE_ID) chosenHomogenization case (HOMOGENIZATION_NONE_ID) chosenHomogenization
Fs(1:3,1:3,1) = subF Fs(1:3,1:3,1) = subF
@ -126,7 +126,7 @@ module subroutine mechanical_partition(subF,ce)
end select chosenHomogenization end select chosenHomogenization
do co = 1,homogenization_Nconstituents(material_homogenizationAt2(ce)) do co = 1,homogenization_Nconstituents(material_homogenizationID(ce))
call phase_mechanical_setF(Fs(1:3,1:3,co),co,ce) call phase_mechanical_setF(Fs(1:3,1:3,co),co,ce)
enddo enddo
@ -143,18 +143,18 @@ module subroutine mechanical_homogenize(dt,ce)
integer, intent(in) :: ce integer, intent(in) :: ce
integer :: co integer :: co
real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt2(ce))) real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationID(ce)))
real(pReal) :: Ps(3,3,homogenization_Nconstituents(material_homogenizationAt2(ce))) real(pReal) :: Ps(3,3,homogenization_Nconstituents(material_homogenizationID(ce)))
chosenHomogenization: select case(homogenization_type(material_homogenizationAt2(ce))) chosenHomogenization: select case(homogenization_type(material_homogenizationID(ce)))
case (HOMOGENIZATION_NONE_ID) chosenHomogenization case (HOMOGENIZATION_NONE_ID) chosenHomogenization
homogenization_P(1:3,1:3,ce) = phase_mechanical_getP(1,ce) homogenization_P(1:3,1:3,ce) = phase_mechanical_getP(1,ce)
homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = phase_mechanical_dPdF(dt,1,ce) homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = phase_mechanical_dPdF(dt,1,ce)
case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) do co = 1, homogenization_Nconstituents(material_homogenizationID(ce))
dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(dt,co,ce) dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(dt,co,ce)
Ps(:,:,co) = phase_mechanical_getP(co,ce) Ps(:,:,co) = phase_mechanical_getP(co,ce)
enddo enddo
@ -162,10 +162,10 @@ module subroutine mechanical_homogenize(dt,ce)
homogenization_P(1:3,1:3,ce), & homogenization_P(1:3,1:3,ce), &
homogenization_dPdF(1:3,1:3,1:3,1:3,ce),& homogenization_dPdF(1:3,1:3,1:3,1:3,ce),&
Ps,dPdFs, & Ps,dPdFs, &
material_homogenizationAt2(ce)) material_homogenizationID(ce))
case (HOMOGENIZATION_RGC_ID) chosenHomogenization case (HOMOGENIZATION_RGC_ID) chosenHomogenization
do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) do co = 1, homogenization_Nconstituents(material_homogenizationID(ce))
dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(dt,co,ce) dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(dt,co,ce)
Ps(:,:,co) = phase_mechanical_getP(co,ce) Ps(:,:,co) = phase_mechanical_getP(co,ce)
enddo enddo
@ -173,7 +173,7 @@ module subroutine mechanical_homogenize(dt,ce)
homogenization_P(1:3,1:3,ce), & homogenization_P(1:3,1:3,ce), &
homogenization_dPdF(1:3,1:3,1:3,1:3,ce),& homogenization_dPdF(1:3,1:3,1:3,1:3,ce),&
Ps,dPdFs, & Ps,dPdFs, &
material_homogenizationAt2(ce)) material_homogenizationID(ce))
end select chosenHomogenization end select chosenHomogenization
@ -195,13 +195,13 @@ module function mechanical_updateState(subdt,subF,ce) result(doneAndHappy)
logical, dimension(2) :: doneAndHappy logical, dimension(2) :: doneAndHappy
integer :: co integer :: co
real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt2(ce))) real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationID(ce)))
real(pReal) :: Fs(3,3,homogenization_Nconstituents(material_homogenizationAt2(ce))) real(pReal) :: Fs(3,3,homogenization_Nconstituents(material_homogenizationID(ce)))
real(pReal) :: Ps(3,3,homogenization_Nconstituents(material_homogenizationAt2(ce))) real(pReal) :: Ps(3,3,homogenization_Nconstituents(material_homogenizationID(ce)))
if (homogenization_type(material_homogenizationAt2(ce)) == HOMOGENIZATION_RGC_ID) then if (homogenization_type(material_homogenizationID(ce)) == HOMOGENIZATION_RGC_ID) then
do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) do co = 1, homogenization_Nconstituents(material_homogenizationID(ce))
dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(subdt,co,ce) dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(subdt,co,ce)
Fs(:,:,co) = phase_mechanical_getF(co,ce) Fs(:,:,co) = phase_mechanical_getF(co,ce)
Ps(:,:,co) = phase_mechanical_getP(co,ce) Ps(:,:,co) = phase_mechanical_getP(co,ce)

View File

@ -204,12 +204,12 @@ module subroutine mechanical_RGC_partitionDeformation(F,avgF,ce)
real(pReal), dimension(3) :: aVect,nVect real(pReal), dimension(3) :: aVect,nVect
integer, dimension(4) :: intFace integer, dimension(4) :: intFace
integer, dimension(3) :: iGrain3 integer, dimension(3) :: iGrain3
integer :: iGrain,iFace,i,j,ho,me integer :: iGrain,iFace,i,j,ho,en
associate(prm => param(material_homogenizationAt2(ce))) associate(prm => param(material_homogenizationID(ce)))
ho = material_homogenizationAt2(ce) ho = material_homogenizationID(ce)
me = material_homogenizationMemberAt2(ce) en = material_homogenizationEntry(ce)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! compute the deformation gradient of individual grains due to relaxations ! compute the deformation gradient of individual grains due to relaxations
F = 0.0_pReal F = 0.0_pReal
@ -217,8 +217,8 @@ module subroutine mechanical_RGC_partitionDeformation(F,avgF,ce)
iGrain3 = grain1to3(iGrain,prm%N_constituents) iGrain3 = grain1to3(iGrain,prm%N_constituents)
do iFace = 1,6 do iFace = 1,6
intFace = getInterface(iFace,iGrain3) ! identifying 6 interfaces of each grain intFace = getInterface(iFace,iGrain3) ! identifying 6 interfaces of each grain
aVect = relaxationVector(intFace,ho,me) ! get the relaxation vectors for each interface from global relaxation vector array aVect = relaxationVector(intFace,ho,en) ! get the relaxation vectors for each interface from global relaxation vector array
nVect = interfaceNormal(intFace,ho,me) nVect = interfaceNormal(intFace,ho,en)
forall (i=1:3,j=1:3) & forall (i=1:3,j=1:3) &
F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! calculating deformation relaxations due to interface relaxation F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! calculating deformation relaxations due to interface relaxation
enddo enddo
@ -247,7 +247,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
integer, dimension(4) :: intFaceN,intFaceP,faceID integer, dimension(4) :: intFaceN,intFaceP,faceID
integer, dimension(3) :: nGDim,iGr3N,iGr3P integer, dimension(3) :: nGDim,iGr3N,iGr3P
integer :: ho,iNum,i,j,nIntFaceTot,iGrN,iGrP,iMun,iFace,k,l,ipert,iGrain,nGrain, me integer :: ho,iNum,i,j,nIntFaceTot,iGrN,iGrP,iMun,iFace,k,l,ipert,nGrain, en
real(pReal), dimension(3,3,size(P,3)) :: R,pF,pR,D,pD real(pReal), dimension(3,3,size(P,3)) :: R,pF,pR,D,pD
real(pReal), dimension(3,size(P,3)) :: NN,devNull real(pReal), dimension(3,size(P,3)) :: NN,devNull
real(pReal), dimension(3) :: normP,normN,mornP,mornN real(pReal), dimension(3) :: normP,normN,mornP,mornN
@ -261,9 +261,9 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
return return
endif zeroTimeStep endif zeroTimeStep
ho = material_homogenizationAt2(ce) ho = material_homogenizationID(ce)
en = material_homogenizationEntry(ce)
me = material_homogenizationMemberAt2(ce)
associate(stt => state(ho), st0 => state0(ho), dst => dependentState(ho), prm => param(ho)) associate(stt => state(ho), st0 => state0(ho), dst => dependentState(ho), prm => param(ho))
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -278,16 +278,16 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
! allocate the size of the global relaxation arrays/jacobian matrices depending on the size of the cluster ! allocate the size of the global relaxation arrays/jacobian matrices depending on the size of the cluster
allocate(resid(3*nIntFaceTot), source=0.0_pReal) allocate(resid(3*nIntFaceTot), source=0.0_pReal)
allocate(tract(nIntFaceTot,3), source=0.0_pReal) allocate(tract(nIntFaceTot,3), source=0.0_pReal)
relax = stt%relaxationVector(:,me) relax = stt%relaxationVector(:,en)
drelax = stt%relaxationVector(:,me) - st0%relaxationVector(:,me) drelax = stt%relaxationVector(:,en) - st0%relaxationVector(:,en)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! computing interface mismatch and stress penalty tensor for all interfaces of all grains ! computing interface mismatch and stress penalty tensor for all interfaces of all grains
call stressPenalty(R,NN,avgF,F,ho,me) call stressPenalty(R,NN,avgF,F,ho,en)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! calculating volume discrepancy and stress penalty related to overall volume discrepancy ! calculating volume discrepancy and stress penalty related to overall volume discrepancy
call volumePenalty(D,dst%volumeDiscrepancy(me),avgF,F,nGrain) call volumePenalty(D,dst%volumeDiscrepancy(en),avgF,F,nGrain)
!------------------------------------------------------------------------------------------------ !------------------------------------------------------------------------------------------------
! computing the residual stress from the balance of traction at all (interior) interfaces ! computing the residual stress from the balance of traction at all (interior) interfaces
@ -299,7 +299,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
iGr3N = faceID(2:4) ! identifying the grain ID in local coordinate system (3-dimensional index) iGr3N = faceID(2:4) ! identifying the grain ID in local coordinate system (3-dimensional index)
iGrN = grain3to1(iGr3N,param(ho)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index) iGrN = grain3to1(iGr3N,param(ho)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index)
intFaceN = getInterface(2*faceID(1),iGr3N) intFaceN = getInterface(2*faceID(1),iGr3N)
normN = interfaceNormal(intFaceN,ho,me) normN = interfaceNormal(intFaceN,ho,en)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! identify the right/up/front grain (+|P) ! identify the right/up/front grain (+|P)
@ -307,7 +307,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identifying the grain ID in local coordinate system (3-dimensional index) iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identifying the grain ID in local coordinate system (3-dimensional index)
iGrP = grain3to1(iGr3P,param(ho)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index) iGrP = grain3to1(iGr3P,param(ho)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index)
intFaceP = getInterface(2*faceID(1)-1,iGr3P) intFaceP = getInterface(2*faceID(1)-1,iGr3P)
normP = interfaceNormal(intFaceP,ho,me) normP = interfaceNormal(intFaceP,ho,en)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! compute the residual of traction at the interface (in local system, 4-dimensional index) ! compute the residual of traction at the interface (in local system, 4-dimensional index)
@ -335,9 +335,9 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
if (residMax < num%rtol*stresMax .or. residMax < num%atol) then if (residMax < num%rtol*stresMax .or. residMax < num%atol) then
doneAndHappy = .true. doneAndHappy = .true.
dst%mismatch(1:3,me) = sum(NN,2)/real(nGrain,pReal) dst%mismatch(1:3,en) = sum(NN,2)/real(nGrain,pReal)
dst%relaxationRate_avg(me) = sum(abs(drelax))/dt/real(3*nIntFaceTot,pReal) dst%relaxationRate_avg(en) = sum(abs(drelax))/dt/real(3*nIntFaceTot,pReal)
dst%relaxationRate_max(me) = maxval(abs(drelax))/dt dst%relaxationRate_max(en) = maxval(abs(drelax))/dt
return return
@ -363,10 +363,10 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
iGr3N = faceID(2:4) ! identifying the grain ID in local coordinate sytem iGr3N = faceID(2:4) ! identifying the grain ID in local coordinate sytem
iGrN = grain3to1(iGr3N,param(ho)%N_constituents) ! translate into global grain ID iGrN = grain3to1(iGr3N,param(ho)%N_constituents) ! translate into global grain ID
intFaceN = getInterface(2*faceID(1),iGr3N) ! identifying the connecting interface in local coordinate system intFaceN = getInterface(2*faceID(1),iGr3N) ! identifying the connecting interface in local coordinate system
normN = interfaceNormal(intFaceN,ho,me) normN = interfaceNormal(intFaceN,ho,en)
do iFace = 1,6 do iFace = 1,6
intFaceN = getInterface(iFace,iGr3N) ! identifying all interfaces that influence relaxation of the above interface intFaceN = getInterface(iFace,iGr3N) ! identifying all interfaces that influence relaxation of the above interface
mornN = interfaceNormal(intFaceN,ho,me) mornN = interfaceNormal(intFaceN,ho,en)
iMun = interface4to1(intFaceN,param(ho)%N_constituents) ! translate the interfaces ID into local 4-dimensional index iMun = interface4to1(intFaceN,param(ho)%N_constituents) ! translate the interfaces ID into local 4-dimensional index
if (iMun > 0) then ! get the corresponding tangent if (iMun > 0) then ! get the corresponding tangent
do i=1,3; do j=1,3; do k=1,3; do l=1,3 do i=1,3; do j=1,3; do k=1,3; do l=1,3
@ -384,10 +384,10 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identifying the grain ID in local coordinate sytem iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identifying the grain ID in local coordinate sytem
iGrP = grain3to1(iGr3P,param(ho)%N_constituents) ! translate into global grain ID iGrP = grain3to1(iGr3P,param(ho)%N_constituents) ! translate into global grain ID
intFaceP = getInterface(2*faceID(1)-1,iGr3P) ! identifying the connecting interface in local coordinate system intFaceP = getInterface(2*faceID(1)-1,iGr3P) ! identifying the connecting interface in local coordinate system
normP = interfaceNormal(intFaceP,ho,me) normP = interfaceNormal(intFaceP,ho,en)
do iFace = 1,6 do iFace = 1,6
intFaceP = getInterface(iFace,iGr3P) ! identifying all interfaces that influence relaxation of the above interface intFaceP = getInterface(iFace,iGr3P) ! identifying all interfaces that influence relaxation of the above interface
mornP = interfaceNormal(intFaceP,ho,me) mornP = interfaceNormal(intFaceP,ho,en)
iMun = interface4to1(intFaceP,param(ho)%N_constituents) ! translate the interfaces ID into local 4-dimensional index iMun = interface4to1(intFaceP,param(ho)%N_constituents) ! translate the interfaces ID into local 4-dimensional index
if (iMun > 0) then ! get the corresponding tangent if (iMun > 0) then ! get the corresponding tangent
do i=1,3; do j=1,3; do k=1,3; do l=1,3 do i=1,3; do j=1,3; do k=1,3; do l=1,3
@ -408,9 +408,9 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
do ipert = 1,3*nIntFaceTot do ipert = 1,3*nIntFaceTot
p_relax = relax p_relax = relax
p_relax(ipert) = relax(ipert) + num%pPert ! perturb the relaxation vector p_relax(ipert) = relax(ipert) + num%pPert ! perturb the relaxation vector
stt%relaxationVector(:,me) = p_relax stt%relaxationVector(:,en) = p_relax
call grainDeformation(pF,avgF,ho,me) ! rain deformation from perturbed state call grainDeformation(pF,avgF,ho,en) ! rain deformation from perturbed state
call stressPenalty(pR,DevNull, avgF,pF,ho,me) ! stress penalty due to interface mismatch from perturbed state call stressPenalty(pR,DevNull, avgF,pF,ho,en) ! stress penalty due to interface mismatch from perturbed state
call volumePenalty(pD,devNull(1,1), avgF,pF,nGrain) ! stress penalty due to volume discrepancy from perturbed state call volumePenalty(pD,devNull(1,1), avgF,pF,nGrain) ! stress penalty due to volume discrepancy from perturbed state
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -424,7 +424,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
iGr3N = faceID(2:4) ! identify the grain ID in local coordinate system (3-dimensional index) iGr3N = faceID(2:4) ! identify the grain ID in local coordinate system (3-dimensional index)
iGrN = grain3to1(iGr3N,param(ho)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index) iGrN = grain3to1(iGr3N,param(ho)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index)
intFaceN = getInterface(2*faceID(1),iGr3N) ! identify the interface ID of the grain intFaceN = getInterface(2*faceID(1),iGr3N) ! identify the interface ID of the grain
normN = interfaceNormal(intFaceN,ho,me) normN = interfaceNormal(intFaceN,ho,en)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! identify the right/up/front grain (+|P) ! identify the right/up/front grain (+|P)
@ -432,7 +432,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identify the grain ID in local coordinate system (3-dimensional index) iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identify the grain ID in local coordinate system (3-dimensional index)
iGrP = grain3to1(iGr3P,param(ho)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index) iGrP = grain3to1(iGr3P,param(ho)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index)
intFaceP = getInterface(2*faceID(1)-1,iGr3P) ! identify the interface ID of the grain intFaceP = getInterface(2*faceID(1)-1,iGr3P) ! identify the interface ID of the grain
normP = interfaceNormal(intFaceP,ho,me) normP = interfaceNormal(intFaceP,ho,en)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! compute the residual stress (contribution of mismatch and volume penalties) from perturbed state ! compute the residual stress (contribution of mismatch and volume penalties) from perturbed state
@ -472,7 +472,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
do i = 1,3*nIntFaceTot;do j = 1,3*nIntFaceTot do i = 1,3*nIntFaceTot;do j = 1,3*nIntFaceTot
drelax(i) = drelax(i) - jnverse(i,j)*resid(j) ! Calculate the correction for the state variable drelax(i) = drelax(i) - jnverse(i,j)*resid(j) ! Calculate the correction for the state variable
enddo; enddo enddo; enddo
stt%relaxationVector(:,me) = relax + drelax ! Updateing the state variable for the next iteration stt%relaxationVector(:,en) = relax + drelax ! Updateing the state variable for the next iteration
if (any(abs(drelax) > num%maxdRelax)) then ! Forcing cutback when the incremental change of relaxation vector becomes too large if (any(abs(drelax) > num%maxdRelax)) then ! Forcing cutback when the incremental change of relaxation vector becomes too large
doneAndHappy = [.true.,.false.] doneAndHappy = [.true.,.false.]
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
@ -488,14 +488,14 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
!------------------------------------------------------------------------------------------------ !------------------------------------------------------------------------------------------------
!> @brief calculate stress-like penalty due to deformation mismatch !> @brief calculate stress-like penalty due to deformation mismatch
!------------------------------------------------------------------------------------------------ !------------------------------------------------------------------------------------------------
subroutine stressPenalty(rPen,nMis,avgF,fDef,ho,me) subroutine stressPenalty(rPen,nMis,avgF,fDef,ho,en)
real(pReal), dimension (:,:,:), intent(out) :: rPen !< stress-like penalty real(pReal), dimension (:,:,:), intent(out) :: rPen !< stress-like penalty
real(pReal), dimension (:,:), intent(out) :: nMis !< total amount of mismatch real(pReal), dimension (:,:), intent(out) :: nMis !< total amount of mismatch
real(pReal), dimension (:,:,:), intent(in) :: fDef !< deformation gradients real(pReal), dimension (:,:,:), intent(in) :: fDef !< deformation gradients
real(pReal), dimension (3,3), intent(in) :: avgF !< initial effective stretch tensor real(pReal), dimension (3,3), intent(in) :: avgF !< initial effective stretch tensor
integer, intent(in) :: ho, me integer, intent(in) :: ho, en
integer, dimension (4) :: intFace integer, dimension (4) :: intFace
integer, dimension (3) :: iGrain3,iGNghb3,nGDim integer, dimension (3) :: iGrain3,iGNghb3,nGDim
@ -515,7 +515,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
! get the correction factor the modulus of penalty stress representing the evolution of area of ! get the correction factor the modulus of penalty stress representing the evolution of area of
! the interfaces due to deformations ! the interfaces due to deformations
surfCorr = surfaceCorrection(avgF,ho,me) surfCorr = surfaceCorrection(avgF,ho,en)
associate(prm => param(ho)) associate(prm => param(ho))
@ -527,7 +527,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
interfaceLoop: do iFace = 1,6 interfaceLoop: do iFace = 1,6
intFace = getInterface(iFace,iGrain3) ! get the 4-dimensional index of the interface in local numbering system of the grain intFace = getInterface(iFace,iGrain3) ! get the 4-dimensional index of the interface in local numbering system of the grain
nVect = interfaceNormal(intFace,ho,me) nVect = interfaceNormal(intFace,ho,en)
iGNghb3 = iGrain3 ! identify the neighboring grain across the interface iGNghb3 = iGrain3 ! identify the neighboring grain across the interface
iGNghb3(abs(intFace(1))) = iGNghb3(abs(intFace(1))) & iGNghb3(abs(intFace(1))) = iGNghb3(abs(intFace(1))) &
+ int(real(intFace(1),pReal)/real(abs(intFace(1)),pReal)) + int(real(intFace(1),pReal)/real(abs(intFace(1)),pReal))
@ -611,14 +611,14 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
!> @brief compute the correction factor accouted for surface evolution (area change) due to !> @brief compute the correction factor accouted for surface evolution (area change) due to
! deformation ! deformation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function surfaceCorrection(avgF,ho,me) function surfaceCorrection(avgF,ho,en)
real(pReal), dimension(3) :: surfaceCorrection real(pReal), dimension(3) :: surfaceCorrection
real(pReal), dimension(3,3), intent(in) :: avgF !< average F real(pReal), dimension(3,3), intent(in) :: avgF !< average F
integer, intent(in) :: & integer, intent(in) :: &
ho, & ho, &
me en
real(pReal), dimension(3,3) :: invC real(pReal), dimension(3,3) :: invC
real(pReal), dimension(3) :: nVect real(pReal), dimension(3) :: nVect
real(pReal) :: detF real(pReal) :: detF
@ -629,7 +629,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
surfaceCorrection = 0.0_pReal surfaceCorrection = 0.0_pReal
do iBase = 1,3 do iBase = 1,3
nVect = interfaceNormal([iBase,1,1,1],ho,me) nVect = interfaceNormal([iBase,1,1,1],ho,en)
do i = 1,3; do j = 1,3 do i = 1,3; do j = 1,3
surfaceCorrection(iBase) = surfaceCorrection(iBase) + invC(i,j)*nVect(i)*nVect(j) ! compute the component of (the inverse of) the stretch in the direction of the normal surfaceCorrection(iBase) = surfaceCorrection(iBase) + invC(i,j)*nVect(i)*nVect(j) ! compute the component of (the inverse of) the stretch in the direction of the normal
enddo; enddo enddo; enddo
@ -651,7 +651,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
real(pReal), dimension(6,6) :: C real(pReal), dimension(6,6) :: C
C = phase_homogenizedC(material_phaseAt2(grainID,ce),material_phaseMemberAt2(grainID,ce)) C = phase_homogenizedC(material_phaseID(grainID,ce),material_phaseEntry(grainID,ce))
equivalentMu = lattice_equivalent_mu(C,'voigt') equivalentMu = lattice_equivalent_mu(C,'voigt')
end function equivalentMu end function equivalentMu
@ -661,14 +661,14 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
!> @brief calculating the grain deformation gradient (the same with !> @brief calculating the grain deformation gradient (the same with
! homogenization_RGC_partitionDeformation, but used only for perturbation scheme) ! homogenization_RGC_partitionDeformation, but used only for perturbation scheme)
!------------------------------------------------------------------------------------------------- !-------------------------------------------------------------------------------------------------
subroutine grainDeformation(F, avgF, ho, me) subroutine grainDeformation(F, avgF, ho, en)
real(pReal), dimension(:,:,:), intent(out) :: F !< partitioned F per grain real(pReal), dimension(:,:,:), intent(out) :: F !< partitioned F per grain
real(pReal), dimension(:,:), intent(in) :: avgF !< averaged F real(pReal), dimension(:,:), intent(in) :: avgF !< averaged F
integer, intent(in) :: & integer, intent(in) :: &
ho, & ho, &
me en
real(pReal), dimension(3) :: aVect,nVect real(pReal), dimension(3) :: aVect,nVect
integer, dimension(4) :: intFace integer, dimension(4) :: intFace
@ -685,8 +685,8 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa
iGrain3 = grain1to3(iGrain,prm%N_constituents) iGrain3 = grain1to3(iGrain,prm%N_constituents)
do iFace = 1,6 do iFace = 1,6
intFace = getInterface(iFace,iGrain3) intFace = getInterface(iFace,iGrain3)
aVect = relaxationVector(intFace,ho,me) aVect = relaxationVector(intFace,ho,en)
nVect = interfaceNormal(intFace,ho,me) nVect = interfaceNormal(intFace,ho,en)
forall (i=1:3,j=1:3) & forall (i=1:3,j=1:3) &
F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! effective relaxations F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! effective relaxations
enddo enddo
@ -753,11 +753,11 @@ end subroutine mechanical_RGC_results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief collect relaxation vectors of an interface !> @brief collect relaxation vectors of an interface
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function relaxationVector(intFace,ho,me) pure function relaxationVector(intFace,ho,en)
real(pReal), dimension (3) :: relaxationVector real(pReal), dimension (3) :: relaxationVector
integer, intent(in) :: ho,me integer, intent(in) :: ho,en
integer, dimension(4), intent(in) :: intFace !< set of interface ID in 4D array (normal and position) integer, dimension(4), intent(in) :: intFace !< set of interface ID in 4D array (normal and position)
integer :: iNum integer :: iNum
@ -770,7 +770,7 @@ pure function relaxationVector(intFace,ho,me)
iNum = interface4to1(intFace,prm%N_constituents) ! identify the position of the interface in global state array iNum = interface4to1(intFace,prm%N_constituents) ! identify the position of the interface in global state array
if (iNum > 0) then if (iNum > 0) then
relaxationVector = stt%relaxationVector((3*iNum-2):(3*iNum),me) relaxationVector = stt%relaxationVector((3*iNum-2):(3*iNum),en)
else else
relaxationVector = 0.0_pReal relaxationVector = 0.0_pReal
endif endif
@ -783,14 +783,14 @@ end function relaxationVector
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief identify the normal of an interface !> @brief identify the normal of an interface
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function interfaceNormal(intFace,ho,me) pure function interfaceNormal(intFace,ho,en)
real(pReal), dimension(3) :: interfaceNormal real(pReal), dimension(3) :: interfaceNormal
integer, dimension(4), intent(in) :: intFace !< interface ID in 4D array (normal and position) integer, dimension(4), intent(in) :: intFace !< interface ID in 4D array (normal and position)
integer, intent(in) :: & integer, intent(in) :: &
ho, & ho, &
me en
integer :: nPos integer :: nPos
associate (dst => dependentState(ho)) associate (dst => dependentState(ho))
@ -801,7 +801,7 @@ pure function interfaceNormal(intFace,ho,me)
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 = matmul(dst%orientation(1:3,1:3,me),interfaceNormal) ! map the normal vector into sample coordinate system (basis) interfaceNormal = matmul(dst%orientation(1:3,1:3,en),interfaceNormal) ! map the normal vector into sample coordinate system (basis)
end associate end associate

View File

@ -44,8 +44,8 @@ module subroutine thermal_init()
allocate(current(configHomogenizations%length)) allocate(current(configHomogenizations%length))
do ho = 1, configHomogenizations%length do ho = 1, configHomogenizations%length
allocate(current(ho)%T(count(material_homogenizationAt2==ho)), source=300.0_pReal) allocate(current(ho)%T(count(material_homogenizationID==ho)), source=300.0_pReal)
allocate(current(ho)%dot_T(count(material_homogenizationAt2==ho)), source=0.0_pReal) allocate(current(ho)%dot_T(count(material_homogenizationID==ho)), source=0.0_pReal)
configHomogenization => configHomogenizations%get(ho) configHomogenization => configHomogenizations%get(ho)
associate(prm => param(ho)) associate(prm => param(ho))
if (configHomogenization%contains('thermal')) then if (configHomogenization%contains('thermal')) then
@ -75,9 +75,9 @@ module subroutine thermal_partition(ce)
integer :: co integer :: co
T = current(material_homogenizationAt2(ce))%T(material_homogenizationMemberAt2(ce)) T = current(material_homogenizationID(ce))%T(material_homogenizationEntry(ce))
dot_T = current(material_homogenizationAt2(ce))%dot_T(material_homogenizationMemberAt2(ce)) dot_T = current(material_homogenizationID(ce))%dot_T(material_homogenizationEntry(ce))
do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) do co = 1, homogenization_Nconstituents(material_homogenizationID(ce))
call phase_thermal_setField(T,dot_T,co,ce) call phase_thermal_setField(T,dot_T,co,ce)
enddo enddo
@ -109,11 +109,11 @@ module function thermal_conduction_getConductivity(ce) result(K)
K = 0.0_pReal K = 0.0_pReal
do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) do co = 1, homogenization_Nconstituents(material_homogenizationID(ce))
K = K + crystallite_push33ToRef(co,ce,lattice_K(:,:,material_phaseAt2(co,ce))) K = K + crystallite_push33ToRef(co,ce,lattice_K(:,:,material_phaseID(co,ce)))
enddo enddo
K = K / real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal) K = K / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
end function thermal_conduction_getConductivity end function thermal_conduction_getConductivity
@ -131,11 +131,11 @@ module function thermal_conduction_getSpecificHeat(ce) result(c_P)
c_P = 0.0_pReal c_P = 0.0_pReal
do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) do co = 1, homogenization_Nconstituents(material_homogenizationID(ce))
c_P = c_P + lattice_c_p(material_phaseAt2(co,ce)) c_P = c_P + lattice_c_p(material_phaseID(co,ce))
enddo enddo
c_P = c_P / real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal) c_P = c_P / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
end function thermal_conduction_getSpecificHeat end function thermal_conduction_getSpecificHeat
@ -153,11 +153,11 @@ module function thermal_conduction_getMassDensity(ce) result(rho)
rho = 0.0_pReal rho = 0.0_pReal
do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) do co = 1, homogenization_Nconstituents(material_homogenizationID(ce))
rho = rho + lattice_rho(material_phaseAt2(co,ce)) rho = rho + lattice_rho(material_phaseID(co,ce))
enddo enddo
rho = rho / real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal) rho = rho / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal)
end function thermal_conduction_getMassDensity end function thermal_conduction_getMassDensity
@ -172,8 +172,8 @@ module subroutine homogenization_thermal_setField(T,dot_T, ce)
real(pReal), intent(in) :: T, dot_T real(pReal), intent(in) :: T, dot_T
current(material_homogenizationAt2(ce))%T(material_homogenizationMemberAt2(ce)) = T current(material_homogenizationID(ce))%T(material_homogenizationEntry(ce)) = T
current(material_homogenizationAt2(ce))%dot_T(material_homogenizationMemberAt2(ce)) = dot_T current(material_homogenizationID(ce))%dot_T(material_homogenizationEntry(ce)) = dot_T
end subroutine homogenization_thermal_setField end subroutine homogenization_thermal_setField
@ -207,7 +207,7 @@ module function homogenization_thermal_T(ce) result(T)
integer, intent(in) :: ce integer, intent(in) :: ce
real(pReal) :: T real(pReal) :: T
T = current(material_homogenizationAt2(ce))%T(material_homogenizationMemberAt2(ce)) T = current(material_homogenizationID(ce))%T(material_homogenizationEntry(ce))
end function homogenization_thermal_T end function homogenization_thermal_T

View File

@ -28,14 +28,14 @@ module material
integer, dimension(:), allocatable, public, protected :: & ! (elem) integer, dimension(:), allocatable, public, protected :: & ! (elem)
material_homogenizationAt, & !< homogenization ID of each element material_homogenizationAt, & !< homogenization ID of each element
material_homogenizationAt2, & !< per cell material_homogenizationID, & !< per cell
material_homogenizationMemberAt2 !< cell material_homogenizationEntry !< cell
integer, dimension(:,:), allocatable :: & ! (ip,elem) integer, dimension(:,:), allocatable :: & ! (ip,elem)
material_homogenizationMemberAt !< position of the element within its homogenization instance material_homogenizationMemberAt !< position of the element within its homogenization instance
integer, dimension(:,:), allocatable, public, protected :: & ! (constituent,elem) integer, dimension(:,:), allocatable, public, protected :: & ! (constituent,elem)
material_phaseAt, & !< phase ID of each element material_phaseAt, & !< phase ID of each element
material_phaseAt2, & !< per constituent,cell material_phaseID, & !< per constituent,cell
material_phaseMemberAt2 !< per constituent, cell material_phaseEntry !< per constituent, cell
integer, dimension(:,:,:), allocatable, public, protected :: & ! (constituent,IP,elem) integer, dimension(:,:,:), allocatable, public, protected :: & ! (constituent,IP,elem)
material_phaseMemberAt !< position of the element within its phase instance material_phaseMemberAt !< position of the element within its phase instance
@ -117,10 +117,10 @@ subroutine material_parseMaterial
allocate(material_phaseMemberAt(homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems),source=0) allocate(material_phaseMemberAt(homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems),source=0)
allocate(material_homogenizationAt2(discretization_nIPs*discretization_Nelems),source=0) allocate(material_homogenizationID(discretization_nIPs*discretization_Nelems),source=0)
allocate(material_homogenizationMemberAt2(discretization_nIPs*discretization_Nelems),source=0) allocate(material_homogenizationEntry(discretization_nIPs*discretization_Nelems),source=0)
allocate(material_phaseAt2(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0) allocate(material_phaseID(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0)
allocate(material_phaseMemberAt2(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0) allocate(material_phaseEntry(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0)
do el = 1, discretization_Nelems do el = 1, discretization_Nelems
material => materials%get(discretization_materialAt(el)) material => materials%get(discretization_materialAt(el))
@ -131,8 +131,8 @@ subroutine material_parseMaterial
ce = (el-1)*discretization_nIPs + ip ce = (el-1)*discretization_nIPs + ip
counterHomogenization(material_homogenizationAt(el)) = counterHomogenization(material_homogenizationAt(el)) + 1 counterHomogenization(material_homogenizationAt(el)) = counterHomogenization(material_homogenizationAt(el)) + 1
material_homogenizationMemberAt(ip,el) = counterHomogenization(material_homogenizationAt(el)) material_homogenizationMemberAt(ip,el) = counterHomogenization(material_homogenizationAt(el))
material_homogenizationAt2(ce) = material_homogenizationAt(el) material_homogenizationID(ce) = material_homogenizationAt(el)
material_homogenizationMemberAt2(ce) = material_homogenizationMemberAt(ip,el) material_homogenizationEntry(ce) = material_homogenizationMemberAt(ip,el)
enddo enddo
frac = 0.0_pReal frac = 0.0_pReal
@ -146,8 +146,8 @@ subroutine material_parseMaterial
counterPhase(material_phaseAt(co,el)) = counterPhase(material_phaseAt(co,el)) + 1 counterPhase(material_phaseAt(co,el)) = counterPhase(material_phaseAt(co,el)) + 1
material_phaseMemberAt(co,ip,el) = counterPhase(material_phaseAt(co,el)) material_phaseMemberAt(co,ip,el) = counterPhase(material_phaseAt(co,el))
material_phaseAt2(co,ce) = material_phaseAt(co,el) material_phaseID(co,ce) = material_phaseAt(co,el)
material_phaseMemberAt2(co,ce) = material_phaseMemberAt(co,ip,el) material_phaseEntry(co,ce) = material_phaseMemberAt(co,ip,el)
enddo enddo
enddo enddo

View File

@ -419,10 +419,10 @@ subroutine phase_restore(ce,includeL)
co co
do co = 1,homogenization_Nconstituents(material_homogenizationAt2(ce)) do co = 1,homogenization_Nconstituents(material_homogenizationID(ce))
if (damageState(material_phaseAt2(co,ce))%sizeState > 0) & if (damageState(material_phaseID(co,ce))%sizeState > 0) &
damageState(material_phaseAt2(co,ce))%state( :,material_phasememberAt2(co,ce)) = & damageState(material_phaseID(co,ce))%state( :,material_phaseEntry(co,ce)) = &
damageState(material_phaseAt2(co,ce))%state0(:,material_phasememberAt2(co,ce)) damageState(material_phaseID(co,ce))%state0(:,material_phaseEntry(co,ce))
enddo enddo
call mechanical_restore(ce,includeL) call mechanical_restore(ce,includeL)
@ -473,7 +473,6 @@ subroutine crystallite_init()
integer :: & integer :: &
ph, & ph, &
me, &
co, & !< counter in integration point component loop co, & !< counter in integration point component loop
ip, & !< counter in integration point loop ip, & !< counter in integration point loop
el, & !< counter in element loop el, & !< counter in element loop
@ -539,12 +538,10 @@ subroutine crystallite_init()
flush(IO_STDOUT) flush(IO_STDOUT)
!$OMP PARALLEL DO PRIVATE(ph,me) !$OMP PARALLEL DO
do el = 1, size(material_phaseMemberAt,3) do el = 1, size(material_phaseMemberAt,3)
do ip = 1, size(material_phaseMemberAt,2) do ip = 1, size(material_phaseMemberAt,2)
do co = 1,homogenization_Nconstituents(material_homogenizationAt(el)) do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
ph = material_phaseAt(co,el)
me = material_phaseMemberAt(co,ip,el)
call crystallite_orientations(co,ip,el) call crystallite_orientations(co,ip,el)
call plastic_dependentState(co,ip,el) ! update dependent state variables to be consistent with basic states call plastic_dependentState(co,ip,el) ! update dependent state variables to be consistent with basic states
enddo enddo
@ -590,11 +587,11 @@ function crystallite_push33ToRef(co,ce, tensor33)
real(pReal), dimension(3,3) :: crystallite_push33ToRef real(pReal), dimension(3,3) :: crystallite_push33ToRef
real(pReal), dimension(3,3) :: T real(pReal), dimension(3,3) :: T
integer :: ph, me integer :: ph, en
ph = material_phaseAt2(co,ce) ph = material_phaseID(co,ce)
me = material_phaseMemberAt2(co,ce) en = material_phaseEntry(co,ce)
T = matmul(material_orientation0(co,ph,me)%asMatrix(),transpose(math_inv33(phase_mechanical_getF(co,ce)))) ! ToDo: initial orientation correct? T = matmul(material_orientation0(co,ph,en)%asMatrix(),transpose(math_inv33(phase_mechanical_getF(co,ce)))) ! ToDo: initial orientation correct?
crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T)) crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T))

View File

@ -151,7 +151,7 @@ module subroutine damage_init
do ph = 1,phases%length do ph = 1,phases%length
Nmembers = count(material_phaseAt2 == ph) Nmembers = count(material_phaseID == ph)
allocate(current(ph)%phi(Nmembers),source=1.0_pReal) allocate(current(ph)%phi(Nmembers),source=1.0_pReal)
allocate(current(ph)%d_phi_d_dot_phi(Nmembers),source=0.0_pReal) allocate(current(ph)%d_phi_d_dot_phi(Nmembers),source=0.0_pReal)
@ -199,9 +199,9 @@ module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi,
phiDot = 0.0_pReal phiDot = 0.0_pReal
dPhiDot_dPhi = 0.0_pReal dPhiDot_dPhi = 0.0_pReal
do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) do co = 1, homogenization_Nconstituents(material_homogenizationID(ce))
ph = material_phaseAt2(co,ce) ph = material_phaseID(co,ce)
me = material_phasememberAt2(co,ce) me = material_phaseEntry(co,ce)
select case(phase_source(ph)) select case(phase_source(ph))
case (DAMAGE_ISOBRITTLE_ID) case (DAMAGE_ISOBRITTLE_ID)
@ -477,7 +477,7 @@ module subroutine phase_damage_set_phi(phi,co,ce)
integer, intent(in) :: ce, co integer, intent(in) :: ce, co
current(material_phaseAt2(co,ce))%phi(material_phaseMemberAt2(co,ce)) = phi current(material_phaseID(co,ce))%phi(material_phaseEntry(co,ce)) = phi
end subroutine phase_damage_set_phi end subroutine phase_damage_set_phi

View File

@ -35,7 +35,7 @@ module function anisoductile_init() result(mySources)
pl, & pl, &
sources, & sources, &
src src
integer :: Ninstances,Nmembers,p integer :: Ninstances,Nmembers,ph
integer, dimension(:), allocatable :: N_sl integer, dimension(:), allocatable :: N_sl
character(len=pStringLen) :: extmsg = '' character(len=pStringLen) :: extmsg = ''
@ -50,15 +50,15 @@ module function anisoductile_init() result(mySources)
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(param(phases%length)) allocate(param(phases%length))
do p = 1, phases%length do ph = 1, phases%length
if(mySources(p)) then if(mySources(ph)) then
phase => phases%get(p) phase => phases%get(ph)
mech => phase%get('mechanical') mech => phase%get('mechanical')
pl => mech%get('plastic') pl => mech%get('plastic')
sources => phase%get('damage') sources => phase%get('damage')
associate(prm => param(p)) associate(prm => param(ph))
src => sources%get(1) src => sources%get(1)
N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray) N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray)
@ -78,10 +78,10 @@ module function anisoductile_init() result(mySources)
if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q' if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q'
if (any(prm%gamma_crit < 0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit' if (any(prm%gamma_crit < 0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit'
Nmembers=count(material_phaseAt2==p) Nmembers=count(material_phaseID==ph)
call phase_allocateState(damageState(p),Nmembers,1,1,0) call phase_allocateState(damageState(ph),Nmembers,1,1,0)
damageState(p)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal) damageState(ph)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal)
if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol' if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol'
end associate end associate

View File

@ -31,7 +31,7 @@ module function isobrittle_init() result(mySources)
phase, & phase, &
sources, & sources, &
src src
integer :: Nmembers,p integer :: Nmembers,ph
character(len=pStringLen) :: extmsg = '' character(len=pStringLen) :: extmsg = ''
@ -45,12 +45,12 @@ module function isobrittle_init() result(mySources)
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(param(phases%length)) allocate(param(phases%length))
do p = 1, phases%length do ph = 1, phases%length
if(mySources(p)) then if(mySources(ph)) then
phase => phases%get(p) phase => phases%get(ph)
sources => phase%get('damage') sources => phase%get('damage')
associate(prm => param(p)) associate(prm => param(ph))
src => sources%get(1) src => sources%get(1)
prm%W_crit = src%get_asFloat('W_crit') prm%W_crit = src%get_asFloat('W_crit')
@ -64,10 +64,10 @@ module function isobrittle_init() result(mySources)
! sanity checks ! sanity checks
if (prm%W_crit <= 0.0_pReal) extmsg = trim(extmsg)//' W_crit' if (prm%W_crit <= 0.0_pReal) extmsg = trim(extmsg)//' W_crit'
Nmembers = count(material_phaseAt2==p) Nmembers = count(material_phaseID==ph)
call phase_allocateState(damageState(p),Nmembers,1,1,1) call phase_allocateState(damageState(ph),Nmembers,1,1,1)
damageState(p)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal) damageState(ph)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal)
if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol' if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol'
end associate end associate

View File

@ -33,7 +33,7 @@ module function isoductile_init() result(mySources)
phase, & phase, &
sources, & sources, &
src src
integer :: Ninstances,Nmembers,p integer :: Ninstances,Nmembers,ph
character(len=pStringLen) :: extmsg = '' character(len=pStringLen) :: extmsg = ''
@ -47,12 +47,12 @@ module function isoductile_init() result(mySources)
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(param(phases%length)) allocate(param(phases%length))
do p = 1, phases%length do ph = 1, phases%length
if(mySources(p)) then if(mySources(ph)) then
phase => phases%get(p) phase => phases%get(ph)
sources => phase%get('damage') sources => phase%get('damage')
associate(prm => param(p)) associate(prm => param(ph))
src => sources%get(1) src => sources%get(1)
prm%q = src%get_asFloat('q') prm%q = src%get_asFloat('q')
@ -68,10 +68,10 @@ module function isoductile_init() result(mySources)
if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q' if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q'
if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit' if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit'
Nmembers=count(material_phaseAt2==p) Nmembers=count(material_phaseID==ph)
call phase_allocateState(damageState(p),Nmembers,1,1,0) call phase_allocateState(damageState(ph),Nmembers,1,1,0)
damageState(p)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal) damageState(ph)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal)
if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol' if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol'
end associate end associate

View File

@ -1217,9 +1217,9 @@ module subroutine mechanical_restore(ce,includeL)
co, ph, me co, ph, me
do co = 1,homogenization_Nconstituents(material_homogenizationAt2(ce)) do co = 1,homogenization_Nconstituents(material_homogenizationID(ce))
ph = material_phaseAt2(co,ce) ph = material_phaseID(co,ce)
me = material_phaseMemberAt2(co,ce) me = material_phaseEntry(co,ce)
if (includeL) then if (includeL) then
phase_mechanical_Lp(ph)%data(1:3,1:3,me) = phase_mechanical_Lp0(ph)%data(1:3,1:3,me) phase_mechanical_Lp(ph)%data(1:3,1:3,me) = phase_mechanical_Lp0(ph)%data(1:3,1:3,me)
phase_mechanical_Li(ph)%data(1:3,1:3,me) = phase_mechanical_Li0(ph)%data(1:3,1:3,me) phase_mechanical_Li(ph)%data(1:3,1:3,me) = phase_mechanical_Li0(ph)%data(1:3,1:3,me)
@ -1267,8 +1267,8 @@ module function phase_mechanical_dPdF(dt,co,ce) result(dPdF)
logical :: error logical :: error
ph = material_phaseAt2(co,ce) ph = material_phaseID(co,ce)
me = material_phaseMemberAt2(co,ce) me = material_phaseEntry(co,ce)
call phase_hooke_SandItsTangents(devNull,dSdFe,dSdFi, & call phase_hooke_SandItsTangents(devNull,dSdFe,dSdFi, &
phase_mechanical_Fe(ph)%data(1:3,1:3,me), & phase_mechanical_Fe(ph)%data(1:3,1:3,me), &
@ -1432,7 +1432,7 @@ module function phase_mechanical_getF(co,ce) result(F)
real(pReal), dimension(3,3) :: F real(pReal), dimension(3,3) :: F
F = phase_mechanical_F(material_phaseAt2(co,ce))%data(1:3,1:3,material_phaseMemberAt2(co,ce)) F = phase_mechanical_F(material_phaseID(co,ce))%data(1:3,1:3,material_phaseEntry(co,ce))
end function phase_mechanical_getF end function phase_mechanical_getF
@ -1461,7 +1461,7 @@ module function phase_mechanical_getP(co,ce) result(P)
real(pReal), dimension(3,3) :: P real(pReal), dimension(3,3) :: P
P = phase_mechanical_P(material_phaseAt2(co,ce))%data(1:3,1:3,material_phaseMemberAt2(co,ce)) P = phase_mechanical_P(material_phaseID(co,ce))%data(1:3,1:3,material_phaseEntry(co,ce))
end function phase_mechanical_getP end function phase_mechanical_getP
@ -1473,7 +1473,7 @@ module subroutine phase_mechanical_setF(F,co,ce)
integer, intent(in) :: co, ce integer, intent(in) :: co, ce
phase_mechanical_F(material_phaseAt2(co,ce))%data(1:3,1:3,material_phaseMemberAt2(co,ce)) = F phase_mechanical_F(material_phaseID(co,ce))%data(1:3,1:3,material_phaseEntry(co,ce)) = F
end subroutine phase_mechanical_setF end subroutine phase_mechanical_setF

View File

@ -220,7 +220,7 @@ module function plastic_dislotungsten_init() result(myPlasticity)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
Nmembers = count(material_phaseAt2 == ph) Nmembers = count(material_phaseID == ph)
sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl
sizeState = sizeDotState sizeState = sizeDotState

View File

@ -406,7 +406,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
Nmembers = count(material_phaseAt2 == ph) Nmembers = count(material_phaseID == ph)
sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl & sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl &
+ size(['f_tw']) * prm%sum_N_tw & + size(['f_tw']) * prm%sum_N_tw &
+ size(['f_tr']) * prm%sum_N_tr + size(['f_tr']) * prm%sum_N_tr

View File

@ -119,7 +119,7 @@ module function plastic_isotropic_init() result(myPlasticity)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
Nmembers = count(material_phaseAt2 == ph) Nmembers = count(material_phaseID == ph)
sizeDotState = size(['xi ','gamma']) sizeDotState = size(['xi ','gamma'])
sizeState = sizeDotState sizeState = sizeDotState

View File

@ -165,7 +165,7 @@ module function plastic_kinehardening_init() result(myPlasticity)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
Nmembers = count(material_phaseAt2 == ph) Nmembers = count(material_phaseID == ph)
sizeDotState = size(['crss ','crss_back', 'accshear ']) * prm%sum_N_sl !ToDo: adjust names like in material.yaml sizeDotState = size(['crss ','crss_back', 'accshear ']) * prm%sum_N_sl !ToDo: adjust names like in material.yaml
sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%sum_N_sl !ToDo: adjust names like in material.yaml sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%sum_N_sl !ToDo: adjust names like in material.yaml
sizeState = sizeDotState + sizeDeltaState sizeState = sizeDotState + sizeDeltaState

View File

@ -30,7 +30,7 @@ module function plastic_none_init() result(myPlasticity)
phases => config_material%get('phase') phases => config_material%get('phase')
do ph = 1, phases%length do ph = 1, phases%length
if(.not. myPlasticity(ph)) cycle if(.not. myPlasticity(ph)) cycle
call phase_allocateState(plasticState(ph),count(material_phaseAt2 == ph),0,0,0) call phase_allocateState(plasticState(ph),count(material_phaseID == ph),0,0,0)
enddo enddo
end function plastic_none_init end function plastic_none_init

View File

@ -398,7 +398,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
Nmembers = count(material_phaseAt2 == ph) Nmembers = count(material_phaseID == ph)
sizeDotState = size([ 'rhoSglEdgePosMobile ','rhoSglEdgeNegMobile ', & sizeDotState = size([ 'rhoSglEdgePosMobile ','rhoSglEdgeNegMobile ', &
'rhoSglScrewPosMobile ','rhoSglScrewNegMobile ', & 'rhoSglScrewPosMobile ','rhoSglScrewNegMobile ', &
'rhoSglEdgePosImmobile ','rhoSglEdgeNegImmobile ', & 'rhoSglEdgePosImmobile ','rhoSglEdgeNegImmobile ', &
@ -527,7 +527,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
if(.not. myPlasticity(ph)) cycle if(.not. myPlasticity(ph)) cycle
phase => phases%get(ph) phase => phases%get(ph)
Nmembers = count(material_phaseAt2 == ph) Nmembers = count(material_phaseID == ph)
l = 0 l = 0
do t = 1,4 do t = 1,4
do s = 1,param(ph)%sum_N_sl do s = 1,param(ph)%sum_N_sl
@ -1824,9 +1824,9 @@ subroutine storeGeometry(ph)
V = reshape(IPvolume,[product(shape(IPvolume))]) V = reshape(IPvolume,[product(shape(IPvolume))])
do ce = 1, size(material_homogenizationMemberAt2,1) do ce = 1, size(material_homogenizationEntry,1)
do co = 1, homogenization_maxNconstituents do co = 1, homogenization_maxNconstituents
if (material_phaseAt2(co,ce) == ph) geom(ph)%V_0(material_phaseMemberAt2(co,ce)) = V(ce) if (material_phaseID(co,ce) == ph) geom(ph)%V_0(material_phaseEntry(co,ce)) = V(ce)
enddo enddo
enddo enddo

View File

@ -223,7 +223,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
Nmembers = count(material_phaseAt2 == ph) Nmembers = count(material_phaseID == ph)
sizeDotState = size(['xi_sl ','gamma_sl']) * prm%sum_N_sl & sizeDotState = size(['xi_sl ','gamma_sl']) * prm%sum_N_sl &
+ size(['xi_tw ','gamma_tw']) * prm%sum_N_tw + size(['xi_tw ','gamma_tw']) * prm%sum_N_tw
sizeState = sizeDotState sizeState = sizeDotState

View File

@ -89,7 +89,7 @@ module subroutine thermal_init(phases)
allocate(thermal_Nsources(phases%length),source = 0) allocate(thermal_Nsources(phases%length),source = 0)
do ph = 1, phases%length do ph = 1, phases%length
Nmembers = count(material_phaseAt2 == ph) Nmembers = count(material_phaseID == ph)
allocate(current(ph)%T(Nmembers),source=300.0_pReal) allocate(current(ph)%T(Nmembers),source=300.0_pReal)
allocate(current(ph)%dot_T(Nmembers),source=0.0_pReal) allocate(current(ph)%dot_T(Nmembers),source=0.0_pReal)
phase => phases%get(ph) phase => phases%get(ph)
@ -268,8 +268,8 @@ module subroutine phase_thermal_setField(T,dot_T, co,ce)
integer, intent(in) :: ce, co integer, intent(in) :: ce, co
current(material_phaseAt2(co,ce))%T(material_phaseMemberAt2(co,ce)) = T current(material_phaseID(co,ce))%T(material_phaseEntry(co,ce)) = T
current(material_phaseAt2(co,ce))%dot_T(material_phaseMemberAt2(co,ce)) = dot_T current(material_phaseID(co,ce))%dot_T(material_phaseEntry(co,ce)) = dot_T
end subroutine phase_thermal_setField end subroutine phase_thermal_setField

View File

@ -54,7 +54,7 @@ module function dissipation_init(source_length) result(mySources)
src => sources%get(so) src => sources%get(so)
prm%kappa = src%get_asFloat('kappa') prm%kappa = src%get_asFloat('kappa')
Nmembers = count(material_phaseAt2 == ph) Nmembers = count(material_phaseID == ph)
call phase_allocateState(thermalState(ph)%p(so),Nmembers,0,0,0) call phase_allocateState(thermalState(ph)%p(so),Nmembers,0,0,0)
end associate end associate

View File

@ -67,7 +67,7 @@ module function externalheat_init(source_length) result(mySources)
prm%f_T = src%get_as1dFloat('f_T',requiredSize = size(prm%t_n)) prm%f_T = src%get_as1dFloat('f_T',requiredSize = size(prm%t_n))
Nmembers = count(material_phaseAt2 == ph) Nmembers = count(material_phaseID == ph)
call phase_allocateState(thermalState(ph)%p(so),Nmembers,1,1,0) call phase_allocateState(thermalState(ph)%p(so),Nmembers,1,1,0)
end associate end associate
endif endif