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

renaming of "instance" and thelike to be consistent with naming scheme introduced in rev2662

This commit is contained in:
Christoph Kords 2013-11-21 13:35:43 +00:00
parent 883669bd77
commit 3cc6781eaf
1 changed files with 200 additions and 200 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

400
code/constitutive_nonlocal.f90
View File

@ -1148,13 +1148,13 @@ do matID = 1_pInt,maxNinstance
do j = 1_pInt,2_pInt do j = 1_pInt,2_pInt
noise(j) = math_sampleGaussVar(0.0_pReal, rhoSglScatter(matID)) noise(j) = math_sampleGaussVar(0.0_pReal, rhoSglScatter(matID))
enddo enddo
state(1,i,e)%p(iRhoU(s,1,matID)) = rhoSglEdgePos0(f, matID) + noise(1) state(1,i,e)%p(iRhoU(s,1,matID)) = rhoSglEdgePos0(f,matID) + noise(1)
state(1,i,e)%p(iRhoU(s,2,matID)) = rhoSglEdgeNeg0(f, matID) + noise(1) state(1,i,e)%p(iRhoU(s,2,matID)) = rhoSglEdgeNeg0(f,matID) + noise(1)
state(1,i,e)%p(iRhoU(s,3,matID)) = rhoSglScrewPos0(f, matID) + noise(2) state(1,i,e)%p(iRhoU(s,3,matID)) = rhoSglScrewPos0(f,matID) + noise(2)
state(1,i,e)%p(iRhoU(s,4,matID)) = rhoSglScrewNeg0(f, matID) + noise(2) state(1,i,e)%p(iRhoU(s,4,matID)) = rhoSglScrewNeg0(f,matID) + noise(2)
enddo enddo
state(1,i,e)%p(iRhoD(from:upto,1,matID)) = rhoDipEdge0(f, matID) state(1,i,e)%p(iRhoD(from:upto,1,matID)) = rhoDipEdge0(f,matID)
state(1,i,e)%p(iRhoD(from:upto,2,matID)) = rhoDipScrew0(f, matID) state(1,i,e)%p(iRhoD(from:upto,2,matID)) = rhoDipScrew0(f,matID)
enddo enddo
endif endif
enddo enddo
@ -1271,7 +1271,7 @@ use lattice, only: &
implicit none implicit none
!*** input variables !*** input variables
integer(pInt), intent(in) :: gr, & ! current grain ID integer(pInt), intent(in) :: gr, & ! current grain ID
ip, & ! current integration point ip, & ! current integration point
el ! current element el ! current element
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
@ -1285,23 +1285,23 @@ type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), in
!*** output variables !*** output variables
!*** local variables !*** local variables
integer(pInt) neighboring_el, & ! element number of neighboring material point integer(pInt) neighbor_el, & ! element number of neighboring material point
neighboring_ip, & ! integration point of neighboring material point neighbor_ip, & ! integration point of neighboring material point
instance, & ! my instance of this plasticity matID, & ! my instance of this plasticity
neighboring_instance, & ! instance of this plasticity of neighboring material point neighbor_matID, & ! instance of this plasticity of neighboring material point
latticeStruct, & ! my lattice structure structID, & ! my lattice structure
neighboring_latticeStruct, & ! lattice structure of neighboring material point neighbor_structID, & ! lattice structure of neighboring material point
phase, & phase, &
neighboring_phase, & neighbor_phaseID, &
ns, & ! total number of active slip systems at my material point ns, & ! total number of active slip systems at my material point
neighboring_ns, & ! total number of active slip systems at neighboring material point neighbor_ns, & ! total number of active slip systems at neighboring material point
c, & ! index of dilsocation character (edge, screw) c, & ! index of dilsocation character (edge, screw)
s, & ! slip system index s, & ! slip system index
t, & ! index of dilsocation type (e+, e-, s+, s-, used e+, used e-, used s+, used s-) t, & ! index of dilsocation type (e+, e-, s+, s-, used e+, used e-, used s+, used s-)
dir, & dir, &
n, & n, &
nRealNeighbors ! number of really existing neighbors nRealNeighbors ! number of really existing neighbors
integer(pInt), dimension(2) :: neighbor integer(pInt), dimension(2) :: neighbors
real(pReal) detFe, & real(pReal) detFe, &
detFp, & detFp, &
FVsize, & FVsize, &
@ -1333,32 +1333,32 @@ real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(gr,ip,
totalNslip(phase_plasticityInstance(material_phase(gr,ip,el)))) :: & totalNslip(phase_plasticityInstance(material_phase(gr,ip,el)))) :: &
myInteractionMatrix ! corrected slip interaction matrix myInteractionMatrix ! corrected slip interaction matrix
real(pReal), dimension(2,maxval(totalNslip),mesh_maxNipNeighbors) :: & real(pReal), dimension(2,maxval(totalNslip),mesh_maxNipNeighbors) :: &
neighboring_rhoExcess, & ! excess density at neighboring material point neighbor_rhoExcess, & ! excess density at neighboring material point
neighboring_rhoTotal ! total density at neighboring material point neighbor_rhoTotal ! total density at neighboring material point
real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(gr,ip,el))),2) :: & real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(gr,ip,el))),2) :: &
m ! direction of dislocation motion m ! direction of dislocation motion
logical inversionError logical inversionError
phase = material_phase(gr,ip,el) phase = material_phase(gr,ip,el)
instance = phase_plasticityInstance(phase) matID = phase_plasticityInstance(phase)
latticeStruct = constitutive_nonlocal_structure(instance) structID = constitutive_nonlocal_structure(matID)
ns = totalNslip(instance) ns = totalNslip(matID)
!*** get basic states !*** get basic states
forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) forall (s = 1_pInt:ns, t = 1_pInt:4_pInt)
rhoSgl(s,t) = max(state(gr,ip,el)%p(iRhoU(s,t,instance)), 0.0_pReal) ! ensure positive single mobile densities rhoSgl(s,t) = max(state(gr,ip,el)%p(iRhoU(s,t,matID)), 0.0_pReal) ! ensure positive single mobile densities
rhoSgl(s,t+4_pInt) = state(gr,ip,el)%p(iRhoB(s,t,instance)) rhoSgl(s,t+4_pInt) = state(gr,ip,el)%p(iRhoB(s,t,matID))
endforall endforall
forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) & forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) &
rhoDip(s,c) = max(state(gr,ip,el)%p(iRhoD(s,c,instance)), 0.0_pReal) ! ensure positive dipole densities rhoDip(s,c) = max(state(gr,ip,el)%p(iRhoD(s,c,matID)), 0.0_pReal) ! ensure positive dipole densities
where (abs(rhoSgl) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(instance) & where (abs(rhoSgl) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(matID) &
.or. abs(rhoSgl) < significantRho(instance)) & .or. abs(rhoSgl) < significantRho(matID)) &
rhoSgl = 0.0_pReal rhoSgl = 0.0_pReal
where (abs(rhoDip) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(instance) & where (abs(rhoDip) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(matID) &
.or. abs(rhoDip) < significantRho(instance)) & .or. abs(rhoDip) < significantRho(matID)) &
rhoDip = 0.0_pReal rhoDip = 0.0_pReal
@ -1367,9 +1367,9 @@ where (abs(rhoDip) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(instance
forall (s = 1_pInt:ns) & forall (s = 1_pInt:ns) &
rhoForest(s) = dot_product((sum(abs(rhoSgl(1:ns,[1,2,5,6])),2) + rhoDip(1:ns,1)), & rhoForest(s) = dot_product((sum(abs(rhoSgl(1:ns,[1,2,5,6])),2) + rhoDip(1:ns,1)), &
forestProjectionEdge(s,1:ns,instance)) & forestProjectionEdge(s,1:ns,matID)) &
+ dot_product((sum(abs(rhoSgl(1:ns,[3,4,7,8])),2) + rhoDip(1:ns,2)), & + dot_product((sum(abs(rhoSgl(1:ns,[3,4,7,8])),2) + rhoDip(1:ns,2)), &
forestProjectionScrew(s,1:ns,instance)) forestProjectionScrew(s,1:ns,matID))
@ -1378,19 +1378,19 @@ forall (s = 1_pInt:ns) &
!*** (see Kubin,Devincre,Hoc; 2008; Modeling dislocation storage rates and mean free paths in face-centered cubic crystals) !*** (see Kubin,Devincre,Hoc; 2008; Modeling dislocation storage rates and mean free paths in face-centered cubic crystals)
myInteractionMatrix = 0.0_pReal myInteractionMatrix = 0.0_pReal
myInteractionMatrix(1:ns,1:ns) = interactionMatrixSlipSlip(1:ns,1:ns,instance) myInteractionMatrix(1:ns,1:ns) = interactionMatrixSlipSlip(1:ns,1:ns,matID)
if (latticeStruct < 3_pInt) then ! only fcc and bcc if (structID < 3_pInt) then ! only fcc and bcc
do s = 1_pInt,ns do s = 1_pInt,ns
myRhoForest = max(rhoForest(s),significantRho(instance)) myRhoForest = max(rhoForest(s),significantRho(matID))
correction = ( 1.0_pReal - linetensionEffect(instance) & correction = ( 1.0_pReal - linetensionEffect(matID) &
+ linetensionEffect(instance) & + linetensionEffect(matID) &
* log(0.35_pReal * burgers(s,instance) * sqrt(myRhoForest)) & * log(0.35_pReal * burgers(s,matID) * sqrt(myRhoForest)) &
/ log(0.35_pReal * burgers(s,instance) * 1e6_pReal)) ** 2.0_pReal / log(0.35_pReal * burgers(s,matID) * 1e6_pReal)) ** 2.0_pReal
myInteractionMatrix(s,1:ns) = correction * myInteractionMatrix(s,1:ns) myInteractionMatrix(s,1:ns) = correction * myInteractionMatrix(s,1:ns)
enddo enddo
endif endif
forall (s = 1_pInt:ns) & forall (s = 1_pInt:ns) &
tauThreshold(s) = mu(instance) * burgers(s,instance) & tauThreshold(s) = mu(matID) * burgers(s,matID) &
* sqrt(dot_product((sum(abs(rhoSgl),2) + sum(abs(rhoDip),2)), myInteractionMatrix(s,1:ns))) * sqrt(dot_product((sum(abs(rhoSgl),2) + sum(abs(rhoDip),2)), myInteractionMatrix(s,1:ns)))
@ -1400,7 +1400,7 @@ forall (s = 1_pInt:ns) &
tauBack = 0.0_pReal tauBack = 0.0_pReal
if (.not. phase_localPlasticity(phase) .and. shortRangeStressCorrection(instance)) then if (.not. phase_localPlasticity(phase) .and. shortRangeStressCorrection(matID)) then
call math_invert33(Fe, invFe, detFe, inversionError) call math_invert33(Fe, invFe, detFe, inversionError)
call math_invert33(Fp, invFp, detFp, inversionError) call math_invert33(Fp, invFp, detFp, inversionError)
rhoExcess(1,1:ns) = rhoSgl(1:ns,1) - rhoSgl(1:ns,2) rhoExcess(1,1:ns) = rhoSgl(1:ns,1) - rhoSgl(1:ns,2)
@ -1410,36 +1410,36 @@ if (.not. phase_localPlasticity(phase) .and. shortRangeStressCorrection(instance
!* loop through my neighborhood and get the connection vectors (in lattice frame) and the excess densities !* loop through my neighborhood and get the connection vectors (in lattice frame) and the excess densities
nRealNeighbors = 0_pInt nRealNeighbors = 0_pInt
neighboring_rhoTotal = 0.0_pReal neighbor_rhoTotal = 0.0_pReal
do n = 1_pInt,FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,el)))) do n = 1_pInt,FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,el))))
neighboring_el = mesh_ipNeighborhood(1,n,ip,el) neighbor_el = mesh_ipNeighborhood(1,n,ip,el)
neighboring_ip = mesh_ipNeighborhood(2,n,ip,el) neighbor_ip = mesh_ipNeighborhood(2,n,ip,el)
if (neighboring_el > 0 .and. neighboring_ip > 0) then if (neighbor_el > 0 .and. neighbor_ip > 0) then
neighboring_phase = material_phase(gr,neighboring_ip,neighboring_el) neighbor_phaseID = material_phase(gr,neighbor_ip,neighbor_el)
neighboring_instance = phase_plasticityInstance(neighboring_phase) neighbor_matID = phase_plasticityInstance(neighbor_phaseID)
neighboring_latticeStruct = constitutive_nonlocal_structure(neighboring_instance) neighbor_structID = constitutive_nonlocal_structure(neighbor_matID)
neighboring_ns = totalNslip(neighboring_instance) neighbor_ns = totalNslip(neighbor_matID)
if (.not. phase_localPlasticity(neighboring_phase) & if (.not. phase_localPlasticity(neighbor_phaseID) &
.and. neighboring_latticeStruct == latticeStruct & .and. neighbor_structID == structID &
.and. neighboring_instance == instance) then .and. neighbor_matID == matID) then
if (neighboring_ns == ns) then if (neighbor_ns == ns) then
nRealNeighbors = nRealNeighbors + 1_pInt nRealNeighbors = nRealNeighbors + 1_pInt
forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) forall (s = 1_pInt:ns, c = 1_pInt:2_pInt)
neighboring_rhoExcess(c,s,n) = & neighbor_rhoExcess(c,s,n) = &
max(state(gr,neighboring_ip,neighboring_el)%p(iRhoU(s,2*c-1,neighboring_instance)), 0.0_pReal) &! positive mobiles max(state(gr,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c-1,neighbor_matID)), 0.0_pReal) &! positive mobiles
- max(state(gr,neighboring_ip,neighboring_el)%p(iRhoU(s,2*c,neighboring_instance)), 0.0_pReal) ! negative mobiles - max(state(gr,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c,neighbor_matID)), 0.0_pReal) ! negative mobiles
neighboring_rhoTotal(c,s,n) = & neighbor_rhoTotal(c,s,n) = &
max(state(gr,neighboring_ip,neighboring_el)%p(iRhoU(s,2*c-1,neighboring_instance)), 0.0_pReal) &! positive mobiles max(state(gr,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c-1,neighbor_matID)), 0.0_pReal) &! positive mobiles
+ max(state(gr,neighboring_ip,neighboring_el)%p(iRhoU(s,2*c,neighboring_instance)), 0.0_pReal) & ! negative mobiles + max(state(gr,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c,neighbor_matID)), 0.0_pReal) & ! negative mobiles
+ abs(state(gr,neighboring_ip,neighboring_el)%p(iRhoB(s,2*c-1,neighboring_instance))) & ! positive deads + abs(state(gr,neighbor_ip,neighbor_el)%p(iRhoB(s,2*c-1,neighbor_matID))) & ! positive deads
+ abs(state(gr,neighboring_ip,neighboring_el)%p(iRhoB(s,2*c,neighboring_instance))) & ! negative deads + abs(state(gr,neighbor_ip,neighbor_el)%p(iRhoB(s,2*c,neighbor_matID))) & ! negative deads
+ max(state(gr,neighboring_ip,neighboring_el)%p(iRhoD(s,c,neighboring_instance)), 0.0_pReal) ! dipoles + max(state(gr,neighbor_ip,neighbor_el)%p(iRhoD(s,c,neighbor_matID)), 0.0_pReal) ! dipoles
endforall endforall
connection_latticeConf(1:3,n) = & connection_latticeConf(1:3,n) = &
math_mul33x3(invFe, mesh_ipCoordinates(1:3,neighboring_ip,neighboring_el) & math_mul33x3(invFe, mesh_ipCoordinates(1:3,neighbor_ip,neighbor_el) &
- mesh_ipCoordinates(1:3,ip,el)) - mesh_ipCoordinates(1:3,ip,el))
normal_latticeConf = math_mul33x3(math_transpose33(invFp), mesh_ipAreaNormal(1:3,n,ip,el)) normal_latticeConf = math_mul33x3(math_transpose33(invFp), mesh_ipAreaNormal(1:3,n,ip,el))
if (math_mul3x3(normal_latticeConf,connection_latticeConf(1:3,n)) < 0.0_pReal) then ! neighbor connection points in opposite direction to face normal: must be periodic image if (math_mul3x3(normal_latticeConf,connection_latticeConf(1:3,n)) < 0.0_pReal) then ! neighboring connection points in opposite direction to face normal: must be periodic image
connection_latticeConf(1:3,n) = normal_latticeConf * mesh_ipVolume(ip,el) & connection_latticeConf(1:3,n) = normal_latticeConf * mesh_ipVolume(ip,el) &
/ mesh_ipArea(n,ip,el) ! instead take the surface normal scaled with the diameter of the cell / mesh_ipArea(n,ip,el) ! instead take the surface normal scaled with the diameter of the cell
endif endif
@ -1450,12 +1450,12 @@ if (.not. phase_localPlasticity(phase) .and. shortRangeStressCorrection(instance
else else
! local neighbor or different lattice structure or different constitution instance -> use central values instead ! local neighbor or different lattice structure or different constitution instance -> use central values instead
connection_latticeConf(1:3,n) = 0.0_pReal connection_latticeConf(1:3,n) = 0.0_pReal
neighboring_rhoExcess(1:2,1:ns,n) = rhoExcess neighbor_rhoExcess(1:2,1:ns,n) = rhoExcess
endif endif
else else
! free surface -> use central values instead ! free surface -> use central values instead
connection_latticeConf(1:3,n) = 0.0_pReal connection_latticeConf(1:3,n) = 0.0_pReal
neighboring_rhoExcess(1:2,1:ns,n) = rhoExcess neighbor_rhoExcess(1:2,1:ns,n) = rhoExcess
endif endif
enddo enddo
@ -1464,8 +1464,8 @@ if (.not. phase_localPlasticity(phase) .and. shortRangeStressCorrection(instance
!* 1. interpolation of the excess density in the neighorhood !* 1. interpolation of the excess density in the neighorhood
!* 2. interpolation of the dead dislocation density in the central volume !* 2. interpolation of the dead dislocation density in the central volume
m(1:3,1:ns,1) = lattice_sd(1:3,slipSystemLattice(1:ns,instance),latticeStruct) m(1:3,1:ns,1) = lattice_sd(1:3,slipSystemLattice(1:ns,matID),structID)
m(1:3,1:ns,2) = -lattice_st(1:3,slipSystemLattice(1:ns,instance),latticeStruct) m(1:3,1:ns,2) = -lattice_st(1:3,slipSystemLattice(1:ns,matID),structID)
do s = 1_pInt,ns do s = 1_pInt,ns
@ -1473,12 +1473,12 @@ if (.not. phase_localPlasticity(phase) .and. shortRangeStressCorrection(instance
do c = 1_pInt,2_pInt do c = 1_pInt,2_pInt
do dir = 1_pInt,3_pInt do dir = 1_pInt,3_pInt
neighbor(1) = 2_pInt * dir - 1_pInt neighbors(1) = 2_pInt * dir - 1_pInt
neighbor(2) = 2_pInt * dir neighbors(2) = 2_pInt * dir
connections(dir,1:3) = connection_latticeConf(1:3,neighbor(1)) & connections(dir,1:3) = connection_latticeConf(1:3,neighbors(1)) &
- connection_latticeConf(1:3,neighbor(2)) - connection_latticeConf(1:3,neighbors(2))
rhoExcessDifferences(dir) = neighboring_rhoExcess(c,s,neighbor(1)) & rhoExcessDifferences(dir) = neighbor_rhoExcess(c,s,neighbors(1)) &
- neighboring_rhoExcess(c,s,neighbor(2)) - neighbor_rhoExcess(c,s,neighbors(2))
enddo enddo
call math_invert33(connections,invConnections,temp,inversionError) call math_invert33(connections,invConnections,temp,inversionError)
if (inversionError) then if (inversionError) then
@ -1499,15 +1499,15 @@ if (.not. phase_localPlasticity(phase) .and. shortRangeStressCorrection(instance
rhoExcessGradient_over_rho = 0.0_pReal rhoExcessGradient_over_rho = 0.0_pReal
forall (c = 1_pInt:2_pInt) & forall (c = 1_pInt:2_pInt) &
rhoTotal(c) = (sum(abs(rhoSgl(s,[2*c-1,2*c,2*c+3,2*c+4]))) + rhoDip(s,c) + sum(neighboring_rhoTotal(c,s,:))) & rhoTotal(c) = (sum(abs(rhoSgl(s,[2*c-1,2*c,2*c+3,2*c+4]))) + rhoDip(s,c) + sum(neighbor_rhoTotal(c,s,:))) &
/ real(1_pInt + nRealNeighbors,pReal) / real(1_pInt + nRealNeighbors,pReal)
forall (c = 1_pInt:2_pInt, rhoTotal(c) > 0.0_pReal) & forall (c = 1_pInt:2_pInt, rhoTotal(c) > 0.0_pReal) &
rhoExcessGradient_over_rho(c) = rhoExcessGradient(c) / rhoTotal(c) rhoExcessGradient_over_rho(c) = rhoExcessGradient(c) / rhoTotal(c)
!* gives the local stress correction when multiplied with a factor !* gives the local stress correction when multiplied with a factor
tauBack(s) = - mu(instance) * burgers(s,instance) / (2.0_pReal * pi) & tauBack(s) = - mu(matID) * burgers(s,matID) / (2.0_pReal * pi) &
* (rhoExcessGradient_over_rho(1) / (1.0_pReal - nu(instance)) + rhoExcessGradient_over_rho(2)) * (rhoExcessGradient_over_rho(1) / (1.0_pReal - nu(matID)) + rhoExcessGradient_over_rho(2))
enddo enddo
endif endif
@ -1515,9 +1515,9 @@ endif
!*** set dependent states !*** set dependent states
state(gr,ip,el)%p(iRhoF(1:ns,instance)) = rhoForest state(gr,ip,el)%p(iRhoF(1:ns,matID)) = rhoForest
state(gr,ip,el)%p(iTauF(1:ns,instance)) = tauThreshold state(gr,ip,el)%p(iTauF(1:ns,matID)) = tauThreshold
state(gr,ip,el)%p(iTauB(1:ns,instance)) = tauBack state(gr,ip,el)%p(iTauB(1:ns,matID)) = tauBack
#ifndef _OPENMP #ifndef _OPENMP
@ -1557,7 +1557,7 @@ use material, only: material_phase, &
implicit none implicit none
!*** input variables !*** input variables
integer(pInt), intent(in) :: ipc, & !< current grain number integer(pInt), intent(in) :: ipc, & !< current grain number
ip, & !< current integration point ip, & !< current integration point
el, & !< current element number el, & !< current element number
c !< dislocation character (1:edge, 2:screw) c !< dislocation character (1:edge, 2:screw)
@ -1574,7 +1574,7 @@ real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip
dv_dtauNS !< velocity derivative with respect to resolved shear stress (including non Schmid contributions) dv_dtauNS !< velocity derivative with respect to resolved shear stress (including non Schmid contributions)
!*** local variables !*** local variables
integer(pInt) :: instance, & !< current instance of this plasticity integer(pInt) :: matID, & !< current instance of this plasticity
ns, & !< short notation for the total number of active slip systems ns, & !< short notation for the total number of active slip systems
s !< index of my current slip system s !< index of my current slip system
real(pReal) tauRel_P, & real(pReal) tauRel_P, &
@ -1600,8 +1600,8 @@ real(pReal) tauRel_P, &
mobility !< dislocation mobility mobility !< dislocation mobility
instance = phase_plasticityInstance(material_phase(ipc,ip,el)) matID = phase_plasticityInstance(material_phase(ipc,ip,el))
ns = totalNslip(instance) ns = totalNslip(matID)
v = 0.0_pReal v = 0.0_pReal
dv_dtau = 0.0_pReal dv_dtau = 0.0_pReal
@ -1617,20 +1617,20 @@ if (Temperature > 0.0_pReal) then
!* The derivative only gives absolute values; the correct sign is taken care of in the formula for the derivative of the velocity !* The derivative only gives absolute values; the correct sign is taken care of in the formula for the derivative of the velocity
tauEff = max(0.0_pReal, abs(tauNS(s)) - tauThreshold(s)) ! ensure that the effective stress is positive tauEff = max(0.0_pReal, abs(tauNS(s)) - tauThreshold(s)) ! ensure that the effective stress is positive
meanfreepath_P = burgers(s,instance) meanfreepath_P = burgers(s,matID)
jumpWidth_P = burgers(s,instance) jumpWidth_P = burgers(s,matID)
activationLength_P = doublekinkwidth(instance) * burgers(s,instance) activationLength_P = doublekinkwidth(matID) * burgers(s,matID)
activationVolume_P = activationLength_P * jumpWidth_P * burgers(s,instance) activationVolume_P = activationLength_P * jumpWidth_P * burgers(s,matID)
criticalStress_P = peierlsStress(s,c,instance) criticalStress_P = peierlsStress(s,c,matID)
activationEnergy_P = criticalStress_P * activationVolume_P activationEnergy_P = criticalStress_P * activationVolume_P
tauRel_P = min(1.0_pReal, tauEff / criticalStress_P) ! ensure that the activation probability cannot become greater than one tauRel_P = min(1.0_pReal, tauEff / criticalStress_P) ! ensure that the activation probability cannot become greater than one
tPeierls = 1.0_pReal / fattack(instance) & tPeierls = 1.0_pReal / fattack(matID) &
* exp(activationEnergy_P / (KB * Temperature) & * exp(activationEnergy_P / (KB * Temperature) &
* (1.0_pReal - tauRel_P**pParam(instance))**qParam(instance)) * (1.0_pReal - tauRel_P**pParam(matID))**qParam(matID))
if (tauEff < criticalStress_P) then if (tauEff < criticalStress_P) then
dtPeierls_dtau = tPeierls * pParam(instance) * qParam(instance) * activationVolume_P / (KB * Temperature) & dtPeierls_dtau = tPeierls * pParam(matID) * qParam(matID) * activationVolume_P / (KB * Temperature) &
* (1.0_pReal - tauRel_P**pParam(instance))**(qParam(instance)-1.0_pReal) & * (1.0_pReal - tauRel_P**pParam(matID))**(qParam(matID)-1.0_pReal) &
* tauRel_P**(pParam(instance)-1.0_pReal) * tauRel_P**(pParam(matID)-1.0_pReal)
else else
dtPeierls_dtau = 0.0_pReal dtPeierls_dtau = 0.0_pReal
endif endif
@ -1640,21 +1640,21 @@ if (Temperature > 0.0_pReal) then
!* The derivative only gives absolute values; the correct sign is taken care of in the formula for the derivative of the velocity !* The derivative only gives absolute values; the correct sign is taken care of in the formula for the derivative of the velocity
tauEff = abs(tau(s)) - tauThreshold(s) tauEff = abs(tau(s)) - tauThreshold(s)
meanfreepath_S = burgers(s,instance) / sqrt(solidSolutionConcentration(instance)) meanfreepath_S = burgers(s,matID) / sqrt(solidSolutionConcentration(matID))
jumpWidth_S = solidSolutionSize(instance) * burgers(s,instance) jumpWidth_S = solidSolutionSize(matID) * burgers(s,matID)
activationLength_S = burgers(s,instance) / sqrt(solidSolutionConcentration(instance)) activationLength_S = burgers(s,matID) / sqrt(solidSolutionConcentration(matID))
activationVolume_S = activationLength_S * jumpWidth_S * burgers(s,instance) activationVolume_S = activationLength_S * jumpWidth_S * burgers(s,matID)
activationEnergy_S = solidSolutionEnergy(instance) activationEnergy_S = solidSolutionEnergy(matID)
criticalStress_S = activationEnergy_S / activationVolume_S criticalStress_S = activationEnergy_S / activationVolume_S
tauRel_S = min(1.0_pReal, tauEff / criticalStress_S) ! ensure that the activation probability cannot become greater than one tauRel_S = min(1.0_pReal, tauEff / criticalStress_S) ! ensure that the activation probability cannot become greater than one
tSolidSolution = 1.0_pReal / fattack(instance) & tSolidSolution = 1.0_pReal / fattack(matID) &
* exp(activationEnergy_S / (KB * Temperature) & * exp(activationEnergy_S / (KB * Temperature) &
* (1.0_pReal - tauRel_S**pParam(instance))**qParam(instance)) * (1.0_pReal - tauRel_S**pParam(matID))**qParam(matID))
if (tauEff < criticalStress_S) then if (tauEff < criticalStress_S) then
dtSolidSolution_dtau = tSolidSolution * pParam(instance) * qParam(instance) & dtSolidSolution_dtau = tSolidSolution * pParam(matID) * qParam(matID) &
* activationVolume_S / (KB * Temperature) & * activationVolume_S / (KB * Temperature) &
* (1.0_pReal - tauRel_S**pParam(instance))**(qParam(instance)-1.0_pReal) & * (1.0_pReal - tauRel_S**pParam(matID))**(qParam(matID)-1.0_pReal) &
* tauRel_S**(pParam(instance)-1.0_pReal) * tauRel_S**(pParam(matID)-1.0_pReal)
else else
dtSolidSolution_dtau = 0.0_pReal dtSolidSolution_dtau = 0.0_pReal
endif endif
@ -1663,7 +1663,7 @@ if (Temperature > 0.0_pReal) then
!* viscous glide velocity !* viscous glide velocity
tauEff = abs(tau(s)) - tauThreshold(s) tauEff = abs(tau(s)) - tauThreshold(s)
mobility = burgers(s,instance) / viscosity(instance) mobility = burgers(s,matID) / viscosity(matID)
vViscous = mobility * tauEff vViscous = mobility * tauEff
@ -1727,7 +1727,7 @@ use mesh, only: mesh_ipVolume
implicit none implicit none
!*** input variables !*** input variables
integer(pInt), intent(in) :: ipc, & !< current grain number integer(pInt), intent(in) :: ipc, & !< current grain number
ip, & !< current integration point ip, & !< current integration point
el !< current element number el !< current element number
real(pReal), intent(in) :: Temperature !< temperature real(pReal), intent(in) :: Temperature !< temperature
@ -1741,8 +1741,8 @@ real(pReal), dimension(3,3), intent(out) :: Lp !< plast
real(pReal), dimension(9,9), intent(out) :: dLp_dTstar99 !< derivative of Lp with respect to Tstar (9x9 matrix) real(pReal), dimension(9,9), intent(out) :: dLp_dTstar99 !< derivative of Lp with respect to Tstar (9x9 matrix)
!*** local variables !*** local variables
integer(pInt) matID, & !< current instance of this plasticity integer(pInt) matID, & !< current instance of this plasticity
structID, & !< current lattice structure structID, & !< current lattice structure
ns, & !< short notation for the total number of active slip systems ns, & !< short notation for the total number of active slip systems
i, & i, &
j, & j, &
@ -1925,7 +1925,7 @@ use material, only: homogenization_maxNgrains, &
implicit none implicit none
!*** input variables !*** input variables
integer(pInt), intent(in) :: ipc, & ! current grain number integer(pInt), intent(in) :: ipc, & ! current grain number
ip, & ! current integration point ip, & ! current integration point
el ! current element number el ! current element number
real(pReal), dimension(6), intent(in) :: Tstar_v ! current 2nd Piola-Kirchhoff stress in Mandel notation real(pReal), dimension(6), intent(in) :: Tstar_v ! current 2nd Piola-Kirchhoff stress in Mandel notation
@ -1938,8 +1938,8 @@ type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), in
type(p_vec), intent(out) :: deltaState ! change of state variables / microstructure type(p_vec), intent(out) :: deltaState ! change of state variables / microstructure
!*** local variables !*** local variables
integer(pInt) matID, & ! current instance of this plasticity integer(pInt) matID, & ! current instance of this plasticity
structID, & ! current lattice structure structID, & ! current lattice structure
ns, & ! short notation for the total number of active slip systems ns, & ! short notation for the total number of active slip systems
c, & ! character of dislocation c, & ! character of dislocation
t, & ! type of dislocation t, & ! type of dislocation
@ -2137,7 +2137,7 @@ use lattice, only: lattice_Sslip_v, &
implicit none implicit none
!*** input variables !*** input variables
integer(pInt), intent(in) :: ipc, & !< current grain number integer(pInt), intent(in) :: ipc, & !< current grain number
ip, & !< current integration point ip, & !< current integration point
el !< current element number el !< current element number
real(pReal), intent(in) :: Temperature, & !< temperature real(pReal), intent(in) :: Temperature, & !< temperature
@ -2159,9 +2159,9 @@ real(pReal), dimension(constitutive_nonlocal_sizeDotState(phase_plasticityInstan
constitutive_nonlocal_dotState !< evolution of state variables / microstructure constitutive_nonlocal_dotState !< evolution of state variables / microstructure
!*** local variables !*** local variables
integer(pInt) matID, & !< current instance of this plasticity integer(pInt) matID, & !< current instance of this plasticity
neighbor_instance, & !< instance of my neighbor's plasticity neighbor_matID, & !< instance of my neighbor's plasticity
structID, & !< current lattice structure structID, & !< current lattice structure
ns, & !< short notation for the total number of active slip systems ns, & !< short notation for the total number of active slip systems
c, & !< character of dislocation c, & !< character of dislocation
n, & !< index of my current neighbor n, & !< index of my current neighbor
@ -2251,12 +2251,12 @@ gdot = 0.0_pReal
forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) forall (s = 1_pInt:ns, t = 1_pInt:4_pInt)
rhoSgl(s,t) = max(state(ipc,ip,el)%p(iRhoU(s,t,matID)), 0.0_pReal) ! ensure positive single mobile densities rhoSgl(s,t) = max(state(ipc,ip,el)%p(iRhoU(s,t,matID)), 0.0_pReal) ! ensure positive single mobile densities
rhoSgl(s,t+4_pInt) = state(ipc,ip,el)%p(iRhoB(s,t,matID)) rhoSgl(s,t+4_pInt) = state(ipc,ip,el)%p(iRhoB(s,t,matID))
v(s,t) = state(ipc,ip,el)%p(iV(s,t,matID)) v(s,t) = state(ipc,ip,el)%p(iV(s,t,matID))
endforall endforall
forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) forall (s = 1_pInt:ns, c = 1_pInt:2_pInt)
rhoDip(s,c) = max(state(ipc,ip,el)%p(iRhoD(s,c,matID)), 0.0_pReal) ! ensure positive dipole densities rhoDip(s,c) = max(state(ipc,ip,el)%p(iRhoD(s,c,matID)), 0.0_pReal) ! ensure positive dipole densities
endforall endforall
rhoForest = state(ipc,ip,el)%p(iRhoF(1:ns,matID)) rhoForest = state(ipc,ip,el)%p(iRhoF(1:ns,matID))
tauThreshold = state(ipc,ip,el)%p(iTauF(1:ns,matID)) tauThreshold = state(ipc,ip,el)%p(iTauF(1:ns,matID))
@ -2336,14 +2336,14 @@ dUpper = max(dUpper,dLower)
!*** calculate dislocation multiplication !*** calculate dislocation multiplication
rhoDotMultiplication = 0.0_pReal rhoDotMultiplication = 0.0_pReal
if (structID == 2_pInt) then ! BCC if (structID == 2_pInt) then ! BCC
forall (s = 1:ns, sum(abs(v(s,1:4))) > 0.0_pReal) forall (s = 1:ns, sum(abs(v(s,1:4))) > 0.0_pReal)
rhoDotMultiplication(s,1:2) = sum(abs(gdot(s,3:4))) / burgers(s,matID) & ! assuming double-cross-slip of screws to be decisive for multiplication rhoDotMultiplication(s,1:2) = sum(abs(gdot(s,3:4))) / burgers(s,matID) & ! assuming double-cross-slip of screws to be decisive for multiplication
* sqrt(rhoForest(s)) / lambda0(s,matID) ! & ! mean free path * sqrt(rhoForest(s)) / lambda0(s,matID) ! & ! mean free path
! * 2.0_pReal * sum(abs(v(s,3:4))) / sum(abs(v(s,1:4))) ! ratio of screw to overall velocity determines edge generation ! * 2.0_pReal * sum(abs(v(s,3:4))) / sum(abs(v(s,1:4))) ! ratio of screw to overall velocity determines edge generation
rhoDotMultiplication(s,3:4) = sum(abs(gdot(s,3:4))) / burgers(s,matID) & ! assuming double-cross-slip of screws to be decisive for multiplication rhoDotMultiplication(s,3:4) = sum(abs(gdot(s,3:4))) / burgers(s,matID) & ! assuming double-cross-slip of screws to be decisive for multiplication
* sqrt(rhoForest(s)) / lambda0(s,matID) ! & ! mean free path * sqrt(rhoForest(s)) / lambda0(s,matID) ! & ! mean free path
! * 2.0_pReal * sum(abs(v(s,1:2))) / sum(abs(v(s,1:4))) ! ratio of edge to overall velocity determines screw generation ! * 2.0_pReal * sum(abs(v(s,1:2))) / sum(abs(v(s,1:4))) ! ratio of edge to overall velocity determines screw generation
endforall endforall
else ! ALL OTHER STRUCTURES else ! ALL OTHER STRUCTURES
@ -2394,7 +2394,7 @@ endif
rhoDotFlux = 0.0_pReal rhoDotFlux = 0.0_pReal
if (.not. phase_localPlasticity(material_phase(ipc,ip,el))) then ! only for nonlocal plasticity if (.not. phase_localPlasticity(material_phase(ipc,ip,el))) then ! only for nonlocal plasticity
!*** check CFL (Courant-Friedrichs-Lewy) condition for flux !*** check CFL (Courant-Friedrichs-Lewy) condition for flux
@ -2440,7 +2440,7 @@ if (.not. phase_localPlasticity(material_phase(ipc,ip,el))) then
opposite_n = mesh_ipNeighborhood(3,opposite_neighbor,ip,el) opposite_n = mesh_ipNeighborhood(3,opposite_neighbor,ip,el)
if (neighbor_n > 0_pInt) then ! if neighbor exists, average deformation gradient if (neighbor_n > 0_pInt) then ! if neighbor exists, average deformation gradient
neighbor_instance = phase_plasticityInstance(material_phase(ipc,neighbor_ip,neighbor_el)) neighbor_matID = phase_plasticityInstance(material_phase(ipc,neighbor_ip,neighbor_el))
neighbor_Fe = Fe(1:3,1:3,ipc,neighbor_ip,neighbor_el) neighbor_Fe = Fe(1:3,1:3,ipc,neighbor_ip,neighbor_el)
neighbor_F = math_mul33x33(neighbor_Fe, Fp(1:3,1:3,ipc,neighbor_ip,neighbor_el)) neighbor_F = math_mul33x33(neighbor_Fe, Fp(1:3,1:3,ipc,neighbor_ip,neighbor_el))
Favg = 0.5_pReal * (my_F + neighbor_F) Favg = 0.5_pReal * (my_F + neighbor_F)
@ -2466,15 +2466,15 @@ if (.not. phase_localPlasticity(material_phase(ipc,ip,el))) then
if (considerEnteringFlux) then if (considerEnteringFlux) then
if(numerics_timeSyncing .and. (subfrac(ipc,neighbor_ip,neighbor_el) /= subfrac(ipc,ip,el))) then ! for timesyncing: in case of a timestep at the interface we have to use "state0" to make sure that fluxes n both sides are equal if(numerics_timeSyncing .and. (subfrac(ipc,neighbor_ip,neighbor_el) /= subfrac(ipc,ip,el))) then ! for timesyncing: in case of a timestep at the interface we have to use "state0" to make sure that fluxes n both sides are equal
forall (s = 1:ns, t = 1_pInt:4_pInt) forall (s = 1:ns, t = 1_pInt:4_pInt)
neighbor_v(s,t) = state0(ipc,neighbor_ip,neighbor_el)%p(iV(s,t,neighbor_instance)) neighbor_v(s,t) = state0(ipc,neighbor_ip,neighbor_el)%p(iV(s,t,neighbor_matID))
neighbor_rhoSgl(s,t) = max(state0(ipc,neighbor_ip,neighbor_el)%p(iRhoU(s,t,neighbor_instance)), 0.0_pReal) neighbor_rhoSgl(s,t) = max(state0(ipc,neighbor_ip,neighbor_el)%p(iRhoU(s,t,neighbor_matID)), 0.0_pReal)
neighbor_rhoSgl(s,t+4_pInt) = state0(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,t,neighbor_instance)) neighbor_rhoSgl(s,t+4_pInt) = state0(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,t,neighbor_matID))
endforall endforall
else else
forall (s = 1:ns, t = 1_pInt:4_pInt) forall (s = 1:ns, t = 1_pInt:4_pInt)
neighbor_v(s,t) = state(ipc,neighbor_ip,neighbor_el)%p(iV(s,t,neighbor_instance)) neighbor_v(s,t) = state(ipc,neighbor_ip,neighbor_el)%p(iV(s,t,neighbor_matID))
neighbor_rhoSgl(s,t) = max(state(ipc,neighbor_ip,neighbor_el)%p(iRhoU(s,t,neighbor_instance)), 0.0_pReal) neighbor_rhoSgl(s,t) = max(state(ipc,neighbor_ip,neighbor_el)%p(iRhoU(s,t,neighbor_matID)), 0.0_pReal)
neighbor_rhoSgl(s,t+4_pInt) = state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,t,neighbor_instance)) neighbor_rhoSgl(s,t+4_pInt) = state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,t,neighbor_matID))
endforall endforall
endif endif
where (abs(neighbor_rhoSgl) * mesh_ipVolume(neighbor_ip,neighbor_el) ** 0.667_pReal & where (abs(neighbor_rhoSgl) * mesh_ipVolume(neighbor_ip,neighbor_el) ** 0.667_pReal &
@ -2748,12 +2748,12 @@ integer(pInt) Nneighbors, & !
n, & ! neighbor index n, & ! neighbor index
neighbor_e, & ! element index of my neighbor neighbor_e, & ! element index of my neighbor
neighbor_i, & ! integration point index of my neighbor neighbor_i, & ! integration point index of my neighbor
my_phase, & phaseID, &
neighbor_phase, & neighbor_phaseID, &
my_texture, & textureID, &
neighbor_texture, & neighbor_textureID, &
my_structure, & ! lattice structure structID, & ! lattice structure
my_instance, & ! instance of plasticity matID, & ! instance of plasticity
ns, & ! number of active slip systems ns, & ! number of active slip systems
s1, & ! slip system index (me) s1, & ! slip system index (me)
s2 ! slip system index (my neighbor) s2 ! slip system index (my neighbor)
@ -2773,13 +2773,13 @@ logical, dimension(totalNslip(phase_plasticityInstance(material_phase(1,i,e))))
Nneighbors = FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,e)))) Nneighbors = FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,e))))
my_phase = material_phase(1,i,e) phaseID = material_phase(1,i,e)
my_texture = material_texture(1,i,e) textureID = material_texture(1,i,e)
my_instance = phase_plasticityInstance(my_phase) matID = phase_plasticityInstance(phaseID)
my_structure = constitutive_nonlocal_structure(my_instance) structID = constitutive_nonlocal_structure(matID)
ns = totalNslip(my_instance) ns = totalNslip(matID)
slipNormal(1:3,1:ns) = lattice_sn(1:3, slipSystemLattice(1:ns,my_instance), my_structure) slipNormal(1:3,1:ns) = lattice_sn(1:3, slipSystemLattice(1:ns,matID), structID)
slipDirection(1:3,1:ns) = lattice_sd(1:3, slipSystemLattice(1:ns,my_instance), my_structure) slipDirection(1:3,1:ns) = lattice_sd(1:3, slipSystemLattice(1:ns,matID), structID)
!*** start out fully compatible !*** start out fully compatible
@ -2801,7 +2801,7 @@ do n = 1_pInt,Nneighbors
if (neighbor_e <= 0_pInt .or. neighbor_i <= 0_pInt) then if (neighbor_e <= 0_pInt .or. neighbor_i <= 0_pInt) then
forall(s1 = 1_pInt:ns) & forall(s1 = 1_pInt:ns) &
my_compatibility(1:2,s1,s1,n) = sqrt(surfaceTransmissivity(my_instance)) my_compatibility(1:2,s1,s1,n) = sqrt(surfaceTransmissivity(matID))
cycle cycle
endif endif
@ -2812,9 +2812,9 @@ do n = 1_pInt,Nneighbors
!* If one of the two "CPFEM" phases has a local plasticity law, !* If one of the two "CPFEM" phases has a local plasticity law,
!* we do not consider this to be a phase boundary, so completely compatible. !* we do not consider this to be a phase boundary, so completely compatible.
neighbor_phase = material_phase(1,neighbor_i,neighbor_e) neighbor_phaseID = material_phase(1,neighbor_i,neighbor_e)
if (neighbor_phase /= my_phase) then if (neighbor_phaseID /= phaseID) then
if (.not. phase_localPlasticity(neighbor_phase) .and. .not. phase_localPlasticity(my_phase)) then if (.not. phase_localPlasticity(neighbor_phaseID) .and. .not. phase_localPlasticity(phaseID)) then
forall(s1 = 1_pInt:ns) & forall(s1 = 1_pInt:ns) &
my_compatibility(1:2,s1,s1,n) = 0.0_pReal ! = sqrt(0.0) my_compatibility(1:2,s1,s1,n) = 0.0_pReal ! = sqrt(0.0)
endif endif
@ -2825,12 +2825,12 @@ do n = 1_pInt,Nneighbors
!* GRAIN BOUNDARY ! !* GRAIN BOUNDARY !
!* fixed transmissivity for adjacent ips with different texture (only if explicitly given in material.config) !* fixed transmissivity for adjacent ips with different texture (only if explicitly given in material.config)
if (grainboundaryTransmissivity(my_instance) >= 0.0_pReal) then if (grainboundaryTransmissivity(matID) >= 0.0_pReal) then
neighbor_texture = material_texture(1,neighbor_i,neighbor_e) neighbor_textureID = material_texture(1,neighbor_i,neighbor_e)
if (neighbor_texture /= my_texture) then if (neighbor_textureID /= textureID) then
if (.not. phase_localPlasticity(neighbor_phase)) then if (.not. phase_localPlasticity(neighbor_phaseID)) then
forall(s1 = 1_pInt:ns) & forall(s1 = 1_pInt:ns) &
my_compatibility(1:2,s1,s1,n) = sqrt(grainboundaryTransmissivity(my_instance)) my_compatibility(1:2,s1,s1,n) = sqrt(grainboundaryTransmissivity(matID))
endif endif
cycle cycle
endif endif
@ -2910,7 +2910,7 @@ implicit none
!*** input variables !*** input variables
integer(pInt), intent(in) :: ipc, & ! current grain ID integer(pInt), intent(in) :: ipc, & ! current grain ID
ip, & ! current integration point ip, & ! current integration point
el ! current element el ! current element
real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
@ -2926,12 +2926,12 @@ real(pReal), dimension(3,3) :: constitutive_nonlocal_dislocationstress
!*** local variables !*** local variables
integer(pInt) neighbor_el, & ! element number of neighbor material point integer(pInt) neighbor_el, & ! element number of neighbor material point
neighbor_ip, & ! integration point of neighbor material point neighbor_ip, & ! integration point of neighbor material point
instance, & ! my instance of this plasticity matID, & ! my instance of this plasticity
neighbor_instance, & ! instance of this plasticity of neighbor material point neighbor_matID, & ! instance of this plasticity of neighbor material point
latticeStruct, & ! my lattice structure structID, & ! my lattice structure
neighbor_latticeStruct, & ! lattice structure of neighbor material point neighbor_structID, & ! lattice structure of neighbor material point
phase, & phase, &
neighbor_phase, & neighbor_phaseID, &
ns, & ! total number of active slip systems at my material point ns, & ! total number of active slip systems at my material point
neighbor_ns, & ! total number of active slip systems at neighbor material point neighbor_ns, & ! total number of active slip systems at neighbor material point
c, & ! index of dilsocation character (edge, screw) c, & ! index of dilsocation character (edge, screw)
@ -2973,17 +2973,17 @@ real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip
logical inversionError logical inversionError
phase = material_phase(ipc,ip,el) phase = material_phase(ipc,ip,el)
instance = phase_plasticityInstance(phase) matID = phase_plasticityInstance(phase)
latticeStruct = constitutive_nonlocal_structure(instance) structID = constitutive_nonlocal_structure(matID)
ns = totalNslip(instance) ns = totalNslip(matID)
!*** get basic states !*** get basic states
forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) forall (s = 1_pInt:ns, t = 1_pInt:4_pInt)
rhoSgl(s,t) = max(state(ipc,ip,el)%p(iRhoU(s,t,instance)), 0.0_pReal) ! ensure positive single mobile densities rhoSgl(s,t) = max(state(ipc,ip,el)%p(iRhoU(s,t,matID)), 0.0_pReal) ! ensure positive single mobile densities
rhoSgl(s,t+4_pInt) = state(ipc,ip,el)%p(iRhoB(s,t,instance)) rhoSgl(s,t+4_pInt) = state(ipc,ip,el)%p(iRhoB(s,t,matID))
endforall endforall
@ -3007,8 +3007,8 @@ if (.not. phase_localPlasticity(phase)) then
periodicImages = 0_pInt periodicImages = 0_pInt
do dir = 1_pInt,3_pInt do dir = 1_pInt,3_pInt
if (mesh_periodicSurface(dir)) then if (mesh_periodicSurface(dir)) then
periodicImages(1,dir) = floor((coords(dir) - cutoffRadius(instance) - minCoord(dir)) / meshSize(dir), pInt) periodicImages(1,dir) = floor((coords(dir) - cutoffRadius(matID) - minCoord(dir)) / meshSize(dir), pInt)
periodicImages(2,dir) = ceiling((coords(dir) + cutoffRadius(instance) - maxCoord(dir)) / meshSize(dir), pInt) periodicImages(2,dir) = ceiling((coords(dir) + cutoffRadius(matID) - maxCoord(dir)) / meshSize(dir), pInt)
endif endif
enddo enddo
@ -3018,20 +3018,20 @@ if (.not. phase_localPlasticity(phase)) then
do neighbor_el = 1_pInt,mesh_NcpElems do neighbor_el = 1_pInt,mesh_NcpElems
ipLoop: do neighbor_ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,neighbor_el))) ipLoop: do neighbor_ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,neighbor_el)))
neighbor_phase = material_phase(ipc,neighbor_ip,neighbor_el) neighbor_phaseID = material_phase(ipc,neighbor_ip,neighbor_el)
if (phase_localPlasticity(neighbor_phase)) then if (phase_localPlasticity(neighbor_phaseID)) then
cycle cycle
endif endif
neighbor_instance = phase_plasticityInstance(neighbor_phase) neighbor_matID = phase_plasticityInstance(neighbor_phaseID)
neighbor_latticeStruct = constitutive_nonlocal_structure(neighbor_instance) neighbor_structID = constitutive_nonlocal_structure(neighbor_matID)
neighbor_ns = totalNslip(neighbor_instance) neighbor_ns = totalNslip(neighbor_matID)
call math_invert33(Fe(1:3,1:3,1,neighbor_ip,neighbor_el), neighbor_invFe, detFe, inversionError) call math_invert33(Fe(1:3,1:3,1,neighbor_ip,neighbor_el), neighbor_invFe, detFe, inversionError)
neighbor_ipVolumeSideLength = mesh_ipVolume(neighbor_ip,neighbor_el) ** (1.0_pReal/3.0_pReal) ! reference volume used here neighbor_ipVolumeSideLength = mesh_ipVolume(neighbor_ip,neighbor_el) ** (1.0_pReal/3.0_pReal) ! reference volume used here
forall (s = 1_pInt:neighbor_ns, c = 1_pInt:2_pInt) forall (s = 1_pInt:neighbor_ns, c = 1_pInt:2_pInt)
neighbor_rhoExcess(c,1,s) = state(ipc,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c-1,neighbor_instance)) & ! positive mobiles neighbor_rhoExcess(c,1,s) = state(ipc,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c-1,neighbor_matID)) & ! positive mobiles
- state(ipc,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c,neighbor_instance)) ! negative mobiles - state(ipc,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c,neighbor_matID)) ! negative mobiles
neighbor_rhoExcess(c,2,s) = abs(state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,2*c-1,neighbor_instance))) & ! positive deads neighbor_rhoExcess(c,2,s) = abs(state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,2*c-1,neighbor_matID))) & ! positive deads
- abs(state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,2*c,neighbor_instance))) ! negative deads - abs(state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,2*c,neighbor_matID))) ! negative deads
endforall endforall
Tdislo_neighborLattice = 0.0_pReal Tdislo_neighborLattice = 0.0_pReal
do deltaX = periodicImages(1,1),periodicImages(2,1) do deltaX = periodicImages(1,1),periodicImages(2,1)
@ -3048,7 +3048,7 @@ ipLoop: do neighbor_ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,neighbor_el))
+ (/real(deltaX,pReal), real(deltaY,pReal), real(deltaZ,pReal)/) * meshSize + (/real(deltaX,pReal), real(deltaY,pReal), real(deltaZ,pReal)/) * meshSize
connection = neighbor_coords - coords connection = neighbor_coords - coords
distance = sqrt(sum(connection * connection)) distance = sqrt(sum(connection * connection))
if (distance > cutoffRadius(instance)) then if (distance > cutoffRadius(matID)) then
cycle cycle
endif endif
@ -3064,15 +3064,15 @@ ipLoop: do neighbor_ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,neighbor_el))
!* and add up the stress contributions from egde and screw excess on these slip systems (if significant) !* and add up the stress contributions from egde and screw excess on these slip systems (if significant)
do s = 1_pInt,neighbor_ns do s = 1_pInt,neighbor_ns
if (all(abs(neighbor_rhoExcess(:,:,s)) < significantRho(instance))) then if (all(abs(neighbor_rhoExcess(:,:,s)) < significantRho(matID))) then
cycle ! not significant cycle ! not significant
endif endif
!* map the connection vector from the lattice into the slip system frame !* map the connection vector from the lattice into the slip system frame
connection_neighborSlip = math_mul33x3(lattice2slip(1:3,1:3,s,neighbor_instance), & connection_neighborSlip = math_mul33x3(lattice2slip(1:3,1:3,s,neighbor_matID), &
connection_neighborLattice) connection_neighborLattice)
!* edge contribution to stress !* edge contribution to stress
@ -3086,12 +3086,12 @@ ipLoop: do neighbor_ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,neighbor_el))
zsquare = z * z zsquare = z * z
do j = 1_pInt,2_pInt do j = 1_pInt,2_pInt
if (abs(neighbor_rhoExcess(1,j,s)) < significantRho(instance)) then if (abs(neighbor_rhoExcess(1,j,s)) < significantRho(matID)) then
cycle cycle
elseif (j > 1_pInt) then elseif (j > 1_pInt) then
x = connection_neighborSlip(1) + sign(0.5_pReal * segmentLength, & x = connection_neighborSlip(1) + sign(0.5_pReal * segmentLength, &
state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,1,neighbor_instance)) & state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,1,neighbor_matID)) &
- state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,2,neighbor_instance))) - state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,2,neighbor_matID)))
xsquare = x * x xsquare = x * x
endif endif
@ -3111,7 +3111,7 @@ ipLoop: do neighbor_ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,neighbor_el))
* (1.0_pReal + xsquare / Rsquare + xsquare / denominator) & * (1.0_pReal + xsquare / Rsquare + xsquare / denominator) &
* neighbor_rhoExcess(1,j,s) * neighbor_rhoExcess(1,j,s)
sigma(2,2) = sigma(2,2) - real(side,pReal) & sigma(2,2) = sigma(2,2) - real(side,pReal) &
* (flipSign * 2.0_pReal * nu(instance) * z / denominator + z * lambda / Rcube) & * (flipSign * 2.0_pReal * nu(matID) * z / denominator + z * lambda / Rcube) &
* neighbor_rhoExcess(1,j,s) * neighbor_rhoExcess(1,j,s)
sigma(3,3) = sigma(3,3) + real(side,pReal) & sigma(3,3) = sigma(3,3) + real(side,pReal) &
* flipSign * z / denominator & * flipSign * z / denominator &
@ -3124,7 +3124,7 @@ ipLoop: do neighbor_ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,neighbor_el))
* (1.0_pReal - zsquare / Rsquare - zsquare / denominator) & * (1.0_pReal - zsquare / Rsquare - zsquare / denominator) &
* neighbor_rhoExcess(1,j,s) * neighbor_rhoExcess(1,j,s)
sigma(2,3) = sigma(2,3) - real(side,pReal) & sigma(2,3) = sigma(2,3) - real(side,pReal) &
* (nu(instance) / R - zsquare / Rcube) * neighbor_rhoExcess(1,j,s) * (nu(matID) / R - zsquare / Rcube) * neighbor_rhoExcess(1,j,s)
enddo enddo
enddo enddo
@ -3132,12 +3132,12 @@ ipLoop: do neighbor_ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,neighbor_el))
x = connection_neighborSlip(1) ! have to restore this value, because position might have been adapted for edge deads before x = connection_neighborSlip(1) ! have to restore this value, because position might have been adapted for edge deads before
do j = 1_pInt,2_pInt do j = 1_pInt,2_pInt
if (abs(neighbor_rhoExcess(2,j,s)) < significantRho(instance)) then if (abs(neighbor_rhoExcess(2,j,s)) < significantRho(matID)) then
cycle cycle
elseif (j > 1_pInt) then elseif (j > 1_pInt) then
y = connection_neighborSlip(2) + sign(0.5_pReal * segmentLength, & y = connection_neighborSlip(2) + sign(0.5_pReal * segmentLength, &
state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,3,neighbor_instance)) & state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,3,neighbor_matID)) &
- state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,4,neighbor_instance))) - state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,4,neighbor_matID)))
ysquare = y * y ysquare = y * y
endif endif
@ -3152,9 +3152,9 @@ ipLoop: do neighbor_ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,neighbor_el))
exit ipLoop exit ipLoop
endif endif
sigma(1,2) = sigma(1,2) - real(side,pReal) * flipSign * z * (1.0_pReal - nu(instance)) / denominator & sigma(1,2) = sigma(1,2) - real(side,pReal) * flipSign * z * (1.0_pReal - nu(matID)) / denominator &
* neighbor_rhoExcess(2,j,s) * neighbor_rhoExcess(2,j,s)
sigma(1,3) = sigma(1,3) + real(side,pReal) * flipSign * y * (1.0_pReal - nu(instance)) / denominator & sigma(1,3) = sigma(1,3) + real(side,pReal) * flipSign * y * (1.0_pReal - nu(matID)) / denominator &
* neighbor_rhoExcess(2,j,s) * neighbor_rhoExcess(2,j,s)
enddo enddo
enddo enddo
@ -3172,12 +3172,12 @@ ipLoop: do neighbor_ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,neighbor_el))
!* scale stresses and map them into the neighbor material point's lattice configuration !* scale stresses and map them into the neighbor material point's lattice configuration
sigma = sigma * mu(neighbor_instance) * burgers(s,neighbor_instance) & sigma = sigma * mu(neighbor_matID) * burgers(s,neighbor_matID) &
/ (4.0_pReal * pi * (1.0_pReal - nu(neighbor_instance))) & / (4.0_pReal * pi * (1.0_pReal - nu(neighbor_matID))) &
* mesh_ipVolume(neighbor_ip,neighbor_el) / segmentLength ! reference volume is used here (according to the segment length calculation) * mesh_ipVolume(neighbor_ip,neighbor_el) / segmentLength ! reference volume is used here (according to the segment length calculation)
Tdislo_neighborLattice = Tdislo_neighborLattice & Tdislo_neighborLattice = Tdislo_neighborLattice &
+ math_mul33x33(math_transpose33(lattice2slip(1:3,1:3,s,neighbor_instance)), & + math_mul33x33(math_transpose33(lattice2slip(1:3,1:3,s,neighbor_matID)), &
math_mul33x33(sigma, lattice2slip(1:3,1:3,s,neighbor_instance))) math_mul33x33(sigma, lattice2slip(1:3,1:3,s,neighbor_matID)))
enddo ! slip system loop enddo ! slip system loop
@ -3190,22 +3190,22 @@ ipLoop: do neighbor_ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,neighbor_el))
else else
forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) & forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) &
rhoExcessDead(c,s) = state(ipc,ip,el)%p(iRhoB(s,2*c-1,instance)) & ! positive deads (here we use symmetry: if this has negative sign it is treated as negative density at positive position instead of positive density at negative position) rhoExcessDead(c,s) = state(ipc,ip,el)%p(iRhoB(s,2*c-1,matID)) & ! positive deads (here we use symmetry: if this has negative sign it is treated as negative density at positive position instead of positive density at negative position)
+ state(ipc,ip,el)%p(iRhoB(s,2*c,instance)) ! negative deads (here we use symmetry: if this has negative sign it is treated as positive density at positive position instead of negative density at negative position) + state(ipc,ip,el)%p(iRhoB(s,2*c,matID)) ! negative deads (here we use symmetry: if this has negative sign it is treated as positive density at positive position instead of negative density at negative position)
do s = 1_pInt,ns do s = 1_pInt,ns
if (all(abs(rhoExcessDead(:,s)) < significantRho(instance))) then if (all(abs(rhoExcessDead(:,s)) < significantRho(matID))) then
cycle ! not significant cycle ! not significant
endif endif
sigma = 0.0_pReal ! all components except for sigma13 are zero sigma = 0.0_pReal ! all components except for sigma13 are zero
sigma(1,3) = - (rhoExcessDead(1,s) + rhoExcessDead(2,s) * (1.0_pReal - nu(instance))) & sigma(1,3) = - (rhoExcessDead(1,s) + rhoExcessDead(2,s) * (1.0_pReal - nu(matID))) &
* neighbor_ipVolumeSideLength * mu(instance) * burgers(s,instance) & * neighbor_ipVolumeSideLength * mu(matID) * burgers(s,matID) &
/ (sqrt(2.0_pReal) * pi * (1.0_pReal - nu(instance))) / (sqrt(2.0_pReal) * pi * (1.0_pReal - nu(matID)))
sigma(3,1) = sigma(1,3) sigma(3,1) = sigma(1,3)
Tdislo_neighborLattice = Tdislo_neighborLattice & Tdislo_neighborLattice = Tdislo_neighborLattice &
+ math_mul33x33(math_transpose33(lattice2slip(1:3,1:3,s,instance)), & + math_mul33x33(math_transpose33(lattice2slip(1:3,1:3,s,matID)), &
math_mul33x33(sigma, lattice2slip(1:3,1:3,s,instance))) math_mul33x33(sigma, lattice2slip(1:3,1:3,s,matID)))
enddo ! slip system loop enddo ! slip system loop