gfortran complaint about implicit casting,floating point comparison, and unused imported variables.

additionally changed use of math_invert33 to math_inv33 if det is not needed
This commit is contained in:
Martin Diehl 2015-04-10 19:09:26 +00:00
parent 5526d2a721
commit d7b36c2c30
6 changed files with 43 additions and 112 deletions

View File

@ -3590,7 +3590,6 @@ logical function crystallite_integrateStress(&
math_mul99x99, & math_mul99x99, &
math_transpose33, & math_transpose33, &
math_inv33, & math_inv33, &
math_invert33, &
math_invert, & math_invert, &
math_det33, & math_det33, &
math_norm33, & math_norm33, &
@ -3654,8 +3653,7 @@ logical function crystallite_integrateStress(&
dLi_dFi3333, & dLi_dFi3333, &
dLp_dT3333, & dLp_dT3333, &
dLi_dT3333 dLi_dT3333
real(pReal) det, & ! determinant real(pReal) detInvFi, & ! determinant of InvFi
detInvFi, &
steplengthLp0, & steplengthLp0, &
steplengthLp, & steplengthLp, &
steplengthLi0, & steplengthLi0, &
@ -3663,7 +3661,6 @@ logical function crystallite_integrateStress(&
dt, & ! time increment dt, & ! time increment
aTolLp, & aTolLp, &
aTolLi aTolLi
logical error ! flag indicating an error
integer(pInt) NiterationStressLp, & ! number of stress integrations integer(pInt) NiterationStressLp, & ! number of stress integrations
NiterationStressLi, & ! number of inner stress integrations NiterationStressLi, & ! number of inner stress integrations
ierr, & ! error indicator for LAPACK ierr, & ! error indicator for LAPACK
@ -3718,7 +3715,7 @@ logical function crystallite_integrateStress(&
!* inversion of Fp_current... !* inversion of Fp_current...
invFp_current = math_inv33(Fp_current) invFp_current = math_inv33(Fp_current)
if (all(invFp_current == 0.0_pReal)) then ! ... failed? if (all(invFp_current <= tiny(0.0_pReal))) then ! math_inv33 returns zero when failed, avoid floating point comparison
#ifndef _OPENMP #ifndef _OPENMP
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3)') '<< CRYST >> integrateStress failed on inversion of Fp_current at el (elFE) ip g ',& write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3)') '<< CRYST >> integrateStress failed on inversion of Fp_current at el (elFE) ip g ',&
@ -3729,12 +3726,12 @@ logical function crystallite_integrateStress(&
#endif #endif
return return
endif endif
A = math_mul33x33(Fg_new,invFp_current) ! intermediate tensor needed later to calculate dFe_dLp A = math_mul33x33(Fg_new,invFp_current) ! intermediate tensor needed later to calculate dFe_dLp
!* inversion of Fi_current... !* inversion of Fi_current...
invFi_current = math_inv33(Fi_current) invFi_current = math_inv33(Fi_current)
if (all(invFi_current == 0.0_pReal)) then ! ... failed? if (all(invFi_current <= tiny(0.0_pReal))) then ! math_inv33 returns zero when failed, avoid floating point comparison
#ifndef _OPENMP #ifndef _OPENMP
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3)') '<< CRYST >> integrateStress failed on inversion of Fi_current at el (elFE) ip g ',& write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3)') '<< CRYST >> integrateStress failed on inversion of Fi_current at el (elFE) ip g ',&
@ -3974,9 +3971,9 @@ logical function crystallite_integrateStress(&
!* calculate new plastic and elastic deformation gradient !* calculate new plastic and elastic deformation gradient
invFp_new = math_mul33x33(invFp_current,B) invFp_new = math_mul33x33(invFp_current,B)
invFp_new = invFp_new / math_det33(invFp_new)**(1.0_pReal/3.0_pReal) ! regularize by det invFp_new = invFp_new / math_det33(invFp_new)**(1.0_pReal/3.0_pReal) ! regularize by det
call math_invert33(invFp_new,Fp_new,det,error) Fp_new = math_inv33(invFp_new)
if (error .or. any(Fp_new /= Fp_new)) then if (all(Fp_new <= tiny(0.0_pReal))) then ! math_inv33 returns zero when failed, avoid floating point comparison
#ifndef _OPENMP #ifndef _OPENMP
if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then if (iand(debug_level(debug_crystallite), debug_levelBasic) /= 0_pInt) then
write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,a,i3)') '<< CRYST >> integrateStress failed on invFp_new inversion at el ip g ',& write(6,'(a,i8,1x,a,i8,a,1x,i2,1x,i3,a,i3)') '<< CRYST >> integrateStress failed on invFp_new inversion at el ip g ',&

View File

@ -1039,8 +1039,7 @@ subroutine homogenization_RGC_stressPenalty(rPen,nMis,avgF,fDef,ip,el,homID)
use constitutive, only: & use constitutive, only: &
constitutive_homogenizedC constitutive_homogenizedC
use math, only: & use math, only: &
math_civita,& math_civita
math_invert33
use material, only: & use material, only: &
homogenization_maxNgrains,& homogenization_maxNgrains,&
homogState, & homogState, &

View File

@ -498,9 +498,8 @@ subroutine mesh_init(ip,el)
#endif #endif
#ifdef Spectral #ifdef Spectral
IO_open_file IO_open_file
use numerics, only: &
divergence_correction
#else #else
use numerics, only: &
IO_open_InputFile IO_open_InputFile
#endif #endif
use debug, only: & use debug, only: &
@ -1547,9 +1546,7 @@ function mesh_regrid(adaptive,resNewInput,minRes)
integer(pInt), dimension(:), allocatable :: indices integer(pInt), dimension(:), allocatable :: indices
real(pReal) :: wgt real(pReal) :: wgt
real(pReal), dimension(3) :: geomSizeNew, geomSize real(pReal), dimension(3) :: geomSizeNew, geomSize
real(pReal), dimension(3,3) :: Favg, Favg_LastInc, & real(pReal), dimension(3,3) :: Favg, Favg_LastInc
FavgNew, Favg_LastIncNew, &
deltaF, deltaF_lastInc
real(pReal), dimension(:,:), allocatable :: & real(pReal), dimension(:,:), allocatable :: &
coordinates, coordinatesNew coordinates, coordinatesNew
real(pReal), dimension(:,:,:), allocatable :: & real(pReal), dimension(:,:,:), allocatable :: &
@ -1559,7 +1556,6 @@ function mesh_regrid(adaptive,resNewInput,minRes)
Tstar, TstarNew, & Tstar, TstarNew, &
stateConst stateConst
real(pReal), dimension(:,:,:,:,:), allocatable :: & real(pReal), dimension(:,:,:,:,:), allocatable :: &
spectralF33, spectralF33New, &
F, FNew, & F, FNew, &
Fp, FpNew, & Fp, FpNew, &
Lp, LpNew, & Lp, LpNew, &

View File

@ -682,8 +682,8 @@ subroutine plastic_dislotwin_init(fileUnit)
if (plastic_dislotwin_Qsd(instance) <= 0.0_pReal) & if (plastic_dislotwin_Qsd(instance) <= 0.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='Qsd ('//PLASTICITY_DISLOTWIN_label//')') call IO_error(211_pInt,el=instance,ext_msg='Qsd ('//PLASTICITY_DISLOTWIN_label//')')
if (sum(plastic_dislotwin_Ntwin(:,instance)) > 0_pInt) then if (sum(plastic_dislotwin_Ntwin(:,instance)) > 0_pInt) then
if (plastic_dislotwin_SFE_0K(instance) == 0.0_pReal .and. & if (plastic_dislotwin_SFE_0K(instance) <= tiny(0.0_pReal) .and. &
plastic_dislotwin_dSFE_dT(instance) == 0.0_pReal .and. & plastic_dislotwin_dSFE_dT(instance) <= tiny(0.0_pReal) .and. &
lattice_structure(phase) == LATTICE_fcc_ID) & lattice_structure(phase) == LATTICE_fcc_ID) &
call IO_error(211_pInt,el=instance,ext_msg='SFE0K ('//PLASTICITY_DISLOTWIN_label//')') call IO_error(211_pInt,el=instance,ext_msg='SFE0K ('//PLASTICITY_DISLOTWIN_label//')')
if (plastic_dislotwin_aTolRho(instance) <= 0.0_pReal) & if (plastic_dislotwin_aTolRho(instance) <= 0.0_pReal) &
@ -708,7 +708,7 @@ subroutine plastic_dislotwin_init(fileUnit)
if (plastic_dislotwin_sbVelocity(instance) > 0.0_pReal .and. & if (plastic_dislotwin_sbVelocity(instance) > 0.0_pReal .and. &
plastic_dislotwin_pShearBand(instance) <= 0.0_pReal) & plastic_dislotwin_pShearBand(instance) <= 0.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='pShearBand ('//PLASTICITY_DISLOTWIN_label//')') call IO_error(211_pInt,el=instance,ext_msg='pShearBand ('//PLASTICITY_DISLOTWIN_label//')')
if (plastic_dislotwin_dipoleFormationFactor(instance) /= 0.0_pReal .and. & if (plastic_dislotwin_dipoleFormationFactor(instance) > tiny(0.0_pReal) .and. &
plastic_dislotwin_dipoleFormationFactor(instance) /= 1.0_pReal) & plastic_dislotwin_dipoleFormationFactor(instance) /= 1.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='dipoleFormationFactor ('//PLASTICITY_DISLOTWIN_label//')') call IO_error(211_pInt,el=instance,ext_msg='dipoleFormationFactor ('//PLASTICITY_DISLOTWIN_label//')')
if (plastic_dislotwin_sbVelocity(instance) > 0.0_pReal .and. & if (plastic_dislotwin_sbVelocity(instance) > 0.0_pReal .and. &
@ -1016,7 +1016,6 @@ subroutine plastic_dislotwin_stateInit(ph,instance)
pi pi
use lattice, only: & use lattice, only: &
lattice_maxNslipFamily, & lattice_maxNslipFamily, &
lattice_structure, &
lattice_mu lattice_mu
use material, only: & use material, only: &
plasticState plasticState
@ -1195,9 +1194,6 @@ function plastic_dislotwin_homogenizedC(ipc,ip,el)
subroutine plastic_dislotwin_microstructure(temperature,ipc,ip,el) subroutine plastic_dislotwin_microstructure(temperature,ipc,ip,el)
use math, only: & use math, only: &
pi pi
use mesh, only: &
mesh_NcpElems, &
mesh_maxNips
use material, only: & use material, only: &
material_phase, & material_phase, &
phase_plasticityInstance, & phase_plasticityInstance, &
@ -1206,9 +1202,7 @@ subroutine plastic_dislotwin_microstructure(temperature,ipc,ip,el)
use lattice, only: & use lattice, only: &
lattice_structure, & lattice_structure, &
lattice_mu, & lattice_mu, &
lattice_nu, & lattice_nu
lattice_maxNslipFamily
implicit none implicit none
integer(pInt), intent(in) :: & integer(pInt), intent(in) :: &
@ -1513,8 +1507,8 @@ subroutine plastic_dislotwin_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,Temperature
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! Shear banding (shearband) part ! Shear banding (shearband) part
if(plastic_dislotwin_sbVelocity(instance) /= 0.0_pReal .and. & if(plastic_dislotwin_sbVelocity(instance) > tiny(0.0_pReal) .and. &
plastic_dislotwin_sbResistance(instance) /= 0.0_pReal) then plastic_dislotwin_sbResistance(instance) > tiny(0.0_pReal)) then
gdot_sb = 0.0_pReal gdot_sb = 0.0_pReal
dgdot_dtausb = 0.0_pReal dgdot_dtausb = 0.0_pReal
call math_spectralDecompositionSym33(math_Mandel6to33(Tstar_v),eigValues,eigVectors, error) call math_spectralDecompositionSym33(math_Mandel6to33(Tstar_v),eigValues,eigVectors, error)
@ -1672,9 +1666,6 @@ subroutine plastic_dislotwin_dotState(Tstar_v,Temperature,ipc,ip,el)
tol_math_check tol_math_check
use math, only: & use math, only: &
pi pi
use mesh, only: &
mesh_NcpElems, &
mesh_maxNips
use material, only: & use material, only: &
material_phase, & material_phase, &
phase_plasticityInstance, & phase_plasticityInstance, &
@ -1696,8 +1687,7 @@ subroutine plastic_dislotwin_dotState(Tstar_v,Temperature,ipc,ip,el)
lattice_fcc_twinNucleationSlipPair, & lattice_fcc_twinNucleationSlipPair, &
lattice_fcc_transNucleationTwinPair, & lattice_fcc_transNucleationTwinPair, &
lattice_fcc_shearCritTrans, & lattice_fcc_shearCritTrans, &
LATTICE_fcc_ID, & LATTICE_fcc_ID
LATTICE_bcc_ID
implicit none implicit none
real(pReal), dimension(6), intent(in):: & real(pReal), dimension(6), intent(in):: &
@ -1802,7 +1792,7 @@ subroutine plastic_dislotwin_dotState(Tstar_v,Temperature,ipc,ip,el)
plastic_dislotwin_CAtomicVolume(instance)*plastic_dislotwin_burgersPerSlipSystem(j,instance)**(3.0_pReal) plastic_dislotwin_CAtomicVolume(instance)*plastic_dislotwin_burgersPerSlipSystem(j,instance)**(3.0_pReal)
VacancyDiffusion = & VacancyDiffusion = &
plastic_dislotwin_D0(instance)*exp(-plastic_dislotwin_Qsd(instance)/(kB*Temperature)) plastic_dislotwin_D0(instance)*exp(-plastic_dislotwin_Qsd(instance)/(kB*Temperature))
if (tau_slip(j) == 0.0_pReal) then if (tau_slip(j) <= tiny(0.0_pReal)) then
DotRhoEdgeDipClimb(j) = 0.0_pReal DotRhoEdgeDipClimb(j) = 0.0_pReal
else else
ClimbVelocity(j) = & ClimbVelocity(j) = &
@ -1912,8 +1902,8 @@ subroutine plastic_dislotwin_dotState(Tstar_v,Temperature,ipc,ip,el)
b = lattice_fcc_transNucleationTwinPair(2,j) b = lattice_fcc_transNucleationTwinPair(2,j)
sa = sign(1_pInt, a) sa = sign(1_pInt, a)
sb = sign(1_pInt, b) sb = sign(1_pInt, b)
ssa = sign(1.0_pReal, shearrate_trans(a)) ssa = int(sign(1.0_pReal, shearrate_trans(a)),pInt)
ssb = sign(1.0_pReal, shearrate_trans(b)) ssb = int(sign(1.0_pReal, shearrate_trans(b)),pInt)
if (sa == ssa .and. sb == ssb) then if (sa == ssa .and. sb == ssb) then
probrate_trans(j) = (abs(shear_trans(a)*shearrate_trans(b)) + abs(shear_trans(b)*shearrate_trans(a))) probrate_trans(j) = (abs(shear_trans(a)*shearrate_trans(b)) + abs(shear_trans(b)*shearrate_trans(a)))
@ -1949,7 +1939,6 @@ function plastic_dislotwin_postResults(Tstar_v,Temperature,ipc,ip,el)
use material, only: & use material, only: &
material_phase, & material_phase, &
phase_plasticityInstance,& phase_plasticityInstance,&
phase_Noutput, &
plasticState, & plasticState, &
mappingConstitutive mappingConstitutive
use lattice, only: & use lattice, only: &
@ -2266,7 +2255,7 @@ function plastic_dislotwin_postResults(Tstar_v,Temperature,ipc,ip,el)
endif endif
!* Stress exponent !* Stress exponent
if (gdot_slip(j)==0.0_pReal) then if (gdot_slip(j)<=tiny(0.0_pReal)) then
plastic_dislotwin_postResults(c+j) = 0.0_pReal plastic_dislotwin_postResults(c+j) = 0.0_pReal
else else
plastic_dislotwin_postResults(c+j) = (tau/gdot_slip(j))*dgdot_dtauslip plastic_dislotwin_postResults(c+j) = (tau/gdot_slip(j))*dgdot_dtauslip

View File

@ -1413,7 +1413,6 @@ use lattice, only: lattice_maxNslipFamily
use math, only: math_sampleGaussVar use math, only: math_sampleGaussVar
use mesh, only: mesh_ipVolume, & use mesh, only: mesh_ipVolume, &
mesh_NcpElems, & mesh_NcpElems, &
mesh_maxNips, &
mesh_element, & mesh_element, &
FE_Nips, & FE_Nips, &
FE_geomtype FE_geomtype
@ -1424,9 +1423,6 @@ use material, only: material_phase, &
material_Nphase, & material_Nphase, &
phase_plasticity ,& phase_plasticity ,&
PLASTICITY_NONLOCAL_ID PLASTICITY_NONLOCAL_ID
use numerics,only: &
numerics_integrator
implicit none implicit none
integer(pInt) :: e, & integer(pInt) :: e, &
@ -1561,20 +1557,16 @@ use math, only: &
math_mul33x3, & math_mul33x3, &
math_mul3x3, & math_mul3x3, &
math_norm3, & math_norm3, &
math_invert33, & math_inv33, &
math_transpose33 math_transpose33
use debug, only: & use debug, only: &
debug_level, & debug_level, &
debug_constitutive, & debug_constitutive, &
debug_levelBasic, &
debug_levelExtensive, & debug_levelExtensive, &
debug_levelSelective, & debug_levelSelective, &
debug_g, &
debug_i, & debug_i, &
debug_e debug_e
use mesh, only: & use mesh, only: &
mesh_NcpElems, &
mesh_maxNips, &
mesh_element, & mesh_element, &
mesh_ipNeighborhood, & mesh_ipNeighborhood, &
mesh_ipCoordinates, & mesh_ipCoordinates, &
@ -1586,7 +1578,6 @@ use mesh, only: &
FE_geomtype, & FE_geomtype, &
FE_celltype FE_celltype
use material, only: & use material, only: &
homogenization_maxNgrains, &
material_phase, & material_phase, &
phase_localPlasticity, & phase_localPlasticity, &
plasticState, & plasticState, &
@ -1629,10 +1620,7 @@ integer(pInt) neighbor_el, & ! element number o
n, & n, &
nRealNeighbors ! number of really existing neighbors nRealNeighbors ! number of really existing neighbors
integer(pInt), dimension(2) :: neighbors integer(pInt), dimension(2) :: neighbors
real(pReal) detFe, & real(pReal) FVsize, &
detFp, &
FVsize, &
temp, &
correction, & correction, &
myRhoForest myRhoForest
real(pReal), dimension(2) :: rhoExcessGradient, & real(pReal), dimension(2) :: rhoExcessGradient, &
@ -1664,7 +1652,6 @@ real(pReal), dimension(2,maxval(totalNslip),mesh_maxNipNeighbors) :: &
neighbor_rhoTotal ! total density at neighboring material point neighbor_rhoTotal ! total density at neighboring material point
real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(1_pInt,ip,el))),2) :: & real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(1_pInt,ip,el))),2) :: &
m ! direction of dislocation motion m ! direction of dislocation motion
logical inversionError
ph = mappingConstitutive(2,1,ip,el) ph = mappingConstitutive(2,1,ip,el)
of = mappingConstitutive(1,1,ip,el) of = mappingConstitutive(1,1,ip,el)
@ -1725,8 +1712,8 @@ forall (s = 1_pInt:ns) &
tauBack = 0.0_pReal tauBack = 0.0_pReal
if (.not. phase_localPlasticity(ph) .and. shortRangeStressCorrection(instance)) then if (.not. phase_localPlasticity(ph) .and. shortRangeStressCorrection(instance)) then
call math_invert33(Fe, invFe, detFe, inversionError) invFe = math_inv33(Fe)
call math_invert33(Fp, invFp, detFp, inversionError) invFp = math_inv33(Fp)
rhoExcess(1,1:ns) = rhoSgl(1:ns,1) - rhoSgl(1:ns,2) rhoExcess(1,1:ns) = rhoSgl(1:ns,1) - rhoSgl(1:ns,2)
rhoExcess(2,1:ns) = rhoSgl(1:ns,3) - rhoSgl(1:ns,4) rhoExcess(2,1:ns) = rhoSgl(1:ns,3) - rhoSgl(1:ns,4)
FVsize = mesh_ipVolume(ip,el) ** (1.0_pReal/3.0_pReal) FVsize = mesh_ipVolume(ip,el) ** (1.0_pReal/3.0_pReal)
@ -1806,10 +1793,9 @@ if (.not. phase_localPlasticity(ph) .and. shortRangeStressCorrection(instance))
rhoExcessDifferences(dir) = neighbor_rhoExcess(c,s,neighbors(1)) & rhoExcessDifferences(dir) = neighbor_rhoExcess(c,s,neighbors(1)) &
- neighbor_rhoExcess(c,s,neighbors(2)) - neighbor_rhoExcess(c,s,neighbors(2))
enddo enddo
call math_invert33(connections,invConnections,temp,inversionError) invConnections = math_inv33(connections)
if (inversionError) then if (all(invConnections <= tiny(0.0_pReal))) & ! check for failed in version (math_inv33 returns 0) and avoid floating point equality comparison
call IO_error(-1_pInt,ext_msg='back stress calculation: inversion error') call IO_error(-1_pInt,ext_msg='back stress calculation: inversion error')
endif
rhoExcessGradient(c) = math_mul3x3(m(1:3,s,c), & rhoExcessGradient(c) = math_mul3x3(m(1:3,s,c), &
math_mul33x3(invConnections,rhoExcessDifferences)) math_mul33x3(invConnections,rhoExcessDifferences))
enddo enddo
@ -1867,10 +1853,8 @@ subroutine plastic_nonlocal_kinetics(v, dv_dtau, dv_dtauNS, tau, tauNS, &
use debug, only: debug_level, & use debug, only: debug_level, &
debug_constitutive, & debug_constitutive, &
debug_levelBasic, &
debug_levelExtensive, & debug_levelExtensive, &
debug_levelSelective, & debug_levelSelective, &
debug_g, &
debug_i, & debug_i, &
debug_e debug_e
use material, only: material_phase, & use material, only: material_phase, &
@ -2029,10 +2013,8 @@ use math, only: math_Plain3333to99, &
math_Mandel6to33 math_Mandel6to33
use debug, only: debug_level, & use debug, only: debug_level, &
debug_constitutive, & debug_constitutive, &
debug_levelBasic, &
debug_levelExtensive, & debug_levelExtensive, &
debug_levelSelective, & debug_levelSelective, &
debug_g, &
debug_i, & debug_i, &
debug_e debug_e
use material, only: material_phase, & use material, only: material_phase, &
@ -2226,7 +2208,6 @@ use debug, only: debug_level, &
debug_levelBasic, & debug_levelBasic, &
debug_levelExtensive, & debug_levelExtensive, &
debug_levelSelective, & debug_levelSelective, &
debug_g, &
debug_i, & debug_i, &
debug_e debug_e
use math, only: pi, & use math, only: pi, &
@ -2234,11 +2215,8 @@ use math, only: pi, &
use lattice, only: lattice_Sslip_v ,& use lattice, only: lattice_Sslip_v ,&
lattice_mu, & lattice_mu, &
lattice_nu lattice_nu
use mesh, only: mesh_NcpElems, & use mesh, only: mesh_ipVolume
mesh_maxNips, & use material, only: material_phase, &
mesh_ipVolume
use material, only: homogenization_maxNgrains, &
material_phase, &
plasticState, & plasticState, &
mappingConstitutive, & mappingConstitutive, &
phase_plasticityInstance phase_plasticityInstance
@ -2360,7 +2338,7 @@ deltaDUpper = dUpper - dUpperOld
!*** dissociation by stress increase !*** dissociation by stress increase
deltaRhoDipole2SingleStress = 0.0_pReal deltaRhoDipole2SingleStress = 0.0_pReal
forall (c=1_pInt:2_pInt, s=1_pInt:ns, deltaDUpper(s,c) < 0.0_pReal .and. & forall (c=1_pInt:2_pInt, s=1_pInt:ns, deltaDUpper(s,c) < 0.0_pReal .and. &
(dUpperOld(s,c) - dLower(s,c))/= 0.0_pReal) & (dUpperOld(s,c) - dLower(s,c)) > tiny(0.0_pReal)) &
deltaRhoDipole2SingleStress(s,8_pInt+c) = rhoDip(s,c) * deltaDUpper(s,c) & deltaRhoDipole2SingleStress(s,8_pInt+c) = rhoDip(s,c) * deltaDUpper(s,c) &
/ (dUpperOld(s,c) - dLower(s,c)) / (dUpperOld(s,c) - dLower(s,c))
@ -2856,11 +2834,11 @@ if (.not. phase_localPlasticity(material_phase(1_pInt,ip,el))) then
my_rhoSgl = rhoSgl my_rhoSgl = rhoSgl
my_v = v my_v = v
if(numerics_timeSyncing) then if(numerics_timeSyncing) then
if (subfrac(1_pInt,ip,el) == 0.0_pReal) then if (subfrac(1_pInt,ip,el) <= tiny(0.0_pReal)) then
my_rhoSgl = rhoSgl0 my_rhoSgl = rhoSgl0
my_v = v0 my_v = v0
elseif (neighbor_n > 0_pInt) then elseif (neighbor_n > 0_pInt) then
if (subfrac(1_pInt,neighbor_ip,neighbor_el) == 0.0_pReal) then if (subfrac(1_pInt,neighbor_ip,neighbor_el) <= tiny(0.0_pReal)) then
my_rhoSgl = rhoSgl0 my_rhoSgl = rhoSgl0
my_v = v0 my_v = v0
endif endif
@ -3206,7 +3184,7 @@ end subroutine plastic_nonlocal_updateCompatibility
function plastic_nonlocal_dislocationstress(Fe, ip, el) function plastic_nonlocal_dislocationstress(Fe, ip, el)
use math, only: math_mul33x33, & use math, only: math_mul33x33, &
math_mul33x3, & math_mul33x3, &
math_invert33, & math_inv33, &
math_transpose33, & math_transpose33, &
pi pi
use mesh, only: mesh_NcpElems, & use mesh, only: mesh_NcpElems, &
@ -3267,8 +3245,7 @@ real(pReal) x, y, z, &
R, Rsquare, Rcube, & R, Rsquare, Rcube, &
denominator, & denominator, &
flipSign, & flipSign, &
neighbor_ipVolumeSideLength, & neighbor_ipVolumeSideLength
detFe
real(pReal), dimension(3) :: connection, & !< connection vector between me and my neighbor in the deformed configuration real(pReal), dimension(3) :: connection, & !< connection vector between me and my neighbor in the deformed configuration
connection_neighborLattice, & !< connection vector between me and my neighbor in the lattice configuration of my neighbor connection_neighborLattice, & !< connection vector between me and my neighbor in the lattice configuration of my neighbor
connection_neighborSlip, & !< connection vector between me and my neighbor in the slip system frame of my neighbor connection_neighborSlip, & !< connection vector between me and my neighbor in the slip system frame of my neighbor
@ -3287,7 +3264,6 @@ real(pReal), dimension(2,maxval(totalNslip)) :: &
rhoExcessDead rhoExcessDead
real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1_pInt,ip,el))),8) :: & real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1_pInt,ip,el))),8) :: &
rhoSgl ! single dislocation density (edge+, edge-, screw+, screw-, used edge+, used edge-, used screw+, used screw-) rhoSgl ! single dislocation density (edge+, edge-, screw+, screw-, used edge+, used edge-, used screw+, used screw-)
logical inversionError
ph = material_phase(1_pInt,ip,el) ph = material_phase(1_pInt,ip,el)
instance = phase_plasticityInstance(ph) instance = phase_plasticityInstance(ph)
@ -3295,8 +3271,6 @@ ns = totalNslip(instance)
p = mappingConstitutive(2,1,ip,el) p = mappingConstitutive(2,1,ip,el)
o = mappingConstitutive(1,1,ip,el) o = mappingConstitutive(1,1,ip,el)
!*** 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)
@ -3312,7 +3286,7 @@ endforall
plastic_nonlocal_dislocationstress = 0.0_pReal plastic_nonlocal_dislocationstress = 0.0_pReal
if (.not. phase_localPlasticity(ph)) then if (.not. phase_localPlasticity(ph)) then
call math_invert33(Fe(1:3,1:3,1_pInt,ip,el), invFe, detFe, inversionError) invFe = math_inv33(Fe(1:3,1:3,1_pInt,ip,el))
!* in case of periodic surfaces we have to find out how many periodic images in each direction we need !* in case of periodic surfaces we have to find out how many periodic images in each direction we need
@ -3335,17 +3309,15 @@ if (.not. phase_localPlasticity(ph)) then
!* but only consider nonlocal neighbors within a certain cutoff radius R !* but only consider nonlocal neighbors within a certain cutoff radius R
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(1_pInt,neighbor_ip,neighbor_el) neighbor_phase = material_phase(1_pInt,neighbor_ip,neighbor_el)
np = mappingConstitutive(2,1,neighbor_ip,neighbor_el) np = mappingConstitutive(2,1,neighbor_ip,neighbor_el)
no = mappingConstitutive(1,1,neighbor_ip,neighbor_el) no = mappingConstitutive(1,1,neighbor_ip,neighbor_el)
if (phase_localPlasticity(neighbor_phase)) then if (phase_localPlasticity(neighbor_phase)) cycle
cycle
endif
neighbor_instance = phase_plasticityInstance(neighbor_phase) neighbor_instance = phase_plasticityInstance(neighbor_phase)
neighbor_ns = totalNslip(neighbor_instance) neighbor_ns = totalNslip(neighbor_instance)
call math_invert33(Fe(1:3,1:3,1,neighbor_ip,neighbor_el), neighbor_invFe, detFe, inversionError) neighbor_invFe = math_inv33(Fe(1:3,1:3,1,neighbor_ip,neighbor_el))
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)
@ -3370,9 +3342,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(instance)) cycle
cycle
endif
!* the segment length is the minimum of the third root of the control volume and the ip distance !* the segment length is the minimum of the third root of the control volume and the ip distance
@ -3424,7 +3394,7 @@ ipLoop: do neighbor_ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,neighbor_el))
Rsquare = R * R Rsquare = R * R
Rcube = Rsquare * R Rcube = Rsquare * R
denominator = R * (R + flipSign * lambda) denominator = R * (R + flipSign * lambda)
if (denominator == 0.0_pReal) exit ipLoop if (denominator <= tiny(0.0_pReal)) exit ipLoop
sigma(1,1) = sigma(1,1) - real(side,pReal) & sigma(1,1) = sigma(1,1) - real(side,pReal) &
* flipSign * z / denominator & * flipSign * z / denominator &
@ -3469,7 +3439,7 @@ ipLoop: do neighbor_ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,neighbor_el))
Rsquare = R * R Rsquare = R * R
Rcube = Rsquare * R Rcube = Rsquare * R
denominator = R * (R + flipSign * lambda) denominator = R * (R + flipSign * lambda)
if (denominator == 0.0_pReal) exit ipLoop if (denominator <= tiny(0.0_pReal)) exit ipLoop
sigma(1,2) = sigma(1,2) - real(side,pReal) * flipSign * z & sigma(1,2) = sigma(1,2) - real(side,pReal) * flipSign * z &
* (1.0_pReal - lattice_nu(ph)) / denominator & * (1.0_pReal - lattice_nu(ph)) / denominator &
@ -3571,8 +3541,7 @@ function plastic_nonlocal_postResults(Tstar_v,Fe,ip,el)
material_phase, & material_phase, &
mappingConstitutive, & mappingConstitutive, &
plasticState, & plasticState, &
phase_plasticityInstance, & phase_plasticityInstance
phase_Noutput
use lattice, only: & use lattice, only: &
lattice_Sslip_v, & lattice_Sslip_v, &
lattice_sd, & lattice_sd, &

