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

Merge branch 'pheno+' into cmake

This commit is contained in:
zhangc43 2016-05-06 18:26:22 -04:00
parent dbaea4838e c6ddfae4e0
commit eaacde48ba
8 changed files with 689 additions and 579 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
CMakeLists.txt
View File

@ -116,7 +116,7 @@ set (DAMASK_VERSION_MINOR ${DAMASK_V})
# Built-in options for DAMASK build system
# -> can be overwritten from commandline/install_script
option(OPENMP "Use OpenMP libaries for DAMASK" ON )
option(OPTIMIZATION "DAMASK optimization level [OFF,DEFENSIVE,AGGRESSIVE]" "DEFENSIVE" )
option(OPTIMIZATION "DAMASK optimization level [OFF,DEFENSIVE,AGGRESSIVE]" "AGGRESSIVE" )
option(SPECTRAL "Build spectral sovler for DAMASAK" OFF )
option(FEM "Build FEM solver for DAMASK" OFF )

2
Makefile
View File

@ -10,7 +10,7 @@ spectral: build/spectral
build/spectral: build
@mkdir build/spectral
@(cd build/spectral; cmake -Wno-dev -DCMAKE_BUILD_TYPE=RELEASE -DDAMASK_DRIVER=SPECTRAL ../..;)
@(cd build/spectral; cmake -Wno-dev -DCMAKE_BUILD_TYPE=RELEASE -DDAMASK_DRIVER=SPECTRAL -DOPTIMIZATION=AGGRESSIVE -DDAMASK_INSTALL=${HOME}/bin ../..;)
build: bin
@mkdir build

1
src/CMakeLists.txt
View File

