* added a new material parameter "surfaceTransmissivity" (default value 1.0) which allows to change the transmissivity of the material surface between 0 and 1

* now complaining when encountering an unknown nonlocal parameter in material.config
* use same error ID for all material parameters out of bounds
* symmetric flux calculation in side dotState can now be omitted (because of new treatment of periodicity)
* switching back to "local flux balance" (add leaving and entering fluxes at central MP, don't touch neighbor) instead of "flux distribution" (subtract leaving fluxes from central MP and add them at neighboring MP). This has the advantage that there is almost no need for CRITICAL statements in parallelization, so hopefully this results in some speed up.
This commit is contained in:
Christoph Kords 2011-02-16 16:35:38 +00:00
parent 8f626c8989
commit 24d33bf2ff
2 changed files with 170 additions and 148 deletions

View File

@ -1216,19 +1216,13 @@ endfunction
case (231) case (231)
msg = 'Non-positive prefactor for self-diffusion coefficient' msg = 'Non-positive prefactor for self-diffusion coefficient'
case (232) case (232)
msg = 'Non-positive activation energy for dislocation climb' msg = 'Non-positive activation energy'
case (233) case (233)
msg = 'Non-positive relevant dislocation density' msg = 'Non-positive relevant dislocation density'
case (234)
msg = 'Negative cutoff radius'
case (235) case (235)
msg = 'Non-positive attack frequency' msg = 'Material parameter in nonlocal constitutive phase out of bounds:'
case (236) case (236)
msg = 'Non-positive wall depth' msg = 'Unknown material parameter in nonlocal constitutive phase:'
case (237)
msg = 'Non-positive obstacle strength'
case (238)
msg = 'Negative scatter for dislocation density'
case (240) case (240)
msg = 'Non-positive Taylor factor' msg = 'Non-positive Taylor factor'
case (241) case (241)

View File

@ -63,7 +63,8 @@ real(pReal), dimension(:), allocatable :: constitutive_nonlocal_
constitutive_nonlocal_tauObs, & ! obstacle strength in Pa constitutive_nonlocal_tauObs, & ! obstacle strength in Pa
constitutive_nonlocal_fattack, & ! attack frequency in Hz constitutive_nonlocal_fattack, & ! attack frequency in Hz
constitutive_nonlocal_vs, & ! maximum dislocation velocity = velocity of sound constitutive_nonlocal_vs, & ! maximum dislocation velocity = velocity of sound
constitutive_nonlocal_rhoSglScatter ! standard deviation of scatter in initial dislocation density constitutive_nonlocal_rhoSglScatter, & ! standard deviation of scatter in initial dislocation density
constitutive_nonlocal_surfaceTransmissivity ! transmissivity at free surface
real(pReal), dimension(:,:,:), allocatable :: constitutive_nonlocal_Cslip_66 ! elasticity matrix in Mandel notation for each instance real(pReal), dimension(:,:,:), allocatable :: constitutive_nonlocal_Cslip_66 ! elasticity matrix in Mandel notation for each instance
real(pReal), dimension(:,:,:,:,:), allocatable :: constitutive_nonlocal_Cslip_3333 ! elasticity matrix for each instance real(pReal), dimension(:,:,:,:,:), allocatable :: constitutive_nonlocal_Cslip_3333 ! elasticity matrix for each instance
real(pReal), dimension(:,:), allocatable :: constitutive_nonlocal_rhoSglEdgePos0, & ! initial edge_pos dislocation density per slip system for each family and instance real(pReal), dimension(:,:), allocatable :: constitutive_nonlocal_rhoSglEdgePos0, & ! initial edge_pos dislocation density per slip system for each family and instance
@ -224,6 +225,7 @@ allocate(constitutive_nonlocal_tauObs(maxNinstance))
allocate(constitutive_nonlocal_vs(maxNinstance)) allocate(constitutive_nonlocal_vs(maxNinstance))
allocate(constitutive_nonlocal_fattack(maxNinstance)) allocate(constitutive_nonlocal_fattack(maxNinstance))
allocate(constitutive_nonlocal_rhoSglScatter(maxNinstance)) allocate(constitutive_nonlocal_rhoSglScatter(maxNinstance))
allocate(constitutive_nonlocal_surfaceTransmissivity(maxNinstance))
constitutive_nonlocal_CoverA = 0.0_pReal constitutive_nonlocal_CoverA = 0.0_pReal
constitutive_nonlocal_C11 = 0.0_pReal constitutive_nonlocal_C11 = 0.0_pReal
constitutive_nonlocal_C12 = 0.0_pReal constitutive_nonlocal_C12 = 0.0_pReal
@ -244,6 +246,7 @@ constitutive_nonlocal_tauObs = 0.0_pReal
constitutive_nonlocal_vs = 0.0_pReal constitutive_nonlocal_vs = 0.0_pReal
constitutive_nonlocal_fattack = 0.0_pReal constitutive_nonlocal_fattack = 0.0_pReal
constitutive_nonlocal_rhoSglScatter = 0.0_pReal constitutive_nonlocal_rhoSglScatter = 0.0_pReal
constitutive_nonlocal_surfaceTransmissivity = 1.0_pReal
allocate(constitutive_nonlocal_rhoSglEdgePos0(lattice_maxNslipFamily, maxNinstance)) allocate(constitutive_nonlocal_rhoSglEdgePos0(lattice_maxNslipFamily, maxNinstance))
allocate(constitutive_nonlocal_rhoSglEdgeNeg0(lattice_maxNslipFamily, maxNinstance)) allocate(constitutive_nonlocal_rhoSglEdgeNeg0(lattice_maxNslipFamily, maxNinstance))
@ -285,12 +288,15 @@ do
if (IO_getTag(line,'[',']') /= '') then ! next section if (IO_getTag(line,'[',']') /= '') then ! next section
section = section + 1 section = section + 1
output = 0 ! reset output counter output = 0 ! reset output counter
cycle
endif endif
if (section > 0 .and. phase_constitution(section) == constitutive_nonlocal_label) then ! one of my sections if (section > 0 .and. phase_constitution(section) == constitutive_nonlocal_label) then ! one of my sections
i = phase_constitutionInstance(section) ! which instance of my constitution is present phase i = phase_constitutionInstance(section) ! which instance of my constitution is present phase
positions = IO_stringPos(line,maxNchunks) positions = IO_stringPos(line,maxNchunks)
tag = IO_lc(IO_stringValue(line,positions,1)) ! extract key tag = IO_lc(IO_stringValue(line,positions,1)) ! extract key
select case(tag) select case(tag)
case('constitution','/nonlocal/')
cycle
case ('(output)') case ('(output)')
output = output + 1 output = output + 1
constitutive_nonlocal_output(output,i) = IO_lc(IO_stringValue(line,positions,2)) constitutive_nonlocal_output(output,i) = IO_lc(IO_stringValue(line,positions,2))
@ -354,6 +360,10 @@ do
constitutive_nonlocal_fattack(i) = IO_floatValue(line,positions,2) constitutive_nonlocal_fattack(i) = IO_floatValue(line,positions,2)
case('rhosglscatter') case('rhosglscatter')
constitutive_nonlocal_rhoSglScatter(i) = IO_floatValue(line,positions,2) constitutive_nonlocal_rhoSglScatter(i) = IO_floatValue(line,positions,2)
case('surfacetransmissivity')
constitutive_nonlocal_surfaceTransmissivity(i) = IO_floatValue(line,positions,2)
case default
call IO_error(236,ext_msg=tag)
end select end select
endif endif
enddo enddo
@ -367,37 +377,39 @@ enddo
!*** sanity checks !*** sanity checks
if (myStructure < 1 .or. myStructure > 3) call IO_error(205) if (myStructure < 1 .or. myStructure > 3) call IO_error(205)
if (sum(constitutive_nonlocal_Nslip(:,i)) <= 0_pInt) call IO_error(225) if (sum(constitutive_nonlocal_Nslip(:,i)) <= 0_pInt) call IO_error(235,ext_msg='Nslip')
do o = 1,maxval(phase_Noutput) do o = 1,maxval(phase_Noutput)
if(len(constitutive_nonlocal_output(o,i)) > 64) call IO_error(666) if(len(constitutive_nonlocal_output(o,i)) > 64) call IO_error(666)
enddo enddo
do f = 1,lattice_maxNslipFamily do f = 1,lattice_maxNslipFamily
if (constitutive_nonlocal_Nslip(f,i) > 0_pInt) then if (constitutive_nonlocal_Nslip(f,i) > 0_pInt) then
if (constitutive_nonlocal_rhoSglEdgePos0(f,i) < 0.0_pReal) call IO_error(220) if (constitutive_nonlocal_rhoSglEdgePos0(f,i) < 0.0_pReal) call IO_error(235,ext_msg='rhoSglEdgePos0')
if (constitutive_nonlocal_rhoSglEdgeNeg0(f,i) < 0.0_pReal) call IO_error(220) if (constitutive_nonlocal_rhoSglEdgeNeg0(f,i) < 0.0_pReal) call IO_error(235,ext_msg='rhoSglEdgeNeg0')
if (constitutive_nonlocal_rhoSglScrewPos0(f,i) < 0.0_pReal) call IO_error(220) if (constitutive_nonlocal_rhoSglScrewPos0(f,i) < 0.0_pReal) call IO_error(235,ext_msg='rhoSglScrewPos0')
if (constitutive_nonlocal_rhoSglScrewNeg0(f,i) < 0.0_pReal) call IO_error(220) if (constitutive_nonlocal_rhoSglScrewNeg0(f,i) < 0.0_pReal) call IO_error(235,ext_msg='rhoSglScrewNeg0')
if (constitutive_nonlocal_rhoDipEdge0(f,i) < 0.0_pReal) call IO_error(220) if (constitutive_nonlocal_rhoDipEdge0(f,i) < 0.0_pReal) call IO_error(235,ext_msg='rhoDipEdge0')
if (constitutive_nonlocal_rhoDipScrew0(f,i) < 0.0_pReal) call IO_error(220) if (constitutive_nonlocal_rhoDipScrew0(f,i) < 0.0_pReal) call IO_error(235,ext_msg='rhoDipScrew0')
if (constitutive_nonlocal_burgersPerSlipFamily(f,i) <= 0.0_pReal) call IO_error(221) if (constitutive_nonlocal_burgersPerSlipFamily(f,i) <= 0.0_pReal) call IO_error(235,ext_msg='burgers')
if (constitutive_nonlocal_lambda0PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(227) if (constitutive_nonlocal_lambda0PerSlipFamily(f,i) <= 0.0_pReal) call IO_error(235,ext_msg='lambda0')
if (constitutive_nonlocal_dLowerEdgePerSlipFamily(f,i) <= 0.0_pReal) call IO_error(228) if (constitutive_nonlocal_dLowerEdgePerSlipFamily(f,i) <= 0.0_pReal) call IO_error(235,ext_msg='dDipMinEdge')
if (constitutive_nonlocal_dLowerScrewPerSlipFamily(f,i) <= 0.0_pReal) call IO_error(228) if (constitutive_nonlocal_dLowerScrewPerSlipFamily(f,i) <= 0.0_pReal) call IO_error(235,ext_msg='dDipMinScrew')
endif endif
enddo enddo
if (any(constitutive_nonlocal_interactionSlipSlip(1:maxval(lattice_interactionSlipSlip(:,:,myStructure)),i) < 0.0_pReal)) & if (any(constitutive_nonlocal_interactionSlipSlip(1:maxval(lattice_interactionSlipSlip(:,:,myStructure)),i) < 0.0_pReal)) &
call IO_error(229) call IO_error(235,ext_msg='interaction_SlipSlip')
if (constitutive_nonlocal_R(i) < 0.0_pReal) call IO_error(234) if (constitutive_nonlocal_R(i) < 0.0_pReal) call IO_error(235,ext_msg='r')
if (constitutive_nonlocal_atomicVolume(i) <= 0.0_pReal) call IO_error(230) if (constitutive_nonlocal_atomicVolume(i) <= 0.0_pReal) call IO_error(235,ext_msg='atomicVolume')
if (constitutive_nonlocal_Dsd0(i) <= 0.0_pReal) call IO_error(231) if (constitutive_nonlocal_Dsd0(i) <= 0.0_pReal) call IO_error(235,ext_msg='Dsd0')
if (constitutive_nonlocal_Qsd(i) <= 0.0_pReal) call IO_error(232) if (constitutive_nonlocal_Qsd(i) <= 0.0_pReal) call IO_error(235,ext_msg='Qsd')
if (constitutive_nonlocal_aTolRho(i) <= 0.0_pReal) call IO_error(233) if (constitutive_nonlocal_aTolRho(i) <= 0.0_pReal) call IO_error(235,ext_msg='aTol_rho')
if (constitutive_nonlocal_d0(i) <= 0.0_pReal) call IO_error(236) if (constitutive_nonlocal_d0(i) <= 0.0_pReal) call IO_error(235,ext_msg='d0')
if (constitutive_nonlocal_tauObs(i) <= 0.0_pReal) call IO_error(237) if (constitutive_nonlocal_tauObs(i) <= 0.0_pReal) call IO_error(235,ext_msg='tauObs')
if (constitutive_nonlocal_vs(i) <= 0.0_pReal) call IO_error(226) if (constitutive_nonlocal_vs(i) <= 0.0_pReal) call IO_error(235,ext_msg='vs')
if (constitutive_nonlocal_fattack(i) <= 0.0_pReal) call IO_error(235) if (constitutive_nonlocal_fattack(i) <= 0.0_pReal) call IO_error(235,ext_msg='fAttack')
if (constitutive_nonlocal_rhoSglScatter(i) < 0.0_pReal) call IO_error(238) if (constitutive_nonlocal_rhoSglScatter(i) < 0.0_pReal) call IO_error(235,ext_msg='rhoSglScatter')
if (constitutive_nonlocal_surfaceTransmissivity(i) < 0.0_pReal &
.or. constitutive_nonlocal_surfaceTransmissivity(i) > 1.0_pReal) call IO_error(235,ext_msg='surfaceTransmissivity')
!*** determine total number of active slip systems !*** determine total number of active slip systems
@ -1356,6 +1368,7 @@ integer(pInt) myInstance, & ! current
n, & ! index of my current neighbor n, & ! index of my current neighbor
neighboring_el, & ! element number of my neighbor neighboring_el, & ! element number of my neighbor
neighboring_ip, & ! integration point of my neighbor neighboring_ip, & ! integration point of my neighbor
neighboring_n, & ! neighbor index pointing to me when looking from my neighbor
opposite_n, & ! index of my opposite neighbor opposite_n, & ! index of my opposite neighbor
opposite_ip, & ! ip of my opposite neighbor opposite_ip, & ! ip of my opposite neighbor
opposite_el, & ! element index of my opposite neighbor opposite_el, & ! element index of my opposite neighbor
@ -1365,52 +1378,57 @@ integer(pInt) myInstance, & ! current
sLattice, & ! index of my current slip system according to lattice order sLattice, & ! index of my current slip system according to lattice order
i i
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),10) :: & real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),10) :: &
rhoDot, & ! density evolution rhoDot, & ! density evolution
rhoDotRemobilization, & ! density evolution by remobilization rhoDotRemobilization, & ! density evolution by remobilization
rhoDotMultiplication, & ! density evolution by multiplication rhoDotMultiplication, & ! density evolution by multiplication
rhoDotFlux, & ! density evolution by flux rhoDotFlux, & ! density evolution by flux
neighboring_rhoDotFlux, & ! density evolution by flux at neighbor neighboring_rhoDotFlux, & ! density evolution by flux at neighbor
rhoDotSingle2DipoleGlide, & ! density evolution by dipole formation (by glide) rhoDotSingle2DipoleGlide, & ! density evolution by dipole formation (by glide)
rhoDotAthermalAnnihilation, & ! density evolution by athermal annihilation rhoDotAthermalAnnihilation, & ! density evolution by athermal annihilation
rhoDotThermalAnnihilation, & ! density evolution by thermal annihilation rhoDotThermalAnnihilation, & ! density evolution by thermal annihilation
rhoDotDipole2SingleStressChange, & ! density evolution by dipole dissociation (by stress increase) rhoDotDipole2SingleStressChange, & ! density evolution by dipole dissociation (by stress increase)
rhoDotSingle2DipoleStressChange ! density evolution by dipole formation (by stress decrease) rhoDotSingle2DipoleStressChange ! density evolution by dipole formation (by stress decrease)
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),8) :: & real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),8) :: &
rhoSgl, & ! current single dislocation densities (positive/negative screw and edge without dipoles) rhoSgl, & ! current single dislocation densities (positive/negative screw and edge without dipoles)
previousRhoSgl ! previous single dislocation densities (positive/negative screw and edge without dipoles) previousRhoSgl ! previous single dislocation densities (positive/negative screw and edge without dipoles)
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),4) :: & real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),4) :: &
fluxdensity, & ! flux density at central material point fluxdensity, & ! flux density at central material point
gdot ! shear rates neighboring_fluxdensity, & ! flux density at neighboring material point
gdot ! shear rates
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: & real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el)))) :: &
rhoForest, & ! forest dislocation density rhoForest, & ! forest dislocation density
tauThreshold, & ! threshold shear stress tauThreshold, & ! threshold shear stress
tau, & ! current resolved shear stress tau, & ! current resolved shear stress
previousTau, & ! previous resolved shear stress previousTau, & ! previous resolved shear stress
invLambda, & ! inverse of mean free path for dislocations invLambda, & ! inverse of mean free path for dislocations
vClimb ! climb velocity of edge dipoles vClimb ! climb velocity of edge dipoles
real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),2) :: & real(pReal), dimension(constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),2) :: &
rhoDip, & ! current dipole dislocation densities (screw and edge dipoles) rhoDip, & ! current dipole dislocation densities (screw and edge dipoles)
previousRhoDip, & ! previous dipole dislocation densities (screw and edge dipoles) previousRhoDip, & ! previous dipole dislocation densities (screw and edge dipoles)
dLower, & ! minimum stable dipole distance for edges and screws dLower, & ! minimum stable dipole distance for edges and screws
dUpper, & ! current maximum stable dipole distance for edges and screws dUpper, & ! current maximum stable dipole distance for edges and screws
previousDUpper, & ! previous maximum stable dipole distance for edges and screws previousDUpper, & ! previous maximum stable dipole distance for edges and screws
dUpperDot ! rate of change of the maximum stable dipole distance for edges and screws dUpperDot ! rate of change of the maximum stable dipole distance for edges and screws
real(pReal), dimension(3,constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),4) :: & real(pReal), dimension(3,constitutive_nonlocal_totalNslip(phase_constitutionInstance(material_phase(g,ip,el))),4) :: &
m ! direction of dislocation motion m ! direction of dislocation motion
real(pReal), dimension(3,3) :: F, & ! total deformation gradient real(pReal), dimension(3,3) :: my_F, & ! my total deformation gradient
neighboring_F, & ! total deformation gradient of my neighbor neighboring_F, & ! total deformation gradient of my neighbor
Favg ! average total deformation gradient of me and my neighbor my_Fe, & ! my elastic deformation gradient
real(pReal), dimension(6) :: Tdislocation_v, & ! current dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress neighboring_Fe, & ! elastic deformation gradient of my neighbor
previousTdislocation_v ! previous dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress Favg ! average total deformation gradient of me and my neighbor
real(pReal), dimension(3) :: surfaceNormal, & ! surface normal in lattice configuration real(pReal), dimension(6) :: Tdislocation_v, & ! current dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress
surfaceNormal_currentconf ! surface normal in current configuration previousTdislocation_v ! previous dislocation stress (resulting from the neighboring excess dislocation densities) as 2nd Piola-Kirchhoff stress
real(pReal) area, & ! area of the current interface real(pReal), dimension(3) :: normal_neighbor2me, & ! interface normal pointing from my neighbor to me in neighbor's lattice configuration
detFe, & ! determinant of elastic defornmation gradient normal_neighbor2me_defConf, & ! interface normal pointing from my neighbor to me in shared deformed configuration
transmissivity, & ! overall transmissivity of dislocation flux to neighboring material point normal_me2neighbor, & ! interface normal pointing from me to my neighbor in my lattice configuration
lineLength, & ! dislocation line length leaving the current interface normal_me2neighbor_defConf ! interface normal pointing from me to my neighbor in shared deformed configuration
D, & ! self diffusion real(pReal) area, & ! area of the current interface
transmissivity, & ! overall transmissivity of dislocation flux to neighboring material point
lineLength, & ! dislocation line length leaving the current interface
D, & ! self diffusion
correction correction
logical, dimension(3) :: periodicSurfaceFlux ! flag indicating periodic fluxes at surfaces when surface normal points mainly in x, y and z direction respectively (in reference configuration) logical considerEnteringFlux, &
considerLeavingFlux
if (verboseDebugger .and. (debug_g==g .and. debug_i==ip .and. debug_e==el)) then if (verboseDebugger .and. (debug_g==g .and. debug_i==ip .and. debug_e==el)) then
!$OMP CRITICAL (write2out) !$OMP CRITICAL (write2out)
@ -1541,7 +1559,6 @@ where (rhoSgl(1:ns,1:2) > 0.0_pReal) &
!*** calculate dislocation fluxes (only for nonlocal constitution) !*** calculate dislocation fluxes (only for nonlocal constitution)
rhoDotFlux = 0.0_pReal rhoDotFlux = 0.0_pReal
periodicSurfaceFlux = (/.true.,.false.,.true./)
if (.not. phase_localConstitution(material_phase(g,ip,el))) then ! only for nonlocal constitution if (.not. phase_localConstitution(material_phase(g,ip,el))) then ! only for nonlocal constitution
@ -1550,100 +1567,111 @@ if (.not. phase_localConstitution(material_phase(g,ip,el))) then
m(1:3,1:ns,3) = lattice_st(1:3, constitutive_nonlocal_slipSystemLattice(1:ns,myInstance), myStructure) m(1:3,1:ns,3) = lattice_st(1:3, constitutive_nonlocal_slipSystemLattice(1:ns,myInstance), myStructure)
m(1:3,1:ns,4) = -lattice_st(1:3, constitutive_nonlocal_slipSystemLattice(1:ns,myInstance), myStructure) m(1:3,1:ns,4) = -lattice_st(1:3, constitutive_nonlocal_slipSystemLattice(1:ns,myInstance), myStructure)
F = math_mul33x33(Fe(1:3,1:3,g,ip,el), Fp(1:3,1:3,g,ip,el)) my_Fe = Fe(1:3,1:3,g,ip,el)
detFe = math_det3x3(Fe(1:3,1:3,g,ip,el)) my_F = math_mul33x33(my_Fe, Fp(1:3,1:3,g,ip,el))
fluxdensity = rhoSgl(1:ns,1:4) * constitutive_nonlocal_v(1:ns,1:4,g,ip,el) fluxdensity = rhoSgl(1:ns,1:4) * constitutive_nonlocal_v(1:ns,1:4,g,ip,el)
!if ((debug_g==g .and. debug_i==ip .and. debug_e==el)) write(6,*) '--> dislocation flux <---'
do n = 1,FE_NipNeighbors(mesh_element(2,el)) ! loop through my neighbors do n = 1,FE_NipNeighbors(mesh_element(2,el)) ! loop through my neighbors
neighboring_el = mesh_ipNeighborhood(1,n,ip,el) neighboring_el = mesh_ipNeighborhood(1,n,ip,el)
neighboring_ip = mesh_ipNeighborhood(2,n,ip,el) neighboring_ip = mesh_ipNeighborhood(2,n,ip,el)
if (neighboring_el > 0_pInt .and. neighboring_ip > 0_pInt) then ! if neighbor exists ...
do neighboring_n = 1,FE_NipNeighbors(mesh_element(2,neighboring_el)) ! find neighboring index that points from my neighbor to myself
if ( el == mesh_ipNeighborhood(1,neighboring_n,neighboring_ip,neighboring_el) &
.and. ip == mesh_ipNeighborhood(2,neighboring_n,neighboring_ip,neighboring_el) ) &
exit
enddo
endif
opposite_n = n + mod(n,2) - mod(n+1,2) opposite_n = n + mod(n,2) - mod(n+1,2)
opposite_el = mesh_ipNeighborhood(1,opposite_n,ip,el) opposite_el = mesh_ipNeighborhood(1,opposite_n,ip,el)
opposite_ip = mesh_ipNeighborhood(2,opposite_n,ip,el) opposite_ip = mesh_ipNeighborhood(2,opposite_n,ip,el)
if ( neighboring_el > 0_pInt .and. neighboring_ip > 0_pInt ) then ! if neighbor exists, average deformation gradient if (neighboring_el > 0_pInt .and. neighboring_ip > 0_pInt) then ! if neighbor exists, average deformation gradient
neighboring_F = math_mul33x33(Fe(1:3,1:3,g,neighboring_ip,neighboring_el), Fp(1:3,1:3,g,neighboring_ip,neighboring_el)) neighboring_Fe = Fe(1:3,1:3,g,neighboring_ip,neighboring_el)
Favg = 0.5_pReal * (F + neighboring_F) neighboring_F = math_mul33x33(neighboring_Fe, Fp(1:3,1:3,g,neighboring_ip,neighboring_el))
Favg = 0.5_pReal * (my_F + neighboring_F)
else ! if no neighbor, take my value as average else ! if no neighbor, take my value as average
Favg = F Favg = my_F
endif endif
surfaceNormal_currentconf = math_det3x3(Favg) * math_mul33x3(math_inv3x3(transpose(Favg)), mesh_ipAreaNormal(1:3,n,ip,el)) ! calculate the normal of the interface in current ...
surfaceNormal = math_mul33x3(transpose(Fe(1:3,1:3,g,ip,el)), surfaceNormal_currentconf) / detFe ! ... and lattice configuration !* FLUX FROM ME TO MY NEIGHBOR
area = mesh_ipArea(n,ip,el) * math_norm3(surfaceNormal) !* This is not considered, if my opposite neighbor has a local constitution.
surfaceNormal = surfaceNormal / math_norm3(surfaceNormal) ! normalize the surface normal to unit length !* Then, we assume, that the opposite(!) neighbor sends an equal amount of dislocations to me.
!* So the net flux in the direction of my neighbor is equal to zero:
!* leaving flux to neighbor == entering flux from opposite neighbor
!* In case of reduced transmissivity, part of the leaving flux is stored as dead dislocation density.
!* That means for an interface of zero transmissivity the leaving flux is fully converted to dead dislocations.
neighboring_rhoDotFlux = 0.0_pReal considerLeavingFlux = .true.
! if ((debug_g==g .and. debug_i==ip .and. debug_e==el)) write(6,'(a,x,i2)') 'neighbor',n if (opposite_el > 0 .and. opposite_ip > 0) then
do s = 1,ns if (phase_localConstitution(material_phase(1,opposite_ip,opposite_el))) &
! if ((debug_g==g .and. debug_i==ip .and. debug_e==el)) write(6,'(a,x,i2)') ' system',s considerLeavingFlux = .false.
do t = 1,4
! if ((debug_g==g .and. debug_i==ip .and. debug_e==el)) write(6,'(a,x,i2)') ' type',t
c = (t + 1) / 2
topp = t + mod(t,2) - mod(t+1,2)
if ( abs(math_mul3x3(m(1:3,s,t),surfaceNormal)) > 0.01_pReal &
.and. fluxdensity(s,t) * math_mul3x3(m(1:3,s,t),surfaceNormal) > 0.0_pReal ) then ! outgoing flux
lineLength = fluxdensity(s,t) * math_mul3x3(m(1:3,s,t),surfaceNormal) * area ! line length that wants to leave thrugh this interface
if (opposite_el > 0 .and. opposite_ip > 0) then ! opposite neighbor is present...
if (.not. phase_localConstitution(material_phase(1,opposite_ip,opposite_el))) then ! and is of nonlocal constitution (if it's not, then we assume, that the neighbor sends an equal amount of dislocations)...
rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / mesh_ipVolume(ip,el) ! subtract dislocation flux from current mobile type
! if ((debug_g==g .and. debug_i==ip .and. debug_e==el)) write(6,'(a,x,e12.5)') ' outgoing flux:', lineLength / mesh_ipVolume(ip,el)
endif
else ! if free surface on opposite surface...
if (.not. all(periodicSurfaceFlux(maxloc(abs(mesh_ipAreaNormal(1:3,opposite_n,ip,el))))) ) then ! has no enforced symmetry...
rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / mesh_ipVolume(ip,el) ! subtract dislocation flux from current mobile type
! if ((debug_g==g .and. debug_i==ip .and. debug_e==el)) write(6,'(a,x,e12.5)') ' outgoing flux:', lineLength / mesh_ipVolume(ip,el)
endif
endif
rhoDotFlux(s,t+4) = rhoDotFlux(s,t+4) + lineLength / mesh_ipVolume(ip,el) &
* (1.0_pReal - sum(constitutive_nonlocal_compatibility(c,1:ns,s,n,ip,el)**2.0_pReal)) &
* sign(1.0_pReal, fluxdensity(s,t)) ! dislocation flux that is not able to leave through interface (because of low transmissivity) will remain as immobile single density at the material point
if (neighboring_el > 0 .and. neighboring_ip > 0) then ! neighbor present...
if ( .not. phase_localConstitution(material_phase(1,neighboring_ip,neighboring_el))) then ! and is of nonlocal constitution
where (constitutive_nonlocal_compatibility(c,1:ns,s,n,ip,el) > 0.0_pReal) & ! ..positive compatibility
neighboring_rhoDotFlux(1:ns,t) = neighboring_rhoDotFlux(1:ns,t) & ! ....transferring to equally signed dislocation type at neighbor
+ lineLength / mesh_ipVolume(neighboring_ip,neighboring_el) &
* constitutive_nonlocal_compatibility(c,1:ns,s,n,ip,el) ** 2.0_pReal
where (constitutive_nonlocal_compatibility(c,1:ns,s,n,ip,el) < 0.0_pReal) & ! ..negative compatibility
neighboring_rhoDotFlux(1:ns,topp) = neighboring_rhoDotFlux(1:ns,topp) & ! ....transferring to opposite signed dislocation type at neighbor
+ lineLength / mesh_ipVolume(neighboring_ip,neighboring_el) &
* constitutive_nonlocal_compatibility(c,1:ns,s,n,ip,el) ** 2.0_pReal
! if ((debug_g==g .and. debug_i==ip .and. debug_e==el)) write(6,'(a,x,e12.5)') ' entering flux at neighbor:', lineLength / mesh_ipVolume(ip,el) &
! * sum(constitutive_nonlocal_compatibility(c,:,s,n,ip,el) ** 2.0_pReal)
endif
endif
endif
enddo ! dislocation type loop
enddo ! slip system loop
if (any(abs(neighboring_rhoDotFlux) > 10.0_pReal)) then ! only significant density change in neighbor is considered
!$OMP CRITICAL (fluxes)
constitutive_nonlocal_rhoDotFlux(1:ns,1:10,g,neighboring_ip,neighboring_el) = &
constitutive_nonlocal_rhoDotFlux(1:ns,1:10,g,neighboring_ip,neighboring_el) + neighboring_rhoDotFlux
dotState(g,neighboring_ip,neighboring_el)%p(1:10*ns) = &
dotState(g,neighboring_ip,neighboring_el)%p(1:10*ns) + reshape(neighboring_rhoDotFlux,(/10*ns/))
!$OMP END CRITICAL (fluxes)
else
neighboring_rhoDotFlux = 0.0_pReal
endif endif
if (considerLeavingFlux) then
normal_me2neighbor_defConf = math_det3x3(Favg) * math_mul33x3(math_inv3x3(transpose(Favg)), mesh_ipAreaNormal(1:3,n,ip,el)) ! calculate the normal of the interface in (average) deformed configuration (pointing from me to my neighbor!!!)
normal_me2neighbor = math_mul33x3(transpose(my_Fe), normal_me2neighbor_defConf) / math_det3x3(my_Fe) ! interface normal in my lattice configuration
area = mesh_ipArea(n,ip,el) * math_norm3(normal_me2neighbor)
normal_me2neighbor = normal_me2neighbor / math_norm3(normal_me2neighbor) ! normalize the surface normal to unit length
do s = 1,ns
do t = 1,4
c = (t + 1) / 2
if (fluxdensity(s,t) * math_mul3x3(m(1:3,s,t), normal_me2neighbor) > 0.0_pReal ) then ! flux from me to my neighbor == leaving flux for me (might also be a pure flux from my mobile density to dead density if interface not at all transmissive)
lineLength = fluxdensity(s,t) * math_mul3x3(m(1:3,s,t), normal_me2neighbor) * area ! positive line length that wants to leave through this interface
transmissivity = sum(constitutive_nonlocal_compatibility(c,1:ns,s,n,ip,el)**2.0_pReal) ! overall transmissivity from this slip system to my neighbor
rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / mesh_ipVolume(ip,el) ! subtract dislocation flux from current mobile type
rhoDotFlux(s,t+4) = rhoDotFlux(s,t+4) + lineLength / mesh_ipVolume(ip,el) * (1.0_pReal - transmissivity) &
* sign(1.0_pReal, fluxdensity(s,t)) ! dislocation flux that is not able to leave through interface (because of low transmissivity) will remain as immobile single density at the material point
endif
enddo
enddo
endif
!* FLUX FROM MY NEIGHBOR TO ME
!* This is only considered, if I have a neighbor of nonlocal constitution that is at least a little bit compatible.
!* If it's not at all compatible, no flux is arriving, because everything is dammed in front of my neighbor's interface.
!* The entering flux from my neighbor will be distributed on my slip systems according to the compatibility
considerEnteringFlux = .false.
if (neighboring_el > 0_pInt .or. neighboring_ip > 0_pInt) then
if (.not. phase_localConstitution(material_phase(1,neighboring_ip,neighboring_el)) &
.and. any(constitutive_nonlocal_compatibility(:,:,:,n,ip,el) > 0.0_pReal)) &
considerEnteringFlux = .true.
endif
if (considerEnteringFlux) then
forall (t = 1:4) &
neighboring_fluxdensity(1:ns,t) = state(g,neighboring_ip,neighboring_el)%p((t-1)*ns+1:t*ns) &
* constitutive_nonlocal_v(1:ns,t,g,neighboring_ip,neighboring_el)
normal_neighbor2me_defConf = math_det3x3(Favg) &
* math_mul33x3(math_inv3x3(transpose(Favg)), mesh_ipAreaNormal(1:3,neighboring_n,neighboring_ip,neighboring_el)) ! calculate the normal of the interface in (average) deformed configuration (now pointing from my neighbor to me!!!)
normal_neighbor2me = math_mul33x3(transpose(neighboring_Fe), normal_neighbor2me_defConf) / math_det3x3(neighboring_Fe) ! interface normal in the lattice configuration of my neighbor
area = mesh_ipArea(neighboring_n,neighboring_ip,neighboring_el) * math_norm3(normal_neighbor2me)
normal_neighbor2me = normal_neighbor2me / math_norm3(normal_neighbor2me) ! normalize the surface normal to unit length
do s = 1,ns
do t = 1,4
c = (t + 1) / 2
topp = t + mod(t,2) - mod(t+1,2)
if (neighboring_fluxdensity(s,t) * math_mul3x3(m(1:3,s,t), normal_neighbor2me) > 0.0_pReal) then ! flux from my neighbor to me == entering flux for me
lineLength = neighboring_fluxdensity(s,t) * math_mul3x3(m(1:3,s,t), normal_neighbor2me) * area ! positive line length that wants to enter through this interface
where (constitutive_nonlocal_compatibility(c,1:ns,s,n,ip,el) > 0.0_pReal) & ! positive compatibility...
rhoDotFlux(1:ns,t) = rhoDotFlux(1:ns,t) + lineLength / mesh_ipVolume(ip,el) & ! ... transferring to equally signed dislocation type
* constitutive_nonlocal_compatibility(c,1:ns,s,n,ip,el) ** 2.0_pReal
where (constitutive_nonlocal_compatibility(c,1:ns,s,n,ip,el) < 0.0_pReal) & ! ..negative compatibility...
rhoDotFlux(1:ns,topp) = rhoDotFlux(1:ns,topp) + lineLength / mesh_ipVolume(ip,el) & ! ... transferring to opposite signed dislocation type
* constitutive_nonlocal_compatibility(c,1:ns,s,n,ip,el) ** 2.0_pReal
endif
enddo
enddo
endif
enddo ! neighbor loop enddo ! neighbor loop
if (any(abs(rhoDotFlux) > 0.0_pReal)) then constitutive_nonlocal_rhoDotFlux(1:ns,1:10,g,ip,el) = constitutive_nonlocal_rhoDotFlux(1:ns,1:10,g,ip,el) + rhoDotFlux
!$OMP CRITICAL (fluxes)
constitutive_nonlocal_rhoDotFlux(1:ns,1:10,g,ip,el) = constitutive_nonlocal_rhoDotFlux(1:ns,1:10,g,ip,el) + rhoDotFlux
!$OMP END CRITICAL (fluxes)
endif
endif endif