View File

@ -209,7 +209,6 @@ subroutine plastic_titanmod_init(fileUnit)
IO_timeStamp, & IO_timeStamp, &
IO_EOF IO_EOF
use material, only: & use material, only: &
homogenization_maxNgrains, &
phase_plasticity, & phase_plasticity, &
phase_plasticityInstance, & phase_plasticityInstance, &
phase_Noutput, & phase_Noutput, &
@ -1074,8 +1073,7 @@ subroutine plastic_titanmod_stateInit(ph,instance)
lattice_mu lattice_mu
use material, only: & use material, only: &
plasticState, & plasticState
mappingConstitutive
implicit none implicit none
integer(pInt), intent(in) :: instance !< number specifying the instance of the plasticity integer(pInt), intent(in) :: instance !< number specifying the instance of the plasticity
@ -1167,11 +1165,7 @@ end subroutine plastic_titanmod_stateInit
!> @brief returns the homogenized elasticity matrix !> @brief returns the homogenized elasticity matrix
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function plastic_titanmod_homogenizedC(ipc,ip,el) function plastic_titanmod_homogenizedC(ipc,ip,el)
use mesh, only: &
mesh_NcpElems, &
mesh_maxNips
use material, only: & use material, only: &
homogenization_maxNgrains, &
material_phase, & material_phase, &
phase_plasticityInstance, & phase_plasticityInstance, &
plasticState, & plasticState, &
@ -1347,11 +1341,7 @@ subroutine plastic_titanmod_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,temperature,
lattice_NtwinSystem, & lattice_NtwinSystem, &
lattice_structure, & lattice_structure, &
LATTICE_hex_ID LATTICE_hex_ID
use mesh, only: &
mesh_NcpElems, &
mesh_maxNips
use material, only: & use material, only: &
homogenization_maxNgrains, &
material_phase, & material_phase, &
phase_plasticityInstance, & phase_plasticityInstance, &
plasticState, & plasticState, &
@ -1662,11 +1652,7 @@ subroutine plastic_titanmod_dotState(Tstar_v,temperature,ipc,ip,el)
lattice_maxNtwinFamily, & lattice_maxNtwinFamily, &
lattice_NslipSystem, & lattice_NslipSystem, &
lattice_NtwinSystem lattice_NtwinSystem
use mesh, only: &
mesh_NcpElems, &
mesh_maxNips
use material, only: & use material, only: &
homogenization_maxNgrains, &
material_phase, & material_phase, &
phase_plasticityInstance, & phase_plasticityInstance, &
plasticState, & plasticState, &
@ -1782,14 +1768,9 @@ end subroutine plastic_titanmod_dotState
!> @brief return array of constitutive results !> @brief return array of constitutive results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function plastic_titanmod_postResults(ipc,ip,el) function plastic_titanmod_postResults(ipc,ip,el)
use mesh, only: &
mesh_NcpElems, &
mesh_maxNips
use material, only: & use material, only: &
homogenization_maxNgrains, &
material_phase, & material_phase, &
phase_plasticityInstance, & phase_plasticityInstance, &
phase_Noutput, &
plasticState, & plasticState, &
mappingConstitutive mappingConstitutive