@ -79,7 +79,6 @@ add_library(DAMASK_DRIVERS ALIAS DAMASK_LATTICE)
add_library (DAMASK_PLASTIC "plastic_dislotwin.f90"
"plastic_disloUCLA.f90"
"plastic_isotropic.f90"
"plastic_j2.f90"
"plastic_phenopowerlaw.f90"
"plastic_titanmod.f90"
"plastic_nonlocal.f90"

14
src/constitutive.f90
View File

@ -423,7 +423,7 @@ end function constitutive_homogenizedC
!--------------------------------------------------------------------------------------------------
!> @brief calls microstructure function of the different constitutive models
!--------------------------------------------------------------------------------------------------
subroutine constitutive_microstructure(orientations, Fe, Fp, ipc, ip, el)
subroutine constitutive_microstructure(orientations, Fe, Fp, ipc, ip, el, F0s,Fes,Fps,Tstar_vs)
use prec, only: &
pReal
use material, only: &
@ -460,7 +460,15 @@ subroutine constitutive_microstructure(orientations, Fe, Fp, ipc, ip, el)
ho, & !< homogenization
tme !< thermal member position
real(pReal), intent(in), dimension(:,:,:,:) :: &
orientations !< crystal orientations as quaternions
orientations
real(pReal), intent(in), dimension(:,:,:,:,:) :: &
F0s, &
Fes, &
Fps
real(pReal), intent(in), dimension(:,:,:,:) :: &
Tstar_vs !< crystal orientations as quaternions
ho = material_homog(ip,el)
tme = thermalMapping(ho)%p(ip,el)
@ -475,7 +483,7 @@ subroutine constitutive_microstructure(orientations, Fe, Fp, ipc, ip, el)
case (PLASTICITY_NONLOCAL_ID) plasticityType
call plastic_nonlocal_microstructure (Fe,Fp,ip,el)
case (PLASTICITY_PHENOPLUS_ID) plasticityType
call plastic_phenoplus_microstructure(orientations,ipc,ip,el)
call plastic_phenoplus_microstructure(orientations,ipc,ip,el,F0s,Fes,Fps,Tstar_vs)
end select plasticityType
end subroutine constitutive_microstructure

42
src/crystallite.f90
View File

@ -438,7 +438,11 @@ subroutine crystallite_init
call constitutive_microstructure(crystallite_orientation, & ! pass orientation to constitutive module
crystallite_Fe(1:3,1:3,c,i,e), &
crystallite_Fp(1:3,1:3,c,i,e), &
c,i,e) ! update dependent state variables to be consistent with basic states
c,i,e, &
crystallite_F0, &
crystallite_Fe, &
crystallite_Fp, &
crystallite_Tstar_v) ! update dependent state variables to be consistent with basic states
enddo
enddo
enddo
@ -1715,7 +1719,11 @@ subroutine crystallite_integrateStateRK4()
call constitutive_microstructure(crystallite_orientation, &
crystallite_Fe(1:3,1:3,g,i,e), &
crystallite_Fp(1:3,1:3,g,i,e), &
g, i, e) ! update dependent state variables to be consistent with basic states
g, i, e, &
crystallite_F0, &
crystallite_Fe, &
crystallite_Fp, &
crystallite_Tstar_v) ! update dependent state variables to be consistent with basic states
enddo; enddo; enddo
!$OMP ENDDO
@ -2041,7 +2049,11 @@ subroutine crystallite_integrateStateRKCK45()
call constitutive_microstructure(crystallite_orientation, &
crystallite_Fe(1:3,1:3,g,i,e), &
crystallite_Fp(1:3,1:3,g,i,e), &
g, i, e) ! update dependent state variables to be consistent with basic states
g, i, e, &
crystallite_F0, &
crystallite_Fe, &
crystallite_Fp, &
crystallite_Tstar_v) ! update dependent state variables to be consistent with basic states
enddo; enddo; enddo
!$OMP ENDDO
@ -2261,7 +2273,11 @@ subroutine crystallite_integrateStateRKCK45()
call constitutive_microstructure(crystallite_orientation, &
crystallite_Fe(1:3,1:3,g,i,e), &
crystallite_Fp(1:3,1:3,g,i,e), &
g, i, e) ! update dependent state variables to be consistent with basic states
g, i, e, &
crystallite_F0, &
crystallite_Fe, &
crystallite_Fp, &
crystallite_Tstar_v) ! update dependent state variables to be consistent with basic states
enddo; enddo; enddo
!$OMP ENDDO
@ -2496,7 +2512,11 @@ subroutine crystallite_integrateStateAdaptiveEuler()
call constitutive_microstructure(crystallite_orientation, &
crystallite_Fe(1:3,1:3,g,i,e), &
crystallite_Fp(1:3,1:3,g,i,e), &
g, i, e) ! update dependent state variables to be consistent with basic states
g, i, e, &
crystallite_F0, &
crystallite_Fe, &
crystallite_Fp, &
crystallite_Tstar_v) ! update dependent state variables to be consistent with basic states
enddo; enddo; enddo
!$OMP ENDDO
!$OMP END PARALLEL
@ -2840,7 +2860,11 @@ eIter = FEsolving_execElem(1:2)
call constitutive_microstructure(crystallite_orientation, &
crystallite_Fe(1:3,1:3,g,i,e), &
crystallite_Fp(1:3,1:3,g,i,e), &
g, i, e) ! update dependent state variables to be consistent with basic states
g, i, e, &
crystallite_F0, &
crystallite_Fe, &
crystallite_Fp, &
crystallite_Tstar_v) ! update dependent state variables to be consistent with basic states
enddo; enddo; enddo
!$OMP ENDDO
!$OMP END PARALLEL
@ -3085,7 +3109,11 @@ subroutine crystallite_integrateStateFPI()
call constitutive_microstructure(crystallite_orientation, &
crystallite_Fe(1:3,1:3,g,i,e), &
crystallite_Fp(1:3,1:3,g,i,e), &
g, i, e) ! update dependent state variables to be consistent with basic states
g, i, e, &
crystallite_F0, &
crystallite_Fe, &
crystallite_Fp, &
crystallite_Tstar_v) ! update dependent state variables to be consistent with basic states
p = phaseAt(g,i,e)
c = phasememberAt(g,i,e)
plasticState(p)%previousDotState2(:,c) = plasticState(p)%previousDotState(:,c)

2
src/lattice.f90
View File

@ -377,6 +377,7 @@ module lattice
LATTICE_bcc_Ntrans = 0_pInt, & !sum(lattice_bcc_NtransSystem), & !< total # of transformation systems for bcc
LATTICE_bcc_Ncleavage = 9_pInt !sum(lattice_bcc_NcleavageSystem) !< total # of cleavage systems for bcc
real(pReal), dimension(3+3,LATTICE_bcc_Nslip), parameter, private :: &
LATTICE_bcc_systemSlip = reshape(real([&
! Slip direction Plane normal
@ -835,6 +836,7 @@ module lattice
],pReal),[ 4_pInt + 4_pInt,LATTICE_hex_Ncleavage])
!--------------------------------------------------------------------------------------------------
! body centered tetragonal
integer(pInt), dimension(LATTICE_maxNslipFamily), parameter, public :: &

239
src/plastic_phenoplus.f90
View File

@ -557,7 +557,7 @@ subroutine plastic_phenoplus_init(fileUnit)
! allocate state arrays
sizeState = plastic_phenoplus_totalNslip(instance) & ! s_slip
+ plastic_phenoplus_totalNtwin(instance) & ! s_twin
+ 2_pInt & ! sum(gamma) + sum(f)
+ 2_pInt & ! sum(gamma) + sum(twinVolFrac)
+ plastic_phenoplus_totalNslip(instance) & ! accshear_slip
+ plastic_phenoplus_totalNtwin(instance) & ! accshear_twin
+ plastic_phenoplus_totalNslip(instance) ! kappa
@ -568,7 +568,7 @@ subroutine plastic_phenoplus_init(fileUnit)
! memory leak issue.
sizeDotState = plastic_phenoplus_totalNslip(instance) & ! s_slip
+ plastic_phenoplus_totalNtwin(instance) & ! s_twin
+ 2_pInt & ! sum(gamma) + sum(f)
+ 2_pInt & ! sum(gamma) + sum(twinVolFrac)
+ plastic_phenoplus_totalNslip(instance) & ! accshear_slip
+ plastic_phenoplus_totalNtwin(instance) ! accshear_twin
@ -739,9 +739,11 @@ end subroutine plastic_phenoplus_aTolState
!--------------------------------------------------------------------------------------------------
!> @brief calculate push-up factors (kappa) for each voxel based on its neighbors
!--------------------------------------------------------------------------------------------------
subroutine plastic_phenoplus_microstructure(orientation,ipc,ip,el)
subroutine plastic_phenoplus_microstructure(orientation,ipc,ip,el,F0,Fe,Fp,Tstar_v)
use math, only: pi, &
math_identity2nd, &
math_mul33x33, &
math_mul33xx33, &
math_mul3x3, &
math_transpose33, &
math_qDot, &
@ -763,7 +765,11 @@ subroutine plastic_phenoplus_microstructure(orientation,ipc,ip,el)
homogenization_maxNgrains, &
plasticState
use lattice, only: lattice_sn, &
use lattice, only: lattice_Sslip_v, &
lattice_maxNslipFamily, &
lattice_NslipSystem, &
lattice_NslipSystem, &
lattice_sn, &
lattice_sd, &
lattice_qDisorientation
@ -772,56 +778,78 @@ subroutine plastic_phenoplus_microstructure(orientation,ipc,ip,el)
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
el
real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
F0, & !< deformation gradient from last increment
Fe, & !< elastic deformation gradient
Fp !< elastic deformation gradient !< element
real(pReal), dimension(4,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
orientation ! crystal orientation in quaternions
orientation !< crystal orientation in quaternions
real(pReal), dimension(6,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: &
Tstar_v !< for calculation of gdot
!***local variables
integer(pInt) instance, & !my instance of this plasticity
ph, & !my phase
of, & !my spatial position in memory (offset)
textureID, & !my texture
index_myFamily, &
Nneighbors, & !number of neighbors (<= 6)
vld_Nneighbors, & !number of my valid neighbors
n, & !neighbor index (for iterating through all neighbors)
n_calcTaylor, & !
n_phasecheck, & !
ns, & !number of slip system
nt, & !number of twin system
me_slip, & !my slip system index
neighbor_el, & !element number of neighboring material point
neighbor_ip, & !integration point of neighboring material point
neighbor_n, & !I have no idea what is this
neighbor_ipc, & !I have no idea what is this
neighbor_of, & !spatial position in memory for this neighbor (offset)
neighbor_ph, & !neighbor's phase
neighbor_instance, & !neighbor's instance of this plasticity
neighbor_tex, & !neighbor's texture ID
ne_slip_ac, & !loop to find neighbor shear
ne_slip, & !slip system index for neighbor
index_kappa, & !index of pushup factors in plasticState
offset_acshear_slip, & !offset in PlasticState for the accumulative shear
j !quickly loop through slip families
j, & !quickly loop through slip families
f,i,& !loop counter for me
f_ne, i_ne !loop counter for neighbor
real(pReal) kappa_max, & !
tmp_myshear_slip, & !temp storage for accumulative shear for me
mprime_cut, & !m' cutoff to consider neighboring effect
avg_acshear_ne, & !the average accumulative shear from my neighbor
tmp_mprime, & !temp holder for m' value
tmp_acshear !temp holder for accumulative shear for m'
real(pReal) mprime_cut, & !m' cutoff to consider neighboring effect
dtaylor_cut, & !threshold for determine high contrast interface using Taylor factor
tau_slip, & !the average accumulative shear from my neighbor
taylor_me, & !Taylor factor for me
taylor_ne, & !Taylor factor for my current neighbor
d_vonstrain, & !von Mises delta strain (temp container)
sum_gdot !total shear rate for given material point
real(pReal), dimension(3,3) :: &
F0_me, & !my deformation gradient from last converged increment
Fe_me, & !my elastic deformation gradient
Fp_me, & !my plastic deformation gradient
dF_me, & !my deformation gradient change (delta)
dE_me, & !my Green Lagrangian strain tensor (delta)
F0_ne, & !
Fe_ne, & !elastic deformation gradient of my current neighbor
Fp_ne, & !plastic deformation gradient of my current neighbor
dF_ne, & !deformation gradient of my current neighbor
dE_ne !delta Green Lagrangian strain tensor
real(pReal), dimension(plastic_phenoplus_totalNslip(phase_plasticityInstance(material_phase(1,ip,el)))) :: &
m_primes, & !m' between me_alpha(one) and neighbor beta(all)
me_acshear, & !temp storage for ac_shear of one particular system for me
ne_acshear !temp storage for ac_shear of one particular system for one of my neighbor
m_primes !m' between me_alpha(one) and neighbor beta(all)
real(pReal), dimension(3,plastic_phenoplus_totalNslip(phase_plasticityInstance(material_phase(1,ip,el)))) :: &
slipNormal, &
slipDirect
real(pReal), dimension(4) :: my_orientation, & !store my orientation
real(pReal), dimension(4) :: &
my_orientation, & !store my orientation
neighbor_orientation, & !store my neighbor orientation
absMisorientation
real(pReal), dimension(FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,el))))) :: &
ne_mprimes !m' between each neighbor
ne_mprimes, & !m' between each neighbor
d_taylors !store (taylor_ne-taylor_me) for each neighbor
!***Get my properties
Nneighbors = FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,el))))
@ -831,91 +859,136 @@ subroutine plastic_phenoplus_microstructure(orientation,ipc,ip,el)
instance = phase_plasticityInstance(ph) !get my instance based on phase ID
ns = plastic_phenoplus_totalNslip(instance)
nt = plastic_phenoplus_totalNtwin(instance)
offset_acshear_slip = ns + nt + 2_pInt
index_kappa = ns + nt + 2_pInt + ns + nt !location of kappa in plasticState
!***init calculation for given voxel
mprime_cut = 0.7_pReal !set by Dr.Bieler
dtaylor_cut = 1.0_pReal !set by Chen, quick test only
!***gather my accumulative shear from palsticState
FINDMYSHEAR: do j = 1_pInt,ns
me_acshear(j) = plasticState(ph)%state(offset_acshear_slip+j, of)
enddo FINDMYSHEAR
!***gather my orientation and slip systems
!***gather my orientation, F and slip systems
my_orientation = orientation(1:4, ipc, ip, el)
F0_me = F0(1:3, 1:3, ipc, ip, el)
Fe_me = Fe(1:3, 1:3, ipc, ip, el)
Fp_me = Fp(1:3, 1:3, ipc, ip, el)
slipNormal(1:3, 1:ns) = lattice_sn(1:3, 1:ns, ph)
slipDirect(1:3, 1:ns) = lattice_sd(1:3, 1:ns, ph)
kappa_max = plastic_phenoplus_kappa_max(instance) !maximum pushups allowed (READIN)
!***calculate kappa between me and all my neighbors
LOOPMYSLIP: DO me_slip=1_pInt,ns
vld_Nneighbors = Nneighbors
tmp_myshear_slip = me_acshear(me_slip)
tmp_mprime = 0.0_pReal !highest m' from all neighbors
tmp_acshear = 0.0_pReal !accumulative shear from highest m'
!***check if all my neighbors have the same phase as me
vld_Nneighbors = 0
PHASECHECK: DO n_phasecheck = 1_pInt, Nneighbors
!******for each of my neighbor
neighbor_el = mesh_ipNeighborhood( 1, n_phasecheck, ip, el )
neighbor_ip = mesh_ipNeighborhood( 2, n_phasecheck, ip, el )
neighbor_ipc = 1
neighbor_of = phasememberAt( neighbor_ipc, neighbor_ip, neighbor_el )
neighbor_ph = phaseAt( neighbor_ipc, neighbor_ip, neighbor_el )
IF (neighbor_ph == ph) THEN
vld_Nneighbors = vld_Nneighbors + 1_pInt
ENDIF
ENDDO PHASECHECK
!***go through my neighbors to find highest m'
LOOPNEIGHBORS: DO n=1_pInt,Nneighbors
neighbor_el = mesh_ipNeighborhood(1,n,ip,el)
neighbor_ip = mesh_ipNeighborhood(2,n,ip,el)
neighbor_n = 1 !It is ipc
neighbor_of = phasememberAt( neighbor_n, neighbor_ip, neighbor_el)
neighbor_ph = phaseAt( neighbor_n, neighbor_ip, neighbor_el)
!***initialize kappa with 1.0 (assume no push-up)
plasticState(ph)%state(index_kappa+1_pInt:index_kappa+ns, of) = 1.0_pReal
!***only calculate kappa for those inside the main phase
IF (vld_Nneighbors == Nneighbors) THEN
!******calculate Taylor factor for me
dF_me = math_mul33x33(Fe_me,Fp_me) - F0_me
dE_me = 0.5*(math_mul33x33(math_transpose33(dF_me), dF_me) - math_identity2nd(3)) !dE = 0.5(dF^tdF-I)
d_vonstrain = SQRT(2.0_pReal/3.0_pReal * math_mul33xx33(dE_me, dE_me))
sum_gdot = 0.0_pReal
!go through my slip system to find the sum of gamma_dot
j = 0_pInt
slipFamilies: DO f = 1_pInt,lattice_maxNslipFamily
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,ph)) !at which index starts my family
slipSystems: DO i = 1_pInt,plastic_phenoplus_Nslip(f,instance)
j = j+1_pInt
tau_slip = dot_product(Tstar_v(1:6, ipc, ip, el),lattice_Sslip_v(1:6,1,index_myFamily+i,ph))
sum_gdot = sum_gdot + &
plastic_phenoplus_gdot0_slip(instance)* &
((abs(tau_slip)/(plasticState(ph)%state(j,of))) &
**plastic_phenoplus_n_slip(instance))*sign(1.0_pReal,tau_slip)
ENDDO slipSystems
ENDDO slipFamilies
taylor_me = d_vonstrain/sum_gdot
!***calculate delta_M (Taylor factor) between each neighbor and me
LOOPCALCTAYLOR: DO n_calcTaylor=1_pInt, Nneighbors
!******for each of my neighbor
neighbor_el = mesh_ipNeighborhood( 1, n_calcTaylor, ip, el )
neighbor_ip = mesh_ipNeighborhood( 2, n_calcTaylor, ip, el )
neighbor_ipc = 1 !It is ipc
neighbor_of = phasememberAt( neighbor_ipc, neighbor_ip, neighbor_el )
neighbor_ph = phaseAt( neighbor_ipc, neighbor_ip, neighbor_el )
neighbor_instance = phase_plasticityInstance( neighbor_ph )
neighbor_tex = material_texture( 1,neighbor_ip, neighbor_el )
neighbor_orientation = orientation(1:4, neighbor_n, neighbor_ip, neighbor_el) !ipc is always 1.
neighbor_orientation = orientation( 1:4, neighbor_ipc, neighbor_ip, neighbor_el ) !ipc is always 1.
Fe_ne = Fe( 1:3, 1:3, neighbor_ipc, neighbor_ip, neighbor_el )
Fp_ne = Fp( 1:3, 1:3, neighbor_ipc, neighbor_ip, neighbor_el )
F0_ne = F0( 1:3, 1:3, neighbor_ipc, neighbor_ip, neighbor_el )
!******calculate the Taylor factor
dF_ne = math_mul33x33(Fe_ne, Fp_ne) - F0_ne
dE_ne = 0.5*(math_mul33x33(math_transpose33(dF_ne), dF_ne) - math_identity2nd(3)) !dE = 0.5(dF^tdF-I)
d_vonstrain = SQRT(2.0_pReal/3.0_pReal * math_mul33xx33(dE_ne, dE_ne))
sum_gdot = 0.0_pReal
!go through my neighbor slip system to calculate sum_gdot
j = 0_pInt
slipFamiliesNeighbor: DO f_ne = 1_pInt,lattice_maxNslipFamily
index_myFamily = sum(lattice_NslipSystem(1:f_ne-1_pInt,neighbor_ph)) ! at which index starts my family
slipSystemsNeighbor: DO i_ne = 1_pInt,plastic_phenoplus_Nslip(f_ne,neighbor_instance)
j = j+1_pInt
tau_slip = dot_product(Tstar_v(1:6, neighbor_ipc, neighbor_ip, neighbor_el), &
lattice_Sslip_v(1:6,1,index_myFamily+i_ne,neighbor_ph))
sum_gdot = sum_gdot &
+plastic_phenoplus_gdot0_slip(neighbor_instance) &
*((abs(tau_slip)/(plasticState(neighbor_ph)%state(j,neighbor_of))) &
**plastic_phenoplus_n_slip(neighbor_instance))*sign(1.0_pReal,tau_slip)
ENDDO slipSystemsNeighbor
ENDDO slipFamiliesNeighbor
taylor_ne = d_vonstrain / sum_gdot
!******calculate Taylor difference
d_taylors(n_calcTaylor) = taylor_ne - taylor_me
ENDDO LOOPCALCTAYLOR
!***Only perform necessary calculation if high contrast interface is detected
IF (maxval(d_taylors) > dtaylor_cut) THEN
!*****calculate kappa per slip system base
LOOPMYSLIP: DO me_slip = 1_pInt, ns
ne_mprimes = 0.0_pReal !initialize max m' to 0 for all neighbors
LOOPMYNEIGHBORS: DO n=1_pInt, Nneighbors
!*******only consider neighbor at the high contrast interface
IF (d_taylors(n) > dtaylor_cut) THEN
neighbor_el = mesh_ipNeighborhood( 1, n_calcTaylor, ip, el )
neighbor_ip = mesh_ipNeighborhood( 2, n_calcTaylor, ip, el )
neighbor_ipc = 1 !It is ipc
neighbor_of = phasememberAt( neighbor_ipc, neighbor_ip, neighbor_el )
neighbor_ph = phaseAt( neighbor_ipc, neighbor_ip, neighbor_el )
neighbor_instance = phase_plasticityInstance( neighbor_ph )
neighbor_tex = material_texture( 1,neighbor_ip, neighbor_el )
neighbor_orientation = orientation( 1:4, neighbor_ipc, neighbor_ip, neighbor_el ) !ipc is always 1.
absMisorientation = lattice_qDisorientation( my_orientation, &
neighbor_orientation, &
0_pInt ) !no need for explicit calculation of symmetry
!***find the accumulative shear for this neighbor
LOOPFINDNEISHEAR: DO ne_slip_ac=1_pInt, ns
ne_acshear(ne_slip_ac) = plasticState(ph)%state(offset_acshear_slip+ne_slip_ac, &
neighbor_of)
ENDDO LOOPFINDNEISHEAR
!***calculate the average accumulative shear and use it as cutoff
avg_acshear_ne = SUM(ne_acshear)/ns
!***
IF (ph==neighbor_ph) THEN
!***walk through all the
!*********go through neighbor slip system to calculate m'
LOOPNEIGHBORSLIP: DO ne_slip=1_pInt,ns
!***only consider slip system that is active (above average accumulative shear)
IF (ne_acshear(ne_slip) > avg_acshear_ne) THEN
m_primes(ne_slip) = abs(math_mul3x3(slipNormal(1:3,me_slip), &
math_qRot(absMisorientation, slipNormal(1:3,ne_slip)))) &
*abs(math_mul3x3(slipDirect(1:3,me_slip), &
math_qRot(absMisorientation, slipDirect(1:3,ne_slip))))
!***find the highest m' and corresponding accumulative shear
IF (m_primes(ne_slip) > tmp_mprime) THEN
tmp_mprime = m_primes(ne_slip)
tmp_acshear = ne_acshear(ne_slip)
ENDIF
ENDIF
ENDDO LOOPNEIGHBORSLIP
ELSE
ne_mprimes(n) = 0.0_pReal
vld_Nneighbors = vld_Nneighbors - 1_pInt
ne_mprimes(n) = maxval(m_primes)
ENDIF
ENDDO LOOPMYNEIGHBORS
ENDDO LOOPNEIGHBORS
!***check if this element close to rim
IF (vld_Nneighbors < Nneighbors) THEN
!***rim voxel, no modification allowed
plasticState(ph)%state(index_kappa+me_slip, of) = 1.0_pReal
ELSE
!***patch voxel, started to calculate push up factor for gamma_dot
IF ((tmp_mprime > mprime_cut) .AND. (tmp_acshear > tmp_myshear_slip)) THEN
plasticState(ph)%state(index_kappa+me_slip, of) = 1.0_pReal / tmp_mprime
ELSE
!***minimum damping factor is 0.5
plasticState(ph)%state(index_kappa+me_slip, of) = 0.5_pReal + tmp_mprime * 0.5_pReal
ENDIF
ENDIF
plasticState(ph)%state(index_kappa+me_slip, of) = &
1.03_pReal + 0.03_pReal*ERF(4.0_pReal * maxval(ne_mprimes) - 4.0_pReal)
ENDDO LOOPMYSLIP
ENDIF
ENDIF
end subroutine plastic_phenoplus_microstructure