diff --git a/code/constitutive.f90 b/code/constitutive.f90 index cf97abf6c..abe6116cb 100644 --- a/code/constitutive.f90 +++ b/code/constitutive.f90 @@ -33,18 +33,16 @@ module constitutive constitutive_sizeDotState, & !< size of dotState array constitutive_sizeState, & !< size of state array per grain constitutive_sizePostResults !< size of postResults array per grain - integer(pInt), public, protected :: & - constitutive_maxSizeDotState, & - constitutive_maxSizePostResults integer(pInt), private :: & constitutive_maxSizeState #else integer(pInt), public, dimension(:,:,:), allocatable :: & constitutive_sizePostResults !< size of postResults array per grain +#endif integer(pInt), public, protected :: & constitutive_maxSizePostResults, & constitutive_maxSizeDotState -#endif + public :: & constitutive_init, & constitutive_homogenizedC, & @@ -445,7 +443,7 @@ subroutine constitutive_init if (nonlocalConstitutionPresent) & #ifdef NEWSTATE - call constitutive_nonlocal_stateInit(mappingConstitutive) + call constitutive_nonlocal_stateInit() #else call constitutive_nonlocal_stateInit(constitutive_state0(1,1:iMax,1:eMax)) #endif @@ -633,7 +631,7 @@ subroutine constitutive_microstructure(temperature, Fe, Fp, ipc, ip, el) #endif case (PLASTICITY_NONLOCAL_ID) #ifdef NEWSTATE - call constitutive_nonlocal_microstructure(mappingConstitutive,Fe,Fp,ipc,ip,el) + call constitutive_nonlocal_microstructure(Fe,Fp,ipc,ip,el) #else call constitutive_nonlocal_microstructure(constitutive_state,Fe,Fp,ipc,ip,el) #endif @@ -730,7 +728,7 @@ subroutine constitutive_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, temperature, ip case (PLASTICITY_NONLOCAL_ID) #ifdef NEWSTATE call constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, & - temperature, mappingConstitutive, ipc,ip,el) + temperature, ipc,ip,el) #else call constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar, Tstar_v, & temperature, constitutive_state(ipc,ip,el), ipc,ip,el) @@ -952,9 +950,10 @@ subroutine constitutive_collectDotState(Tstar_v, FeArray, FpArray, Temperature, #endif case (PLASTICITY_NONLOCAL_ID) #ifdef NEWSTATE +!* plasticState(mappingConstitutive(2,ipc,ip,el))%dotState(:,mappingConstitutive(1,ipc,ip,el)) = & constitutive_nonlocal_dotState(Tstar_v, FeArray, FpArray, & - Temperature, mappingConstitutive, subdt, & + Temperature, subdt, & subfracArray, ipc, ip, el) #else @@ -1019,7 +1018,7 @@ logical function constitutive_collectDeltaState(Tstar_v, ipc, ip, el) case (PLASTICITY_NONLOCAL_ID) constitutive_collectDeltaState = .true. #ifdef NEWSTATE - call constitutive_nonlocal_deltaState(mappingConstitutive, Tstar_v,ipc,ip,el) + call constitutive_nonlocal_deltaState(Tstar_v,ip,el) #else call constitutive_nonlocal_deltaState(constitutive_deltaState(ipc,ip,el),& constitutive_state(ipc,ip,el), Tstar_v,ipc,ip,el) diff --git a/code/constitutive_dislotwin.f90 b/code/constitutive_dislotwin.f90 index 87a9d1c90..3e49accb2 100644 --- a/code/constitutive_dislotwin.f90 +++ b/code/constitutive_dislotwin.f90 @@ -338,93 +338,93 @@ subroutine constitutive_dislotwin_init(fileUnit) case ('(output)') select case(IO_lc(IO_stringValue(line,positions,2_pInt))) case ('edge_density') - constitutive_dislotwin_outputID(constitutive_dislotwin_Noutput(instance),instance) = edge_density_ID constitutive_dislotwin_Noutput(instance) = constitutive_dislotwin_Noutput(instance) + 1_pInt + constitutive_dislotwin_outputID(constitutive_dislotwin_Noutput(instance),instance) = edge_density_ID constitutive_dislotwin_output(constitutive_dislotwin_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('dipole_density') - constitutive_dislotwin_outputID(constitutive_dislotwin_Noutput(instance),instance) = dipole_density_ID constitutive_dislotwin_Noutput(instance) = constitutive_dislotwin_Noutput(instance) + 1_pInt + constitutive_dislotwin_outputID(constitutive_dislotwin_Noutput(instance),instance) = dipole_density_ID constitutive_dislotwin_output(constitutive_dislotwin_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('shear_rate_slip') - constitutive_dislotwin_outputID(constitutive_dislotwin_Noutput(instance),instance) = shear_rate_slip_ID constitutive_dislotwin_Noutput(instance) = constitutive_dislotwin_Noutput(instance) + 1_pInt + constitutive_dislotwin_outputID(constitutive_dislotwin_Noutput(instance),instance) = shear_rate_slip_ID constitutive_dislotwin_output(constitutive_dislotwin_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('accumulated_shear_slip') - constitutive_dislotwin_outputID(constitutive_dislotwin_Noutput(instance),instance) = accumulated_shear_slip_ID constitutive_dislotwin_Noutput(instance) = constitutive_dislotwin_Noutput(instance) + 1_pInt + constitutive_dislotwin_outputID(constitutive_dislotwin_Noutput(instance),instance) = accumulated_shear_slip_ID constitutive_dislotwin_output(constitutive_dislotwin_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('mfp_slip') - constitutive_dislotwin_outputID(constitutive_dislotwin_Noutput(instance),instance) = mfp_slip_ID constitutive_dislotwin_Noutput(instance) = constitutive_dislotwin_Noutput(instance) + 1_pInt + constitutive_dislotwin_outputID(constitutive_dislotwin_Noutput(instance),instance) = mfp_slip_ID constitutive_dislotwin_output(constitutive_dislotwin_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('resolved_stress_slip') - constitutive_dislotwin_outputID(constitutive_dislotwin_Noutput(instance),instance) = resolved_stress_slip_ID constitutive_dislotwin_Noutput(instance) = constitutive_dislotwin_Noutput(instance) + 1_pInt + constitutive_dislotwin_outputID(constitutive_dislotwin_Noutput(instance),instance) = resolved_stress_slip_ID constitutive_dislotwin_output(constitutive_dislotwin_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('edge_dipole_distance') - constitutive_dislotwin_outputID(constitutive_dislotwin_Noutput(instance),instance) = edge_dipole_distance_ID constitutive_dislotwin_Noutput(instance) = constitutive_dislotwin_Noutput(instance) + 1_pInt + constitutive_dislotwin_outputID(constitutive_dislotwin_Noutput(instance),instance) = edge_dipole_distance_ID constitutive_dislotwin_output(constitutive_dislotwin_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('stress_exponent') - constitutive_dislotwin_outputID(constitutive_dislotwin_Noutput(instance),instance) = stress_exponent_ID constitutive_dislotwin_Noutput(instance) = constitutive_dislotwin_Noutput(instance) + 1_pInt + constitutive_dislotwin_outputID(constitutive_dislotwin_Noutput(instance),instance) = stress_exponent_ID constitutive_dislotwin_output(constitutive_dislotwin_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('twin_fraction') - constitutive_dislotwin_outputID(constitutive_dislotwin_Noutput(instance),instance) = twin_fraction_ID constitutive_dislotwin_Noutput(instance) = constitutive_dislotwin_Noutput(instance) + 1_pInt + constitutive_dislotwin_outputID(constitutive_dislotwin_Noutput(instance),instance) = twin_fraction_ID constitutive_dislotwin_output(constitutive_dislotwin_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('shear_rate_twin') - constitutive_dislotwin_outputID(constitutive_dislotwin_Noutput(instance),instance) = shear_rate_twin_ID constitutive_dislotwin_Noutput(instance) = constitutive_dislotwin_Noutput(instance) + 1_pInt + constitutive_dislotwin_outputID(constitutive_dislotwin_Noutput(instance),instance) = shear_rate_twin_ID constitutive_dislotwin_output(constitutive_dislotwin_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('accumulated_shear_twin') - constitutive_dislotwin_outputID(constitutive_dislotwin_Noutput(instance),instance) = accumulated_shear_twin_ID constitutive_dislotwin_Noutput(instance) = constitutive_dislotwin_Noutput(instance) + 1_pInt + constitutive_dislotwin_outputID(constitutive_dislotwin_Noutput(instance),instance) = accumulated_shear_twin_ID constitutive_dislotwin_output(constitutive_dislotwin_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('mfp_twin') - constitutive_dislotwin_outputID(constitutive_dislotwin_Noutput(instance),instance) = mfp_twin_ID constitutive_dislotwin_Noutput(instance) = constitutive_dislotwin_Noutput(instance) + 1_pInt + constitutive_dislotwin_outputID(constitutive_dislotwin_Noutput(instance),instance) = mfp_twin_ID constitutive_dislotwin_output(constitutive_dislotwin_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('resolved_stress_twin') - constitutive_dislotwin_outputID(constitutive_dislotwin_Noutput(instance),instance) = resolved_stress_twin_ID constitutive_dislotwin_Noutput(instance) = constitutive_dislotwin_Noutput(instance) + 1_pInt + constitutive_dislotwin_outputID(constitutive_dislotwin_Noutput(instance),instance) = resolved_stress_twin_ID constitutive_dislotwin_output(constitutive_dislotwin_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('threshold_stress_twin') - constitutive_dislotwin_outputID(constitutive_dislotwin_Noutput(instance),instance) = threshold_stress_twin_ID constitutive_dislotwin_Noutput(instance) = constitutive_dislotwin_Noutput(instance) + 1_pInt + constitutive_dislotwin_outputID(constitutive_dislotwin_Noutput(instance),instance) = threshold_stress_twin_ID constitutive_dislotwin_output(constitutive_dislotwin_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('resolved_stress_shearband') - constitutive_dislotwin_outputID(constitutive_dislotwin_Noutput(instance),instance) = resolved_stress_shearband_ID constitutive_dislotwin_Noutput(instance) = constitutive_dislotwin_Noutput(instance) + 1_pInt + constitutive_dislotwin_outputID(constitutive_dislotwin_Noutput(instance),instance) = resolved_stress_shearband_ID constitutive_dislotwin_output(constitutive_dislotwin_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('shear_rate_shearband') - constitutive_dislotwin_outputID(constitutive_dislotwin_Noutput(instance),instance) = shear_rate_shearband_ID constitutive_dislotwin_Noutput(instance) = constitutive_dislotwin_Noutput(instance) + 1_pInt + constitutive_dislotwin_outputID(constitutive_dislotwin_Noutput(instance),instance) = shear_rate_shearband_ID constitutive_dislotwin_output(constitutive_dislotwin_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('sb_eigenvalues') - constitutive_dislotwin_outputID(constitutive_dislotwin_Noutput(instance),instance) = sb_eigenvalues_ID constitutive_dislotwin_Noutput(instance) = constitutive_dislotwin_Noutput(instance) + 1_pInt + constitutive_dislotwin_outputID(constitutive_dislotwin_Noutput(instance),instance) = sb_eigenvalues_ID constitutive_dislotwin_output(constitutive_dislotwin_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('sb_eigenvectors') - constitutive_dislotwin_outputID(constitutive_dislotwin_Noutput(instance),instance) = sb_eigenvectors_ID constitutive_dislotwin_Noutput(instance) = constitutive_dislotwin_Noutput(instance) + 1_pInt + constitutive_dislotwin_outputID(constitutive_dislotwin_Noutput(instance),instance) = sb_eigenvectors_ID constitutive_dislotwin_output(constitutive_dislotwin_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) end select diff --git a/code/constitutive_j2.f90 b/code/constitutive_j2.f90 index 9d8b72e5c..75be972fe 100644 --- a/code/constitutive_j2.f90 +++ b/code/constitutive_j2.f90 @@ -225,8 +225,8 @@ subroutine constitutive_j2_init(fileUnit) case ('(output)') select case(IO_lc(IO_stringValue(line,positions,2_pInt))) case ('flowstress') - constitutive_j2_outputID(constitutive_j2_Noutput(instance),instance) = flowstress_ID constitutive_j2_Noutput(instance) = constitutive_j2_Noutput(instance) + 1_pInt + constitutive_j2_outputID(constitutive_j2_Noutput(instance),instance) = flowstress_ID constitutive_j2_output(constitutive_j2_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) #ifdef HDF @@ -235,8 +235,8 @@ subroutine constitutive_j2_init(fileUnit) constitutive_j2_Output2(instance)%flowstressActive = .true. #endif case ('strainrate') - constitutive_j2_outputID(constitutive_j2_Noutput(instance),instance) = strainrate_ID constitutive_j2_Noutput(instance) = constitutive_j2_Noutput(instance) + 1_pInt + constitutive_j2_outputID(constitutive_j2_Noutput(instance),instance) = strainrate_ID constitutive_j2_output(constitutive_j2_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) #ifdef HDF diff --git a/code/constitutive_nonlocal.f90 b/code/constitutive_nonlocal.f90 index 2b01c07fb..64a7cb698 100644 --- a/code/constitutive_nonlocal.f90 +++ b/code/constitutive_nonlocal.f90 @@ -437,340 +437,340 @@ allocate(nonSchmidCoeff(lattice_maxNnonSchmid,maxNinstances), s case ('(output)') select case(IO_lc(IO_stringValue(line,positions,2_pInt))) case('rho') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('delta') - constitutive_nonlocal_outputID(Noutput(instance),instance) = delta_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = delta_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_edge') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_edge_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_edge_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_screw') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_screw_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_screw_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_sgl') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('delta_sgl') - constitutive_nonlocal_outputID(Noutput(instance),instance) = delta_sgl_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = delta_sgl_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_sgl_edge') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_edge_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_edge_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_sgl_edge_pos') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_edge_pos_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_edge_pos_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_sgl_edge_neg') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_edge_neg_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_edge_neg_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_sgl_screw') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_screw_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_screw_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_sgl_screw_pos') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_screw_pos_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_screw_pos_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_sgl_screw_neg') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_screw_neg_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_screw_neg_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_sgl_mobile') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_mobile_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_mobile_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_sgl_edge_mobile') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_edge_mobile_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_edge_mobile_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_sgl_edge_pos_mobile') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_edge_pos_mobile_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_edge_pos_mobile_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_sgl_edge_neg_mobile') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_edge_neg_mobile_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_edge_neg_mobile_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_sgl_screw_mobile') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_screw_mobile_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_screw_mobile_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_sgl_screw_pos_mobile') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_screw_pos_mobile_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_screw_pos_mobile_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_sgl_screw_neg_mobile') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_screw_neg_mobile_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_screw_neg_mobile_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_sgl_immobile') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_immobile_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_immobile_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_sgl_edge_immobile') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_edge_immobile_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_edge_immobile_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_sgl_edge_pos_immobile') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_edge_pos_immobile_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_edge_pos_immobile_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_sgl_edge_neg_immobile') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_edge_neg_immobile_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_edge_neg_immobile_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_sgl_screw_immobile') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_screw_immobile_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_screw_immobile_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_sgl_screw_pos_immobile') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_screw_pos_immobile_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_screw_pos_immobile_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_sgl_screw_neg_immobile') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_screw_neg_immobile_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_sgl_screw_neg_immobile_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_dip') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dip_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dip_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('delta_dip') - constitutive_nonlocal_outputID(Noutput(instance),instance) = delta_dip_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = delta_dip_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_dip_edge') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dip_edge_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dip_edge_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_dip_screw') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dip_screw_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dip_screw_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('excess_rho') - constitutive_nonlocal_outputID(Noutput(instance),instance) = excess_rho_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = excess_rho_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('excess_rho_edge') - constitutive_nonlocal_outputID(Noutput(instance),instance) = excess_rho_edge_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = excess_rho_edge_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('excess_rho_screw') - constitutive_nonlocal_outputID(Noutput(instance),instance) = excess_rho_screw_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = excess_rho_screw_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_forest') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_forest_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_forest_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('shearrate') - constitutive_nonlocal_outputID(Noutput(instance),instance) = shearrate_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = shearrate_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('resolvedstress') - constitutive_nonlocal_outputID(Noutput(instance),instance) = resolvedstress_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = resolvedstress_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('resolvedstress_external') - constitutive_nonlocal_outputID(Noutput(instance),instance) = resolvedstress_external_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = resolvedstress_external_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('resolvedstress_back') - constitutive_nonlocal_outputID(Noutput(instance),instance) = resolvedstress_back_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = resolvedstress_back_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('resistance') - constitutive_nonlocal_outputID(Noutput(instance),instance) = resistance_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = resistance_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_dot') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dot_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dot_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_dot_sgl') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dot_sgl_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dot_sgl_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_dot_sgl_mobile') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dot_sgl_mobile_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dot_sgl_mobile_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_dot_dip') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dot_dip_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dot_dip_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_dot_gen') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dot_gen_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dot_gen_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_dot_gen_edge') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dot_gen_edge_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dot_gen_edge_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_dot_gen_screw') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dot_gen_screw_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dot_gen_screw_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_dot_sgl2dip') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dot_sgl2dip_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dot_sgl2dip_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_dot_sgl2dip_edge') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dot_sgl2dip_edge_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dot_sgl2dip_edge_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_dot_sgl2dip_screw') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dot_sgl2dip_screw_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dot_sgl2dip_screw_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_dot_ann_ath') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dot_ann_ath_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dot_ann_ath_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_dot_ann_the') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dot_ann_the_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dot_ann_the_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_dot_ann_the_edge') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dot_ann_the_edge_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dot_ann_the_edge_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_dot_ann_the_screw') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dot_ann_the_screw_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dot_ann_the_screw_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_dot_edgejogs') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dot_edgejogs_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dot_edgejogs_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_dot_flux') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dot_flux_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dot_flux_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_dot_flux_mobile') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dot_flux_mobile_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dot_flux_mobile_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_dot_flux_edge') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dot_flux_edge_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dot_flux_edge_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('rho_dot_flux_screw') - constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dot_flux_screw_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = rho_dot_flux_screw_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('velocity_edge_pos') - constitutive_nonlocal_outputID(Noutput(instance),instance) = velocity_edge_pos_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = velocity_edge_pos_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('velocity_edge_neg') - constitutive_nonlocal_outputID(Noutput(instance),instance) = velocity_edge_neg_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = velocity_edge_neg_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('velocity_screw_pos') - constitutive_nonlocal_outputID(Noutput(instance),instance) = velocity_screw_pos_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = velocity_screw_pos_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('velocity_screw_neg') - constitutive_nonlocal_outputID(Noutput(instance),instance) = velocity_screw_neg_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = velocity_screw_neg_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('slipdirection.x') - constitutive_nonlocal_outputID(Noutput(instance),instance) = slipdirectionx_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = slipdirectionx_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('slipdirection.y') - constitutive_nonlocal_outputID(Noutput(instance),instance) = slipdirectiony_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = slipdirectiony_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('slipdirection.z') - constitutive_nonlocal_outputID(Noutput(instance),instance) = slipdirectionz_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = slipdirectionz_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('slipnormal.x') - constitutive_nonlocal_outputID(Noutput(instance),instance) = slipnormalx_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = slipnormalx_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('slipnormal.y') - constitutive_nonlocal_outputID(Noutput(instance),instance) = slipnormaly_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = slipnormaly_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('slipnormal.z') - constitutive_nonlocal_outputID(Noutput(instance),instance) = slipnormalz_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = slipnormalz_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('fluxdensity_edge_pos.x') - constitutive_nonlocal_outputID(Noutput(instance),instance) = fluxdensity_edge_posx_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = fluxdensity_edge_posx_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('fluxdensity_edge_pos.y') - constitutive_nonlocal_outputID(Noutput(instance),instance) = fluxdensity_edge_posy_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = fluxdensity_edge_posy_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('fluxdensity_edge_pos.z') - constitutive_nonlocal_outputID(Noutput(instance),instance) = fluxdensity_edge_posz_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = fluxdensity_edge_posz_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('fluxdensity_edge_neg.x') - constitutive_nonlocal_outputID(Noutput(instance),instance) = fluxdensity_edge_negx_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = fluxdensity_edge_negx_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('fluxdensity_edge_neg.y') - constitutive_nonlocal_outputID(Noutput(instance),instance) = fluxdensity_edge_negy_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = fluxdensity_edge_negy_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('fluxdensity_edge_neg.z') - constitutive_nonlocal_outputID(Noutput(instance),instance) = fluxdensity_edge_negz_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = fluxdensity_edge_negz_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('fluxdensity_screw_pos.x') - constitutive_nonlocal_outputID(Noutput(instance),instance) = fluxdensity_screw_posx_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = fluxdensity_screw_posx_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('fluxdensity_screw_pos.y') - constitutive_nonlocal_outputID(Noutput(instance),instance) = fluxdensity_screw_posy_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = fluxdensity_screw_posy_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('fluxdensity_screw_pos.z') - constitutive_nonlocal_outputID(Noutput(instance),instance) = fluxdensity_screw_posz_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = fluxdensity_screw_posz_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('fluxdensity_screw_neg.x') - constitutive_nonlocal_outputID(Noutput(instance),instance) = fluxdensity_screw_negx_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = fluxdensity_screw_negx_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('fluxdensity_screw_neg.y') - constitutive_nonlocal_outputID(Noutput(instance),instance) = fluxdensity_screw_negy_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = fluxdensity_screw_negy_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('fluxdensity_screw_neg.z') - constitutive_nonlocal_outputID(Noutput(instance),instance) = fluxdensity_screw_negz_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = fluxdensity_screw_negz_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('maximumdipoleheight_edge') - constitutive_nonlocal_outputID(Noutput(instance),instance) = maximumdipoleheight_edge_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = maximumdipoleheight_edge_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('maximumdipoleheight_screw') - constitutive_nonlocal_outputID(Noutput(instance),instance) = maximumdipoleheight_screw_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = maximumdipoleheight_screw_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('accumulatedshear') - constitutive_nonlocal_outputID(Noutput(instance),instance) = accumulatedshear_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = accumulatedshear_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) case('dislocationstress') - constitutive_nonlocal_outputID(Noutput(instance),instance) = dislocationstress_ID Noutput(instance) = Noutput(instance) + 1_pInt + constitutive_nonlocal_outputID(Noutput(instance),instance) = dislocationstress_ID constitutive_nonlocal_output(Noutput(instance),instance) = IO_lc(IO_stringValue(line,positions,2_pInt)) end select case ('nslip') @@ -1068,11 +1068,10 @@ allocate(nonSchmidProjection(3,3,4,maxTotalNslip,maxNinstances), #ifdef NEWSTATE ! Determine size of state array ! plasticState(phase)%nonlocal = .true. - sizeDotState = int(size(BASICSTATES),pInt) * ns - sizeDependentState = int(size(DEPENDENTSTATES),pInt) * ns - sizeState = sizeDotState & - + constitutive_nonlocal_sizeDependentState(instance) & - + int(size(OTHERSTATES),pInt) * ns + sizeDotState = int(size(BASICSTATES),pInt) * ns + sizeDependentState = int(size(DEPENDENTSTATES),pInt) * ns + sizeState = sizeDotState + sizeDependentState & + + int(size(OTHERSTATES),pInt) * ns plasticState(phase)%sizeState = sizeState plasticState(phase)%sizeDotState = sizeDotState @@ -1323,34 +1322,45 @@ allocate(nonSchmidProjection(3,3,4,maxTotalNslip,maxNinstances), enddo initializeInstances end subroutine constitutive_nonlocal_init -#ifdef NEWSTATE + !-------------------------------------------------------------------------------------------------- !> @brief sets the initial microstructural state for a given instance of this plasticity !-------------------------------------------------------------------------------------------------- -subroutine constitutive_nonlocal_stateInit(mapping) +#ifdef NEWSTATE +subroutine constitutive_nonlocal_stateInit() +#else +subroutine constitutive_nonlocal_stateInit(state) +#endif use IO, only: IO_error use lattice, only: lattice_maxNslipFamily use math, only: math_sampleGaussVar use mesh, only: mesh_ipVolume, & mesh_NcpElems, & - mesh_maxNips, & - mesh_element, & - FE_Nips, & + mesh_maxNips, & + mesh_element, & + FE_Nips, & FE_geomtype use material, only: material_phase, & phase_plasticityInstance, & +#ifdef NEWSTATE plasticState, & + mappingConstitutive, & material_Nphase, & - homogenization_maxNgrains, & +#endif phase_plasticity ,& PLASTICITY_NONLOCAL_ID -use numerics,only: & +#ifdef NEWSTATE + use numerics,only: & numerics_integrator +#endif implicit none -integer(pInt),dimension(2,homogenization_maxngrains,mesh_maxNips,mesh_ncpelems), intent(in) :: & - mapping +#ifndef NEWSTATE +type(p_vec), dimension(1,mesh_maxNips,mesh_NcpElems), intent(inout) :: & + state ! microstructural state +#endif + integer(pInt) el, & ip, & e, & @@ -1363,6 +1373,7 @@ integer(pInt) el, & t, & j, & instance, & + phase, & maxNinstances real(pReal), dimension(2) :: noise real(pReal), dimension(4) :: rnd @@ -1373,6 +1384,16 @@ real(pReal) meanDensity, & maxNinstances = int(count(phase_plasticity == PLASTICITY_NONLOCAL_ID),pInt) +#ifndef NEWSTATE +! ititalize all states to zero (already done for new state) +do e = 1_pInt,mesh_NcpElems + do i = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,e))) + if (PLASTICITY_NONLOCAL_ID == phase_plasticity(material_phase(1,i,e))) & + state(1,i,e)%p = 0.0_pReal + enddo +enddo +#endif + do instance = 1_pInt,maxNinstances ns = totalNslip(instance) @@ -1406,18 +1427,13 @@ do instance = 1_pInt,maxNinstances s = nint(rnd(3)*real(ns,pReal)+0.5_pReal,pInt) t = nint(rnd(4)*4.0_pReal+0.5_pReal,pInt) meanDensity = meanDensity + densityBinning * mesh_ipVolume(ip,el) / totalVolume - plasticState(mapping(2,1,ip,el))%state0(iRhoU(s,t,instance),mapping(1,1,ip,el)) = & - plasticState(mapping(2,1,ip,el))%state0(iRhoU(s,t,instance),mapping(1,1,ip,el)) & - + densityBinning - - plasticState(mapping(2,1,ip,el))%state(iRhoU(s,t,instance),mapping(1,1,ip,el)) = & - plasticState(mapping(2,1,ip,el))%state(iRhoU(s,t,instance),mapping(1,1,ip,el)) & - + densityBinning - - plasticState(mapping(2,1,ip,el))%partionedState0(iRhoU(s,t,instance),mapping(1,1,ip,el)) = & - plasticState(mapping(2,1,ip,el))%partionedState0(iRhoU(s,t,instance),mapping(1,1,ip,el)) & - + densityBinning - +#ifdef NEWSTATE + plasticState(mappingConstitutive(2,1,ip,el))%state0(iRhoU(s,t,instance),mappingConstitutive(2,1,ip,el)) = & + plasticState(mappingConstitutive(2,1,ip,el))%state0(iRhoU(s,t,instance),mappingConstitutive(2,1,ip,el)) & + + densityBinning +#else + state(1,ip,el)%p(iRhoU(s,t,instance)) = state(1,ip,el)%p(iRhoU(s,t,instance)) + densityBinning +#endif endif enddo ! homogeneous distribution of density with some noise @@ -1433,75 +1449,42 @@ do instance = 1_pInt,maxNinstances do j = 1_pInt,2_pInt noise(j) = math_sampleGaussVar(0.0_pReal, rhoSglScatter(instance)) enddo +#ifdef NEWSTATE + plasticState(mappingConstitutive(2,1,i,e))%state0(iRhoU(s,1,instance),mappingConstitutive(2,1,i,e)) = & + rhoSglEdgePos0(f,instance) + noise(1) + plasticState(mappingConstitutive(2,1,i,e))%state0(iRhoU(s,2,instance),mappingConstitutive(2,1,i,e)) = & + rhoSglEdgeNeg0(f,instance) + noise(1) + plasticState(mappingConstitutive(2,1,i,e))%state0(iRhoU(s,3,instance),mappingConstitutive(2,1,i,e)) = & + rhoSglScrewPos0(f,instance) + noise(2) + plasticState(mappingConstitutive(2,1,i,e))%state0(iRhoU(s,4,instance),mappingConstitutive(2,1,i,e)) = & + rhoSglScrewNeg0(f,instance) + noise(2) +#else - plasticState(mapping(2,1,i,e))%state0(iRhoU(s,1,instance),mapping(1,1,i,e)) = & - rhoSglEdgePos0(f,instance) + noise(1) - plasticState(mapping(2,1,i,e))%state0(iRhoU(s,2,instance),mapping(1,1,i,e)) = & - rhoSglEdgeNeg0(f,instance) + noise(1) - plasticState(mapping(2,1,i,e))%state0(iRhoU(s,3,instance),mapping(1,1,i,e)) = & - rhoSglScrewPos0(f,instance) + noise(2) - plasticState(mapping(2,1,i,e))%state0(iRhoU(s,4,instance),mapping(1,1,i,e)) = & - rhoSglScrewNeg0(f,instance) + noise(2) - -! plasticState(mapping(2,1,i,e))%state(iRhoU(s,1,instance),mapping(1,1,i,e)) = & -! rhoSglEdgePos0(f,instance) + noise(1) -! plasticState(mapping(2,1,i,e))%state(iRhoU(s,2,instance),mapping(1,1,i,e)) = & -! rhoSglEdgeNeg0(f,instance) + noise(1) -! plasticState(mapping(2,1,i,e))%state(iRhoU(s,3,instance),mapping(1,1,i,e)) = & -! rhoSglScrewPos0(f,instance) + noise(2) -! plasticState(mapping(2,1,i,e))%state(iRhoU(s,4,instance),mapping(1,1,i,e)) = & -! rhoSglScrewNeg0(f,instance) + noise(2) -! -! plasticState(mapping(2,1,i,e))%partionedState0(iRhoU(s,1,instance),mapping(1,1,i,e)) = & -! rhoSglEdgePos0(f,instance) + noise(1) -! plasticState(mapping(2,1,i,e))%partionedState0(iRhoU(s,2,instance),mapping(1,1,i,e)) = & -! rhoSglEdgeNeg0(f,instance) + noise(1) -! plasticState(mapping(2,1,i,e))%partionedState0(iRhoU(s,3,instance),mapping(1,1,i,e)) = & -! rhoSglScrewPos0(f,instance) + noise(2) -! plasticState(mapping(2,1,i,e))%partionedState0(iRhoU(s,4,instance),mapping(1,1,i,e)) = & -! rhoSglScrewNeg0(f,instance) + noise(2) - - - enddo - - plasticState(mapping(2,1,i,e))%state0(iRhoD(from:upto,1,instance),mapping(1,1,i,e)) = & - rhoDipEdge0(f,instance) - plasticState(mapping(2,1,i,e))%state0(iRhoD(from:upto,2,instance),mapping(1,1,i,e)) = & - rhoDipScrew0(f,instance) - -! plasticState(mapping(2,1,i,e))%state(iRhoD(from:upto,1,instance),mapping(1,1,i,e)) = & -! rhoDipEdge0(f,instance) -! plasticState(mapping(2,1,i,e))%state(iRhoD(from:upto,2,instance),mapping(1,1,i,e)) = & -! rhoDipScrew0(f,instance) -! -! plasticState(mapping(2,1,i,e))%partionedState0(iRhoD(from:upto,1,instance),mapping(1,1,i,e)) = & -! rhoDipEdge0(f,instance) -! plasticState(mapping(2,1,i,e))%partionedState0(iRhoD(from:upto,2,instance),mapping(1,1,i,e)) = & -! rhoDipScrew0(f,instance) - - + state(1,i,e)%p(iRhoU(s,1,instance)) = rhoSglEdgePos0(f,instance) + noise(1) + state(1,i,e)%p(iRhoU(s,2,instance)) = rhoSglEdgeNeg0(f,instance) + noise(1) + state(1,i,e)%p(iRhoU(s,3,instance)) = rhoSglScrewPos0(f,instance) + noise(2) + state(1,i,e)%p(iRhoU(s,4,instance)) = rhoSglScrewNeg0(f,instance) + noise(2) +#endif + enddo +#ifdef NEWSTATE + plasticState(mappingConstitutive(2,1,i,e))%state0(iRhoD(from:upto,1,instance),mappingConstitutive(2,1,i,e)) = & + rhoDipEdge0(f,instance) + plasticState(mappingConstitutive(2,1,i,e))%state0(iRhoD(from:upto,2,instance),mappingConstitutive(2,1,i,e)) = & + rhoDipScrew0(f,instance) +#else + state(1,i,e)%p(iRhoD(from:upto,1,instance)) = rhoDipEdge0(f,instance) + state(1,i,e)%p(iRhoD(from:upto,2,instance)) = rhoDipScrew0(f,instance) +#endif enddo endif enddo enddo -!do e = 1_pInt,mesh_NcpElems -! do i = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,e))) -! if(plasticState(mapping(2,1,i,e))%nonlocal == .true.) then -! if (PLASTICITY_NONLOCAL_ID == phase_plasticity(material_phase(1,i,e)) & -! .and. instance == phase_plasticityInstance(material_phase(1,i,e))) then -! plasticState(mapping(2,1,i,e))%partionedState0(:,mapping(1,1,i,e)) = & -! plasticState(mapping(2,1,i,e))%state0(:,mapping(1,1,i,e)) -! plasticState(mapping(2,1,i,e))%State(:,mapping(1,1,i,e)) = & -! plasticState(mapping(2,1,i,e))%state0(:,mapping(1,1,i,e)) -! endif -! enddo -!enddo - endif enddo end subroutine constitutive_nonlocal_stateInit +#ifdef NEWSTATE !-------------------------------------------------------------------------------------------------- !> @brief sets the relevant state values for a given instance of this plasticity !-------------------------------------------------------------------------------------------------- @@ -1532,140 +1515,7 @@ forall (c = 1_pInt:2_pInt) & tempTol(1:plasticState(phase)%sizeDotState) end subroutine constitutive_nonlocal_aTolState - #else -!-------------------------------------------------------------------------------------------------- -!> @brief sets the initial microstructural state for a given instance of this plasticity -!-------------------------------------------------------------------------------------------------- -subroutine constitutive_nonlocal_stateInit(state) -use IO, only: IO_error -use lattice, only: lattice_maxNslipFamily -use math, only: math_sampleGaussVar -use mesh, only: mesh_ipVolume, & - mesh_NcpElems, & - mesh_maxNips, & - mesh_element, & - FE_Nips, & - FE_geomtype -use material, only: material_phase, & - phase_plasticityInstance, & -#ifdef NEWSTATE - plasticState, & - material_Nphase, & -#endif - phase_plasticity ,& - PLASTICITY_NONLOCAL_ID -#ifdef NEWSTATE - use numerics,only: & - numerics_integrator -#endif - -implicit none - -!*** input/output variables -type(p_vec), dimension(1,mesh_maxNips,mesh_NcpElems), intent(inout) :: & - state ! microstructural state - -!*** local variables -integer(pInt) el, & - ip, & - e, & - i, & - ns, & ! short notation for total number of active slip systems - f, & ! index of lattice family - from, & - upto, & - s, & ! index of slip system - t, & - j, & - instance, & - phase, & - maxNinstances -real(pReal), dimension(2) :: noise -real(pReal), dimension(4) :: rnd -real(pReal) meanDensity, & - totalVolume, & - densityBinning, & - minimumIpVolume - -maxNinstances = int(count(phase_plasticity == PLASTICITY_NONLOCAL_ID),pInt) - -! ititalize all states to zero - -do e = 1_pInt,mesh_NcpElems - do i = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,e))) - if (PLASTICITY_NONLOCAL_ID == phase_plasticity(material_phase(1,i,e))) & - state(1,i,e)%p = 0.0_pReal - enddo -enddo - - -do instance = 1_pInt,maxNinstances - ns = totalNslip(instance) - - ! randomly distribute dislocation segments on random slip system and of random type in the volume - if (rhoSglRandom(instance) > 0.0_pReal) then - - ! get the total volume of the instance - - minimumIpVolume = 1e99_pReal !use huge here? - totalVolume = 0.0_pReal - do e = 1_pInt,mesh_NcpElems - do i = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,e))) - if (PLASTICITY_NONLOCAL_ID == phase_plasticity(material_phase(1,i,e)) & - .and. instance == phase_plasticityInstance(material_phase(1,i,e))) then - totalVolume = totalVolume + mesh_ipVolume(i,e) - minimumIpVolume = min(minimumIpVolume, mesh_ipVolume(i,e)) - endif - enddo - enddo - densityBinning = rhoSglRandomBinning(instance) / minimumIpVolume ** (2.0_pReal / 3.0_pReal) - - ! subsequently fill random ips with dislocation segments until we reach the desired overall density - - meanDensity = 0.0_pReal - do while(meanDensity < rhoSglRandom(instance)) - call random_number(rnd) - el = nint(rnd(1)*real(mesh_NcpElems,pReal)+0.5_pReal,pInt) - ip = nint(rnd(2)*real(FE_Nips(FE_geomtype(mesh_element(2,el))),pReal)+0.5_pReal,pInt) - if (PLASTICITY_NONLOCAL_ID == phase_plasticity(material_phase(1,ip,el)) & - .and. instance == phase_plasticityInstance(material_phase(1,ip,el))) then - s = nint(rnd(3)*real(ns,pReal)+0.5_pReal,pInt) - t = nint(rnd(4)*4.0_pReal+0.5_pReal,pInt) - meanDensity = meanDensity + densityBinning * mesh_ipVolume(ip,el) / totalVolume - state(1,ip,el)%p(iRhoU(s,t,instance)) = state(1,ip,el)%p(iRhoU(s,t,instance)) + densityBinning - endif - enddo - ! homogeneous distribution of density with some noise - else - do e = 1_pInt,mesh_NcpElems - do i = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,e))) - if (PLASTICITY_NONLOCAL_ID == phase_plasticity(material_phase(1,i,e)) & - .and. instance == phase_plasticityInstance(material_phase(1,i,e))) then - do f = 1_pInt,lattice_maxNslipFamily - from = 1_pInt + sum(Nslip(1:f-1_pInt,instance)) - upto = sum(Nslip(1:f,instance)) - do s = from,upto - do j = 1_pInt,2_pInt - noise(j) = math_sampleGaussVar(0.0_pReal, rhoSglScatter(instance)) - enddo - state(1,i,e)%p(iRhoU(s,1,instance)) = rhoSglEdgePos0(f,instance) + noise(1) - state(1,i,e)%p(iRhoU(s,2,instance)) = rhoSglEdgeNeg0(f,instance) + noise(1) - state(1,i,e)%p(iRhoU(s,3,instance)) = rhoSglScrewPos0(f,instance) + noise(2) - state(1,i,e)%p(iRhoU(s,4,instance)) = rhoSglScrewNeg0(f,instance) + noise(2) - enddo - state(1,i,e)%p(iRhoD(from:upto,1,instance)) = rhoDipEdge0(f,instance) - state(1,i,e)%p(iRhoD(from:upto,2,instance)) = rhoDipScrew0(f,instance) - enddo - endif - enddo - enddo - endif -enddo - -end subroutine constitutive_nonlocal_stateInit - - !-------------------------------------------------------------------------------------------------- !> @brief sets the relevant state values for a given instance of this plasticity !-------------------------------------------------------------------------------------------------- @@ -1693,339 +1543,18 @@ forall (c = 1_pInt:2_pInt) & constitutive_nonlocal_aTolState(iGamma(1:ns,instance)) = aTolShear(instance) end function constitutive_nonlocal_aTolState - #endif + +!-------------------------------------------------------------------------------------------------- +!> @brief calculates quantities characterizing the microstructure +!-------------------------------------------------------------------------------------------------- #ifdef NEWSTATE -!-------------------------------------------------------------------------------------------------- -!> @brief calculates quantities characterizing the microstructure -!-------------------------------------------------------------------------------------------------- -subroutine constitutive_nonlocal_microstructure(mapping, Fe, Fp, ipc, ip, el) - -use IO, only: & - IO_error -use math, only: & - pi, & - math_mul33x3, & - math_mul3x3, & - math_norm3, & - math_invert33, & - math_transpose33 -use debug, only: & - debug_level, & - debug_constitutive, & - debug_levelBasic, & - debug_levelExtensive, & - debug_levelSelective, & - debug_g, & - debug_i, & - debug_e -use mesh, only: & - mesh_NcpElems, & - mesh_maxNips, & - mesh_element, & - mesh_ipNeighborhood, & - mesh_ipCoordinates, & - mesh_ipVolume, & - mesh_ipAreaNormal, & - mesh_ipArea, & - FE_NipNeighbors, & - mesh_maxNipNeighbors, & - FE_geomtype, & - FE_celltype -use material, only: & - homogenization_maxNgrains, & - material_phase, & - plasticState ,& - phase_localPlasticity, & - phase_plasticityInstance -use lattice, only: & - lattice_sd, & - lattice_st, & - lattice_mu, & - lattice_nu, & - lattice_structure, & - LATTICE_bcc_ID, & - LATTICE_fcc_ID - -implicit none - -!*** input variables -integer(pInt), intent(in) :: ipc, & ! current grain ID - ip, & ! current integration point - el ! current element -real(pReal), dimension(3,3), intent(in) :: & - Fe, & ! elastic deformation gradient - Fp ! elastic deformation gradient - - -!*** output variables - -!*** local variables -integer(pInt) neighbor_el, & ! element number of neighboring material point - neighbor_ip, & ! integration point of neighboring material point - instance, & ! my instance of this plasticity - neighbor_instance, & ! instance of this plasticity of neighboring material point - phase, & - neighbor_phase, & - ns, & ! total number of active slip systems at my material point - neighbor_ns, & ! total number of active slip systems at neighboring material point - c, & ! index of dilsocation character (edge, screw) - s, & ! slip system index - t, & ! index of dilsocation type (e+, e-, s+, s-, used e+, used e-, used s+, used s-) - dir, & - n, & - nRealNeighbors ! number of really existing neighbors -integer(pInt), dimension(2) :: neighbors -real(pReal) detFe, & - detFp, & - FVsize, & - temp, & - correction, & - myRhoForest -real(pReal), dimension(2) :: rhoExcessGradient, & - rhoExcessGradient_over_rho, & - rhoTotal -real(pReal), dimension(3) :: rhoExcessDifferences, & - normal_latticeConf -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & - rhoForest, & ! forest dislocation density - tauBack, & ! back stress from pileup on same slip system - tauThreshold ! threshold shear stress -real(pReal), dimension(3,3) :: invFe, & ! inverse of elastic deformation gradient - invFp, & ! inverse of plastic deformation gradient - connections, & - invConnections -real(pReal), dimension(3,mesh_maxNipNeighbors) :: & - connection_latticeConf -real(pReal), dimension(2,totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & - rhoExcess -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),2) :: & - rhoDip ! dipole dislocation density (edge, screw) -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),8) :: & - rhoSgl ! single dislocation density (edge+, edge-, screw+, screw-, used edge+, used edge-, used screw+, used screw-) -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))), & - totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & - myInteractionMatrix ! corrected slip interaction matrix -real(pReal), dimension(2,maxval(totalNslip),mesh_maxNipNeighbors) :: & - neighbor_rhoExcess, & ! excess density at neighboring material point - neighbor_rhoTotal ! total density at neighboring material point -real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),2) :: & - m ! direction of dislocation motion -logical inversionError -integer(pInt),dimension(2,homogenization_maxngrains,mesh_maxNips,mesh_ncpelems), intent(in) :: & - mapping - -phase = material_phase(ipc,ip,el) -instance = phase_plasticityInstance(phase) -ns = totalNslip(instance) - -!*** get basic states - -forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) - - rhoSgl(s,t) = max(plasticState(mapping(2,1,ip,el))%state(iRhoU(s,t,instance),mapping(1,1,ip,el)), & - 0.0_pReal) ! ensure positive single mobile densities - rhoSgl(s,t+4_pInt) = plasticState(mapping(2,1,ip,el))%state(iRhoB(s,t,instance),mapping(1,1,ip,el)) -endforall -forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) & - rhoDip(s,c) = max(plasticState(mapping(2,1,ip,el))%state(iRhoD(s,c,instance),mapping(1,1,ip,el)), 0.0_pReal) ! ensure positive dipole densities -where (abs(rhoSgl) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(instance) & - .or. abs(rhoSgl) < significantRho(instance)) & - rhoSgl = 0.0_pReal -where (abs(rhoDip) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(instance) & - .or. abs(rhoDip) < significantRho(instance)) & - rhoDip = 0.0_pReal - - -!*** calculate the forest dislocation density -!*** (= projection of screw and edge dislocations) - -forall (s = 1_pInt:ns) & - rhoForest(s) = dot_product((sum(abs(rhoSgl(1:ns,[1,2,5,6])),2) + rhoDip(1:ns,1)), & - forestProjectionEdge(s,1:ns,instance)) & - + dot_product((sum(abs(rhoSgl(1:ns,[3,4,7,8])),2) + rhoDip(1:ns,2)), & - forestProjectionScrew(s,1:ns,instance)) - -!*** calculate the threshold shear stress for dislocation slip -!*** coefficients are corrected for the line tension effect -!*** (see Kubin,Devincre,Hoc; 2008; Modeling dislocation storage rates and mean free paths in face-centered cubic crystals) - -myInteractionMatrix = 0.0_pReal -myInteractionMatrix(1:ns,1:ns) = interactionMatrixSlipSlip(1:ns,1:ns,instance) -if (lattice_structure(phase) == LATTICE_bcc_ID .or. lattice_structure(phase) == LATTICE_fcc_ID) then ! only fcc and bcc - do s = 1_pInt,ns - myRhoForest = max(rhoForest(s),significantRho(instance)) - correction = ( 1.0_pReal - linetensionEffect(instance) & - + linetensionEffect(instance) & - * log(0.35_pReal * burgers(s,instance) * sqrt(myRhoForest)) & - / log(0.35_pReal * burgers(s,instance) * 1e6_pReal)) ** 2.0_pReal - myInteractionMatrix(s,1:ns) = correction * myInteractionMatrix(s,1:ns) - enddo -endif -forall (s = 1_pInt:ns) & - tauThreshold(s) = lattice_mu(phase) * burgers(s,instance) & - * sqrt(dot_product((sum(abs(rhoSgl),2) + sum(abs(rhoDip),2)), myInteractionMatrix(s,1:ns))) - -!*** calculate the dislocation stress of the neighboring excess dislocation densities -!*** zero for material points of local plasticity - -tauBack = 0.0_pReal - -if (.not. phase_localPlasticity(phase) .and. shortRangeStressCorrection(instance)) then - call math_invert33(Fe, invFe, detFe, inversionError) - call math_invert33(Fp, invFp, detFp, inversionError) - rhoExcess(1,1:ns) = rhoSgl(1:ns,1) - rhoSgl(1:ns,2) - rhoExcess(2,1:ns) = rhoSgl(1:ns,3) - rhoSgl(1:ns,4) - FVsize = mesh_ipVolume(ip,el) ** (1.0_pReal/3.0_pReal) - - !* loop through my neighborhood and get the connection vectors (in lattice frame) and the excess densities - - nRealNeighbors = 0_pInt - neighbor_rhoTotal = 0.0_pReal - do n = 1_pInt,FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,el)))) - neighbor_el = mesh_ipNeighborhood(1,n,ip,el) - neighbor_ip = mesh_ipNeighborhood(2,n,ip,el) - if (neighbor_el > 0 .and. neighbor_ip > 0) then - neighbor_phase = material_phase(ipc,neighbor_ip,neighbor_el) - neighbor_instance = phase_plasticityInstance(neighbor_phase) - neighbor_ns = totalNslip(neighbor_instance) - if (.not. phase_localPlasticity(neighbor_phase) & - .and. neighbor_instance == instance) then ! same instance should be same structure - if (neighbor_ns == ns) then - nRealNeighbors = nRealNeighbors + 1_pInt - forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) - neighbor_rhoExcess(c,s,n) = & - max(plasticState(mapping(2,ipc,neighbor_ip,neighbor_el))% & - state(iRhoU(s,2*c-1,neighbor_instance), & - mapping(1,ipc,neighbor_ip,neighbor_el)),0.0_pReal) &! positive mobiles - - max(plasticState(mapping(2,ipc,neighbor_ip,neighbor_el))% & - state(iRhoU(s,2*c,neighbor_instance), & - mapping(1,ipc,neighbor_ip,neighbor_el)),0.0_pReal) ! negative mobiles - - neighbor_rhoTotal(c,s,n) = & - max(plasticState(mapping(2,ipc,neighbor_ip,neighbor_el))% & - state(iRhoU(s,2*c-1,neighbor_instance),mapping(1,ipc,neighbor_ip,neighbor_el)),& - 0.0_pReal) &! positive mobiles - + max(plasticState(mapping(2,ipc,neighbor_ip,neighbor_el))% & - state(iRhoU(s,2*c,neighbor_instance),mapping(1,ipc,neighbor_ip,neighbor_el)), & - 0.0_pReal) & ! negative mobiles - + abs(plasticState(mapping(2,ipc,neighbor_ip,neighbor_el))% & - state(iRhoB(s,2*c-1,neighbor_instance),mapping(1,ipc,neighbor_ip,neighbor_el)))& ! positive deads - + abs(plasticState(mapping(2,ipc,neighbor_ip,neighbor_el))% & - state(iRhoB(s,2*c,neighbor_instance),mapping(1,ipc,neighbor_ip,neighbor_el))) & ! negative deads - + max(plasticState(mapping(2,ipc,neighbor_ip,neighbor_el))% & - state(iRhoD(s,c,neighbor_instance),mapping(1,ipc,neighbor_ip,neighbor_el)), & - 0.0_pReal) ! dipoles - - endforall - connection_latticeConf(1:3,n) = & - math_mul33x3(invFe, mesh_ipCoordinates(1:3,neighbor_ip,neighbor_el) & - - mesh_ipCoordinates(1:3,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 ! neighboring connection points in opposite direction to face normal: must be periodic image - 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 - endif - else - ! different number of active slip systems - call IO_error(-1_pInt,ext_msg='different number of active slip systems in neighboring IPs of same crystal structure') - endif - else - ! local neighbor or different lattice structure or different constitution instance -> use central values instead - connection_latticeConf(1:3,n) = 0.0_pReal - neighbor_rhoExcess(1:2,1:ns,n) = rhoExcess - endif - else - ! free surface -> use central values instead - connection_latticeConf(1:3,n) = 0.0_pReal - neighbor_rhoExcess(1:2,1:ns,n) = rhoExcess - endif - enddo - - - !* loop through the slip systems and calculate the dislocation gradient by - !* 1. interpolation of the excess density in the neighorhood - !* 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),phase) - m(1:3,1:ns,2) = -lattice_st(1:3,slipSystemLattice(1:ns,instance),phase) - - do s = 1_pInt,ns - - !* gradient from interpolation of neighboring excess density - - do c = 1_pInt,2_pInt - do dir = 1_pInt,3_pInt - neighbors(1) = 2_pInt * dir - 1_pInt - neighbors(2) = 2_pInt * dir - connections(dir,1:3) = connection_latticeConf(1:3,neighbors(1)) & - - connection_latticeConf(1:3,neighbors(2)) - rhoExcessDifferences(dir) = neighbor_rhoExcess(c,s,neighbors(1)) & - - neighbor_rhoExcess(c,s,neighbors(2)) - enddo - call math_invert33(connections,invConnections,temp,inversionError) - if (inversionError) then - call IO_error(-1_pInt,ext_msg='back stress calculation: inversion error') - endif - rhoExcessGradient(c) = math_mul3x3(m(1:3,s,c), & - math_mul33x3(invConnections,rhoExcessDifferences)) - enddo - - !* plus gradient from deads - - do t = 1_pInt,4_pInt - c = (t - 1_pInt) / 2_pInt + 1_pInt - rhoExcessGradient(c) = rhoExcessGradient(c) + rhoSgl(s,t+4_pInt) / FVsize - enddo - - !* normalized with the total density - - rhoExcessGradient_over_rho = 0.0_pReal - 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(neighbor_rhoTotal(c,s,:))) / real(1_pInt + nRealNeighbors,pReal) - forall (c = 1_pInt:2_pInt, rhoTotal(c) > 0.0_pReal) & - rhoExcessGradient_over_rho(c) = rhoExcessGradient(c) / rhoTotal(c) - - !* gives the local stress correction when multiplied with a factor - - tauBack(s) = - lattice_mu(phase) * burgers(s,instance) / (2.0_pReal * pi) & - * (rhoExcessGradient_over_rho(1) / (1.0_pReal - lattice_nu(phase)) & - + rhoExcessGradient_over_rho(2)) - - enddo -endif - - -!*** set dependent states - - plasticState(mapping(2,ipc,ip,el))%state(iRhoF(1:ns,instance),mapping(1,ipc,ip,el)) = rhoForest - plasticState(mapping(2,ipc,ip,el))%state(iTauF(1:ns,instance),mapping(1,ipc,ip,el)) = tauThreshold - plasticState(mapping(2,ipc,ip,el))%state(iTauB(1:ns,instance),mapping(1,ipc,ip,el)) = tauBack - -#ifndef _OPENMP - if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt & - .and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)& - .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then - write(6,*) - write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_microstructure at el ip g',el,ip,ipc - write(6,*) - write(6,'(a,/,12x,12(e10.3,1x))') '<< CONST >> rhoForest', rhoForest - write(6,'(a,/,12x,12(f10.5,1x))') '<< CONST >> tauThreshold / MPa', tauThreshold/1e6 - write(6,'(a,/,12x,12(f10.5,1x))') '<< CONST >> tauBack / MPa', tauBack/1e6 - write(6,*) - endif +subroutine constitutive_nonlocal_microstructure(Fe, Fp, ipc, ip, el) +#else +subroutine constitutive_nonlocal_microstructure(state, Fe, Fp, ipc, ip, el) #endif -end subroutine constitutive_nonlocal_microstructure - -#else -!-------------------------------------------------------------------------------------------------- -!> @brief calculates quantities characterizing the microstructure -!-------------------------------------------------------------------------------------------------- -subroutine constitutive_nonlocal_microstructure(state, Fe, Fp, ipc, ip, el) - use IO, only: & IO_error use math, only: & @@ -2061,6 +1590,10 @@ use material, only: & homogenization_maxNgrains, & material_phase, & phase_localPlasticity, & +#ifdef NEWSTATE + plasticState, & + mappingConstitutive, & +#endif phase_plasticityInstance use lattice, only: & lattice_sd, & @@ -2081,11 +1614,16 @@ real(pReal), dimension(3,3), intent(in) :: & Fe, & ! elastic deformation gradient Fp ! elastic deformation gradient -!*** input/output variables +#ifndef NEWSTATE type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(inout) :: & state ! microstructural state - -!*** output variables +#else + integer(pInt) :: & + p, & !< phase + o, & !< offset + np, & !< neighbor phase + no !< nieghbor offset +#endif !*** local variables integer(pInt) neighbor_el, & ! element number of neighboring material point @@ -2145,13 +1683,23 @@ instance = phase_plasticityInstance(phase) ns = totalNslip(instance) !*** get basic states - +#ifdef NEWSTATE +p = mappingConstitutive(2,ipc,ip,el) +o = mappingConstitutive(1,ipc,ip,el) +forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) + rhoSgl(s,t) = max(plasticState(p)%state(iRhoU(s,t,instance),o), 0.0_pReal) ! ensure positive single mobile densities + rhoSgl(s,t+4_pInt) = plasticState(p)%state(iRhoB(s,t,instance),o) +endforall +forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) & + rhoDip(s,c) = max(plasticState(p)%state(iRhoD(s,c,instance),o), 0.0_pReal) ! ensure positive dipole densities +#else 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+4_pInt) = state(ipc,ip,el)%p(iRhoB(s,t,instance)) endforall forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) & rhoDip(s,c) = max(state(ipc,ip,el)%p(iRhoD(s,c,instance)), 0.0_pReal) ! ensure positive dipole densities +#endif where (abs(rhoSgl) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(instance) & .or. abs(rhoSgl) < significantRho(instance)) & rhoSgl = 0.0_pReal @@ -2159,7 +1707,6 @@ where (abs(rhoDip) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(instance .or. abs(rhoDip) < significantRho(instance)) & rhoDip = 0.0_pReal - !*** calculate the forest dislocation density !*** (= projection of screw and edge dislocations) @@ -2210,6 +1757,10 @@ if (.not. phase_localPlasticity(phase) .and. shortRangeStressCorrection(instance do n = 1_pInt,FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,el)))) neighbor_el = mesh_ipNeighborhood(1,n,ip,el) neighbor_ip = mesh_ipNeighborhood(2,n,ip,el) +#ifdef NEWSTATE + np = mappingConstitutive(2,ipc,neighbor_ip,neighbor_el) + no = mappingConstitutive(1,ipc,neighbor_ip,neighbor_el) +#endif if (neighbor_el > 0 .and. neighbor_ip > 0) then neighbor_phase = material_phase(ipc,neighbor_ip,neighbor_el) neighbor_instance = phase_plasticityInstance(neighbor_phase) @@ -2219,6 +1770,17 @@ if (.not. phase_localPlasticity(phase) .and. shortRangeStressCorrection(instance if (neighbor_ns == ns) then nRealNeighbors = nRealNeighbors + 1_pInt forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) +#ifdef NEWSTATE + neighbor_rhoExcess(c,s,n) = & + max(plasticState(np)%state(iRhoU(s,2*c-1,neighbor_instance),no), 0.0_pReal) & ! positive mobiles + - max(plasticState(np)%state(iRhoU(s,2*c,neighbor_instance), no), 0.0_pReal) ! negative mobiles + neighbor_rhoTotal(c,s,n) = & + max(plasticState(np)%state(iRhoU(s,2*c-1,neighbor_instance),no), 0.0_pReal) & ! positive mobiles + + max(plasticState(np)%state(iRhoU(s,2*c,neighbor_instance), no), 0.0_pReal) & ! negative mobiles + + abs(plasticState(np)%state(iRhoB(s,2*c-1,neighbor_instance),no)) & ! positive deads + + abs(plasticState(np)%state(iRhoB(s,2*c,neighbor_instance), no)) & ! negative deads + + max(plasticState(np)%state(iRhoD(s,c,neighbor_instance), no), 0.0_pReal) ! dipoles +#else neighbor_rhoExcess(c,s,n) = & max(state(ipc,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c-1,neighbor_instance)), 0.0_pReal) &! positive mobiles - max(state(ipc,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c,neighbor_instance)), 0.0_pReal) ! negative mobiles @@ -2228,6 +1790,7 @@ if (.not. phase_localPlasticity(phase) .and. shortRangeStressCorrection(instance + abs(state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,2*c-1,neighbor_instance))) & ! positive deads + abs(state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,2*c,neighbor_instance))) & ! negative deads + max(state(ipc,neighbor_ip,neighbor_el)%p(iRhoD(s,c,neighbor_instance)), 0.0_pReal) ! dipoles +#endif endforall connection_latticeConf(1:3,n) = & math_mul33x3(invFe, mesh_ipCoordinates(1:3,neighbor_ip,neighbor_el) & @@ -2309,11 +1872,15 @@ endif !*** set dependent states - +#ifdef NEWSTATE +plasticState(p)%state(iRhoF(1:ns,instance),o) = rhoForest +plasticState(p)%state(iTauF(1:ns,instance),o) = tauThreshold +plasticState(p)%state(iTauB(1:ns,instance),o) = tauBack +#else state(ipc,ip,el)%p(iRhoF(1:ns,instance)) = rhoForest state(ipc,ip,el)%p(iTauF(1:ns,instance)) = tauThreshold state(ipc,ip,el)%p(iTauB(1:ns,instance)) = tauBack - +#endif #ifndef _OPENMP if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt & @@ -2331,7 +1898,6 @@ state(ipc,ip,el)%p(iTauB(1:ns,instance)) = tauBack end subroutine constitutive_nonlocal_microstructure -#endif !-------------------------------------------------------------------------------------------------- !> @brief calculates kinetics @@ -2495,217 +2061,14 @@ endif end subroutine constitutive_nonlocal_kinetics +!-------------------------------------------------------------------------------------------------- +!> @brief calculates plastic velocity gradient and its tangent +!-------------------------------------------------------------------------------------------------- #ifdef NEWSTATE -!-------------------------------------------------------------------------------------------------- -!> @brief calculates plastic velocity gradient and its tangent -!-------------------------------------------------------------------------------------------------- -subroutine constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar99, Tstar_v, Temperature, mapping, ipc, ip, el) - -use math, only: math_Plain3333to99, & - math_mul6x6, & - math_mul33xx33, & - math_Mandel6to33 -use debug, only: debug_level, & - debug_constitutive, & - debug_levelBasic, & - debug_levelExtensive, & - debug_levelSelective, & - debug_g, & - debug_i, & - debug_e -use material, only: material_phase, & - phase_plasticityInstance, & - plasticState, & - homogenization_maxngrains -use lattice, only: lattice_Sslip, & - lattice_Sslip_v, & - lattice_NnonSchmid -use mesh, only: mesh_ipVolume, & - mesh_maxNips, & - mesh_ncpelems - -implicit none - -!*** input variables -integer(pInt), intent(in) :: ipc, & !< current grain number - ip, & !< current integration point - el !< current element number -real(pReal), intent(in) :: Temperature !< temperature -real(pReal), dimension(6), intent(in) :: Tstar_v !< 2nd Piola-Kirchhoff stress in Mandel notation - -!*** input/output variables - -!*** output variables -real(pReal), dimension(3,3), intent(out) :: Lp !< plastic velocity gradient -real(pReal), dimension(9,9), intent(out) :: dLp_dTstar99 !< derivative of Lp with respect to Tstar (9x9 matrix) - -!*** local variables -integer(pInt) instance, & !< current instance of this plasticity - phase, & !< phase - ns, & !< short notation for the total number of active slip systems - i, & - j, & - k, & - l, & - t, & !< dislocation type - s, & !< index of my current slip system - sLattice !< index of my current slip system according to lattice order -real(pReal), dimension(3,3,3,3) :: dLp_dTstar3333 !< derivative of Lp with respect to Tstar (3x3x3x3 matrix) -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),8) :: & - rhoSgl !< single dislocation densities (including blocked) -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),4) :: & - v, & !< velocity - tauNS, & !< resolved shear stress including non Schmid and backstress terms - dv_dtau, & !< velocity derivative with respect to the shear stress - dv_dtauNS !< velocity derivative with respect to the shear stress -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & - tau, & !< resolved shear stress including backstress terms - gdotTotal, & !< shear rate - tauBack, & !< back stress from dislocation gradients on same slip system - tauThreshold !< threshold shear stress - -integer(pInt),dimension(2,homogenization_maxngrains,mesh_maxNips,mesh_ncpelems), intent(in) :: & - mapping - -!*** initialize local variables - -Lp = 0.0_pReal -dLp_dTstar3333 = 0.0_pReal - -phase = material_phase(ipc,ip,el) -instance = phase_plasticityInstance(phase) -ns = totalNslip(instance) - - -!*** shortcut to state variables - - -forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) - rhoSgl(s,t) = max(plasticState(mapping(2,ipc,ip,el))%state(iRhoU(s,t,instance),mapping(1,ipc,ip,el)), & - 0.0_pReal) ! ensure positive single mobile densities - rhoSgl(s,t+4_pInt) = plasticState(mapping(2,ipc,ip,el))%state(iRhoB(s,t,instance),mapping(1,ipc,ip,el)) -endforall -where (abs(rhoSgl) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(instance) & - .or. abs(rhoSgl) < significantRho(instance)) & - rhoSgl = 0.0_pReal - -tauBack = plasticState(mapping(2,ipc,ip,el))%state(iTauB(1:ns,instance),mapping(1,ipc,ip,el)) -tauThreshold = plasticState(mapping(2,ipc,ip,el))%state(iTauF(1:ns,instance),mapping(1,ipc,ip,el)) - - -!*** get resolved shear stress -!*** for screws possible non-schmid contributions are also taken into account - -do s = 1_pInt,ns - sLattice = slipSystemLattice(s,instance) - tau(s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,1,sLattice,phase)) - tauNS(s,1) = tau(s) - tauNS(s,2) = tau(s) - if (tau(s) > 0.0_pReal) then - tauNS(s,3) = math_mul33xx33(math_Mandel6to33(Tstar_v), nonSchmidProjection(1:3,1:3,1,s,instance)) - tauNS(s,4) = math_mul33xx33(math_Mandel6to33(Tstar_v), nonSchmidProjection(1:3,1:3,3,s,instance)) - else - tauNS(s,3) = math_mul33xx33(math_Mandel6to33(Tstar_v), nonSchmidProjection(1:3,1:3,2,s,instance)) - tauNS(s,4) = math_mul33xx33(math_Mandel6to33(Tstar_v), nonSchmidProjection(1:3,1:3,4,s,instance)) - endif -enddo -forall (t = 1_pInt:4_pInt) & - tauNS(1:ns,t) = tauNS(1:ns,t) + tauBack ! add backstress -tau = tau + tauBack ! add backstress - - -!*** get dislocation velocity and its tangent and store the velocity in the state array - -! edges -call constitutive_nonlocal_kinetics(v(1:ns,1), dv_dtau(1:ns,1), dv_dtauNS(1:ns,1), & - tau(1:ns), tauNS(1:ns,1), tauThreshold(1:ns), & - 1_pInt, Temperature, ipc, ip, el) -v(1:ns,2) = v(1:ns,1) -dv_dtau(1:ns,2) = dv_dtau(1:ns,1) -dv_dtauNS(1:ns,2) = dv_dtauNS(1:ns,1) - -!screws -if (lattice_NnonSchmid(phase) == 0_pInt) then ! no non-Schmid contributions - forall(t = 3_pInt:4_pInt) - v(1:ns,t) = v(1:ns,1) - dv_dtau(1:ns,t) = dv_dtau(1:ns,1) - dv_dtauNS(1:ns,t) = dv_dtauNS(1:ns,1) - endforall -else ! take non-Schmid contributions into account - do t = 3_pInt,4_pInt - call constitutive_nonlocal_kinetics(v(1:ns,t), dv_dtau(1:ns,t), dv_dtauNS(1:ns,t), & - tau(1:ns), tauNS(1:ns,t), tauThreshold(1:ns), & - 2_pInt , Temperature, ipc, ip, el) - enddo -endif - - -!*** store velocity in state - -forall (t = 1_pInt:4_pInt) & - plasticState(mapping(2,ipc,ip,el))%state(iV(1:ns,t,instance),mapping(1,ipc,ip,el)) = v(1:ns,t) - - -!*** Bauschinger effect - -forall (s = 1_pInt:ns, t = 5_pInt:8_pInt, rhoSgl(s,t) * v(s,t-4_pInt) < 0.0_pReal) & - rhoSgl(s,t-4_pInt) = rhoSgl(s,t-4_pInt) + abs(rhoSgl(s,t)) - - -!*** Calculation of Lp and its tangent - -gdotTotal = sum(rhoSgl(1:ns,1:4) * v, 2) * burgers(1:ns,instance) - -do s = 1_pInt,ns - sLattice = slipSystemLattice(s,instance) - Lp = Lp + gdotTotal(s) * lattice_Sslip(1:3,1:3,1,sLattice,phase) - - ! Schmid contributions to tangent - forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt,k=1_pInt:3_pInt,l=1_pInt:3_pInt) & - dLp_dTstar3333(i,j,k,l) = dLp_dTstar3333(i,j,k,l) & - + lattice_Sslip(i,j,1,sLattice,phase) * lattice_Sslip(k,l,1,sLattice,phase) & - * sum(rhoSgl(s,1:4) * dv_dtau(s,1:4)) * burgers(s,instance) - - ! non Schmid contributions to tangent - if (tau(s) > 0.0_pReal) then - forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt,k=1_pInt:3_pInt,l=1_pInt:3_pInt) & - dLp_dTstar3333(i,j,k,l) = dLp_dTstar3333(i,j,k,l) & - + lattice_Sslip(i,j,1,sLattice,phase) & - * ( nonSchmidProjection(k,l,1,s,instance) * rhoSgl(s,3) * dv_dtauNS(s,3) & - + nonSchmidProjection(k,l,3,s,instance) * rhoSgl(s,4) * dv_dtauNS(s,4) ) & - * burgers(s,instance) - else - forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt,k=1_pInt:3_pInt,l=1_pInt:3_pInt) & - dLp_dTstar3333(i,j,k,l) = dLp_dTstar3333(i,j,k,l) & - + lattice_Sslip(i,j,1,sLattice,phase) & - * ( nonSchmidProjection(k,l,2,s,instance) * rhoSgl(s,3) * dv_dtauNS(s,3) & - + nonSchmidProjection(k,l,4,s,instance) * rhoSgl(s,4) * dv_dtauNS(s,4) ) & - * burgers(s,instance) - endif -enddo -dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333) - - -#ifndef _OPENMP - if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt & - .and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)& - .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then - write(6,*) - write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_LpandItsTangent at el ip ipc ',el,ip,ipc - write(6,*) - write(6,'(a,/,12x,12(f12.5,1x))') '<< CONST >> gdot total / 1e-3',gdotTotal*1e3_pReal - write(6,'(a,/,3(12x,3(f12.7,1x),/))') '<< CONST >> Lp',transpose(Lp) - endif -#endif - -end subroutine constitutive_nonlocal_LpAndItsTangent - +subroutine constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar99, Tstar_v, Temperature,ipc, ip, el) #else -!-------------------------------------------------------------------------------------------------- -!> @brief calculates plastic velocity gradient and its tangent -!-------------------------------------------------------------------------------------------------- subroutine constitutive_nonlocal_LpAndItsTangent(Lp, dLp_dTstar99, Tstar_v, Temperature, state, ipc, ip, el) - +#endif use math, only: math_Plain3333to99, & math_mul6x6, & math_mul33xx33, & @@ -2719,6 +2082,10 @@ use debug, only: debug_level, & debug_i, & debug_e use material, only: material_phase, & +#ifdef NEWSTATE + plasticState, & + mappingConstitutive,& +#endif phase_plasticityInstance use lattice, only: lattice_Sslip, & lattice_Sslip_v, & @@ -2728,45 +2095,51 @@ use mesh, only: mesh_ipVolume implicit none !*** input variables -integer(pInt), intent(in) :: ipc, & !< current grain number - ip, & !< current integration point - el !< current element number -real(pReal), intent(in) :: Temperature !< temperature -real(pReal), dimension(6), intent(in) :: Tstar_v !< 2nd Piola-Kirchhoff stress in Mandel notation +integer(pInt), intent(in) :: ipc, & !< current grain number + ip, & !< current integration point + el !< current element number +real(pReal), intent(in) :: Temperature !< temperature +real(pReal), dimension(6), intent(in) :: Tstar_v !< 2nd Piola-Kirchhoff stress in Mandel notation !*** input/output variables -type(p_vec), intent(inout) :: state !< microstructural state - +#ifndef NEWSTATE +type(p_vec), intent(inout) :: state !< microstructural state +#endif !*** output variables -real(pReal), dimension(3,3), intent(out) :: Lp !< plastic velocity gradient -real(pReal), dimension(9,9), intent(out) :: dLp_dTstar99 !< derivative of Lp with respect to Tstar (9x9 matrix) +real(pReal), dimension(3,3), intent(out) :: Lp !< plastic velocity gradient +real(pReal), dimension(9,9), intent(out) :: dLp_dTstar99 !< derivative of Lp with respect to Tstar (9x9 matrix) !*** local variables -integer(pInt) instance, & !< current instance of this plasticity - phase, & !< phase - ns, & !< short notation for the total number of active slip systems +integer(pInt) instance, & !< current instance of this plasticity + phase, & !< phase + ns, & !< short notation for the total number of active slip systems i, & j, & k, & l, & - t, & !< dislocation type - s, & !< index of my current slip system - sLattice !< index of my current slip system according to lattice order -real(pReal), dimension(3,3,3,3) :: dLp_dTstar3333 !< derivative of Lp with respect to Tstar (3x3x3x3 matrix) + p, & !phase number + o, & !offset + t, & !< dislocation type + s, & !< index of my current slip system + sLattice !< index of my current slip system according to lattice order +real(pReal), dimension(3,3,3,3) :: dLp_dTstar3333 !< derivative of Lp with respect to Tstar (3x3x3x3 matrix) real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),8) :: & - rhoSgl !< single dislocation densities (including blocked) + rhoSgl !< single dislocation densities (including blocked) real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),4) :: & - v, & !< velocity - tauNS, & !< resolved shear stress including non Schmid and backstress terms - dv_dtau, & !< velocity derivative with respect to the shear stress - dv_dtauNS !< velocity derivative with respect to the shear stress + v, & !< velocity + tauNS, & !< resolved shear stress including non Schmid and backstress terms + dv_dtau, & !< velocity derivative with respect to the shear stress + dv_dtauNS !< velocity derivative with respect to the shear stress real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & - tau, & !< resolved shear stress including backstress terms - gdotTotal, & !< shear rate - tauBack, & !< back stress from dislocation gradients on same slip system - tauThreshold !< threshold shear stress - - + tau, & !< resolved shear stress including backstress terms + gdotTotal, & !< shear rate + tauBack, & !< back stress from dislocation gradients on same slip system + tauThreshold !< threshold shear stress +#ifdef NEWSTATE +!*** shortcut for mapping +p = mappingConstitutive(2,ipc,ip,el) +o = mappingConstitutive(1,ipc,ip,el) +#endif !*** initialize local variables Lp = 0.0_pReal @@ -2781,16 +2154,25 @@ ns = totalNslip(instance) forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) - rhoSgl(s,t) = max(state%p(iRhoU(s,t,instance)), 0.0_pReal) ! ensure positive single mobile densities +#ifdef NEWSTATE + rhoSgl(s,t) = max(plasticState(p)%state(iRhoU(s,t,instance),o), 0.0_pReal) ! ensure positive single mobile densities + rhoSgl(s,t+4_pInt) = plasticState(p)%state(iRhoB(s,t,instance),o) +#else + rhoSgl(s,t) = max(state%p(iRhoU(s,t,instance)), 0.0_pReal) ! ensure positive single mobile densities rhoSgl(s,t+4_pInt) = state%p(iRhoB(s,t,instance)) +#endif endforall where (abs(rhoSgl) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(instance) & .or. abs(rhoSgl) < significantRho(instance)) & rhoSgl = 0.0_pReal - + +#ifdef NEWSTATE +tauBack = plasticState(p)%state(iTauB(1:ns,instance),o) +tauThreshold = plasticState(p)%state(iTauF(1:ns,instance),o) +#else tauBack = state%p(iTauB(1:ns,instance)) tauThreshold = state%p(iTauF(1:ns,instance)) - +#endif !*** get resolved shear stress !*** for screws possible non-schmid contributions are also taken into account @@ -2842,8 +2224,11 @@ endif !*** store velocity in state forall (t = 1_pInt:4_pInt) & +#ifdef NEWSTATE + plasticState(p)%state(iV(1:ns,t,instance),o) = v(1:ns,t) +#else state%p(iV(1:ns,t,instance)) = v(1:ns,t) - +#endif !*** Bauschinger effect @@ -2898,205 +2283,17 @@ dLp_dTstar99 = math_Plain3333to99(dLp_dTstar3333) #endif end subroutine constitutive_nonlocal_LpAndItsTangent -#endif + + +!-------------------------------------------------------------------------------------------------- +!> @brief (instantaneous) incremental change of microstructure +!-------------------------------------------------------------------------------------------------- #ifdef NEWSTATE -!-------------------------------------------------------------------------------------------------- -!> @brief (instantaneous) incremental change of microstructure -!-------------------------------------------------------------------------------------------------- -subroutine constitutive_nonlocal_deltaState(mapping, Tstar_v, ipc,ip,el) - -use debug, only: debug_level, & - debug_constitutive, & - debug_levelBasic, & - debug_levelExtensive, & - debug_levelSelective, & - debug_g, & - debug_i, & - debug_e -use math, only: pi, & - math_mul6x6 -use lattice, only: lattice_Sslip_v ,& - lattice_mu, & - lattice_nu -use mesh, only: mesh_NcpElems, & - mesh_maxNips, & - mesh_ipVolume -use material, only: homogenization_maxNgrains, & - material_phase, & - phase_plasticityInstance, & - plasticState - -implicit none - -!*** input variables -integer(pInt), intent(in) :: ipc, & ! current grain number - ip, & ! current integration point - el ! current element number -real(pReal), dimension(6), intent(in) :: Tstar_v ! current 2nd Piola-Kirchhoff stress in Mandel notation - -!*** local variables -integer(pInt) phase, & - instance, & ! current instance of this plasticity - ns, & ! short notation for the total number of active slip systems - c, & ! character of dislocation - t, & ! type of dislocation - s, & ! index of my current slip system - sLattice ! index of my current slip system according to lattice order -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),10) :: & - deltaRho, & ! density increment - deltaRhoRemobilization, & ! density increment by remobilization - deltaRhoDipole2SingleStress ! density increment by dipole dissociation (by stress change) -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),8) :: & - rhoSgl ! current single dislocation densities (positive/negative screw and edge without dipoles) -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),4) :: & - v ! dislocation glide velocity -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & - tau, & ! current resolved shear stress - tauBack ! current back stress from pileups on same slip system -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),2) :: & - rhoDip, & ! current dipole dislocation densities (screw and edge dipoles) - dLower, & ! minimum stable dipole distance for edges and screws - dUpper, & ! current maximum stable dipole distance for edges and screws - dUpperOld, & ! old maximum stable dipole distance for edges and screws - deltaDUpper ! change in maximum stable dipole distance for edges and screws -integer(pInt),dimension(2,homogenization_maxngrains,mesh_maxNips,mesh_ncpelems), intent(in) :: & - mapping -integer(pInt) :: sizeDot -#ifndef _OPENMP - if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt & - .and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)& - .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then - write(6,*) - write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_deltaState at el ip ipc ',el,ip,ipc - write(6,*) - endif -#endif - -phase = material_phase(ipc,ip,el) -instance = phase_plasticityInstance(phase) -ns = totalNslip(instance) - - -!*** shortcut to state variables - -forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) - rhoSgl(s,t) = max(plasticState(mapping(2,ipc,ip,el))%state(iRhoU(s,t,instance), & - mapping(1,ipc,ip,el)),0.0_pReal) ! ensure positive single mobile densities - rhoSgl(s,t+4_pInt) = plasticState(mapping(2,ipc,ip,el))%state(iRhoB(s,t,instance), & - mapping(1,ipc,ip,el)) - v(s,t) = plasticState(mapping(2,ipc,ip,el))%state(iV(s,t,instance),mapping(1,ipc,ip,el)) -endforall -forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) - rhoDip(s,c) = max(plasticState(mapping(2,ipc,ip,el))%state(iRhoD(s,c,instance),mapping(1,ipc,ip,el)), 0.0_pReal) ! ensure positive dipole densities - dUpperOld(s,c) = plasticState(mapping(2,ipc,ip,el))%state(iD(s,c,instance),mapping(1,ipc,ip,el)) - -endforall -tauBack = plasticState(mapping(2,ipc,ip,el))%state(iTauB(1:ns,instance),mapping(1,ipc,ip,el)) - -where (abs(rhoSgl) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(instance) & - .or. abs(rhoSgl) < significantRho(instance)) & - rhoSgl = 0.0_pReal -where (abs(rhoDip) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(instance) & - .or. abs(rhoDip) < significantRho(instance)) & - rhoDip = 0.0_pReal - - -!**************************************************************************** -!*** dislocation remobilization (bauschinger effect) - -deltaRhoRemobilization = 0.0_pReal -do t = 1_pInt,4_pInt - do s = 1_pInt,ns - if (rhoSgl(s,t+4_pInt) * v(s,t) < 0.0_pReal) then - deltaRhoRemobilization(s,t) = abs(rhoSgl(s,t+4_pInt)) - rhoSgl(s,t) = rhoSgl(s,t) + abs(rhoSgl(s,t+4_pInt)) - deltaRhoRemobilization(s,t+4_pInt) = - rhoSgl(s,t+4_pInt) - rhoSgl(s,t+4_pInt) = 0.0_pReal - endif - enddo -enddo - - -!**************************************************************************** -!*** calculate dipole formation and dissociation by stress change - -!*** calculate limits for stable dipole height - -do s = 1_pInt,ns - sLattice = slipSystemLattice(s,instance) - tau(s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,1,sLattice,phase)) + tauBack(s) - if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal -enddo -dLower = minDipoleHeight(1:ns,1:2,instance) -dUpper(1:ns,1) = lattice_mu(phase) * burgers(1:ns,instance) & - / (8.0_pReal * pi * (1.0_pReal - lattice_nu(phase)) * abs(tau)) -dUpper(1:ns,2) = lattice_mu(phase) * burgers(1:ns,instance) / (4.0_pReal * pi * abs(tau)) - - -!in the test there is an FPE exception. Division by zero? -forall (c = 1_pInt:2_pInt) & - dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) & - + abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), & - dUpper(1:ns,c)) -dUpper = max(dUpper,dLower) -deltaDUpper = dUpper - dUpperOld - -!*** dissociation by stress increase - -deltaRhoDipole2SingleStress = 0.0_pReal -forall (c=1_pInt:2_pInt, s=1_pInt:ns, deltaDUpper(s,c) < 0.0_pReal) & - deltaRhoDipole2SingleStress(s,8_pInt+c) = rhoDip(s,c) * deltaDUpper(s,c) & - / (dUpperOld(s,c) - dLower(s,c)) - -forall (t=1_pInt:4_pInt) & - deltaRhoDipole2SingleStress(1_pInt:ns,t) = -0.5_pReal & - * deltaRhoDipole2SingleStress(1_pInt:ns,(t-1_pInt)/2_pInt+9_pInt) - - - -!*** store new maximum dipole height in state - -forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) & - - plasticState(mapping(2,ipc,ip,el))%state(iD(s,c,instance),mapping(1,ipc,ip,el)) = dUpper(s,c) - -!**************************************************************************** -!*** assign the changes in the dislocation densities to deltaState - -deltaRho = deltaRhoRemobilization & - + deltaRhoDipole2SingleStress - -sizeDot = plasticState(mapping(2,ipc,ip,el))%sizeDotState -!deltaState%p = 0.0_pReal -plasticState(mapping(2,ipc,ip,el))%deltaState(1:sizeDot,mapping(1,ipc,ip,el)) = 0.0_pReal -forall (s = 1:ns, t = 1_pInt:4_pInt) - plasticState(mapping(2,ipc,ip,el))%deltaState(iRhoU(s,t,instance),mapping(1,ipc,ip,el)) = & - deltaRho(s,t) - plasticState(mapping(2,ipc,ip,el))%deltaState(iRhoB(s,t,instance),mapping(1,ipc,ip,el)) = & - deltaRho(s,t+4_pInt) -endforall -forall (s = 1:ns, c = 1_pInt:2_pInt) & - plasticState(mapping(2,ipc,ip,el))%deltaState(iRhoD(s,c,instance),mapping(1,ipc,ip,el)) = & - deltaRho(s,c+8_pInt) - -#ifndef _OPENMP - if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt & - .and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)& - .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then - write(6,'(a,/,8(12x,12(e12.5,1x),/))') '<< CONST >> dislocation remobilization', deltaRhoRemobilization(1:ns,1:8) - write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> dipole dissociation by stress increase', deltaRhoDipole2SingleStress - write(6,*) - endif -#endif - -end subroutine constitutive_nonlocal_deltaState - +subroutine constitutive_nonlocal_deltaState(Tstar_v,ip,el) #else -!-------------------------------------------------------------------------------------------------- -!> @brief (instantaneous) incremental change of microstructure -!-------------------------------------------------------------------------------------------------- subroutine constitutive_nonlocal_deltaState(deltaState, state, Tstar_v, ipc,ip,el) +#endif use debug, only: debug_level, & debug_constitutive, & @@ -3116,22 +2313,36 @@ use mesh, only: mesh_NcpElems, & mesh_ipVolume use material, only: homogenization_maxNgrains, & material_phase, & +#ifdef NEWSTATE + plasticState, & + mappingConstitutive, & +#endif phase_plasticityInstance implicit none !*** input variables +#ifndef NEWSTATE integer(pInt), intent(in) :: ipc, & ! current grain number ip, & ! current integration point +#else +integer(pInt), intent(in) :: ip, & ! current grain number +#endif el ! current element number real(pReal), dimension(6), intent(in) :: Tstar_v ! current 2nd Piola-Kirchhoff stress in Mandel notation !*** input/output variables +#ifndef NEWSTATE type(p_vec), intent(inout) :: & state ! current microstructural state - -!*** output variables type(p_vec), intent(out) :: deltaState ! change of state variables / microstructure +#else + integer(pInt) :: & + p, & !< phase + o !< offset +#endif +!*** output variables + !*** local variables integer(pInt) phase, & @@ -3141,18 +2352,18 @@ integer(pInt) phase, & t, & ! type of dislocation s, & ! index of my current slip system sLattice ! index of my current slip system according to lattice order -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),10) :: & +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),10) :: & deltaRho, & ! density increment deltaRhoRemobilization, & ! density increment by remobilization deltaRhoDipole2SingleStress ! density increment by dipole dissociation (by stress change) -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),8) :: & +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),8) :: & rhoSgl ! current single dislocation densities (positive/negative screw and edge without dipoles) -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),4) :: & +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),4) :: & v ! dislocation glide velocity -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el)))) :: & tau, & ! current resolved shear stress tauBack ! current back stress from pileups on same slip system -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),2) :: & +real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(1,ip,el))),2) :: & rhoDip, & ! current dipole dislocation densities (screw and edge dipoles) dLower, & ! minimum stable dipole distance for edges and screws dUpper, & ! current maximum stable dipole distance for edges and screws @@ -3162,21 +2373,35 @@ real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip #ifndef _OPENMP if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt & - .and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)& + .and. ((debug_e == el .and. debug_i == ip .and. debug_g == 1)& .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then write(6,*) - write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_deltaState at el ip ipc ',el,ip,ipc + write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_deltaState at el ip ipc ',el,ip,1 write(6,*) endif #endif -phase = material_phase(ipc,ip,el) +phase = material_phase(1,ip,el) instance = phase_plasticityInstance(phase) ns = totalNslip(instance) !*** shortcut to state variables - +#ifdef NEWSTATE + p = mappingConstitutive(2,1,ip,el) + o = mappingConstitutive(1,1,ip,el) + + forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) + rhoSgl(s,t) = max(plasticState(p)%state(iRhoU(s,t,instance),o), 0.0_pReal) ! ensure positive single mobile densities + rhoSgl(s,t+4_pInt) = plasticState(p)%state(iRhoB(s,t,instance),o) + v(s,t) = plasticState(p)%state(iV(s,t,instance),o) +endforall +forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) + rhoDip(s,c) = max(plasticState(p)%state(iRhoD(s,c,instance),o), 0.0_pReal) ! ensure positive dipole densities + dUpperOld(s,c) = plasticState(p)%state(iD(s,c,instance),o) +endforall + tauBack = plasticState(p)%state(iTauB(1:ns,instance),o) +#else forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) rhoSgl(s,t) = max(state%p(iRhoU(s,t,instance)), 0.0_pReal) ! ensure positive single mobile densities rhoSgl(s,t+4_pInt) = state%p(iRhoB(s,t,instance)) @@ -3187,6 +2412,7 @@ forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) dUpperOld(s,c) = state%p(iD(s,c,instance)) endforall tauBack = state%p(iTauB(1:ns,instance)) +#endif where (abs(rhoSgl) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(instance) & .or. abs(rhoSgl) < significantRho(instance)) & @@ -3256,8 +2482,11 @@ forall (t=1_pInt:4_pInt) & !*** store new maximum dipole height in state forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) & +#ifdef NEWSTATE + plasticState(p)%state(iD(s,c,instance),o) = dUpper(s,c) +#else state%p(iD(s,c,instance)) = dUpper(s,c) - +#endif !**************************************************************************** @@ -3266,6 +2495,15 @@ forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) & deltaRho = deltaRhoRemobilization & + deltaRhoDipole2SingleStress +#ifdef NEWSTATE +plasticState(p)%deltaState(:,o) = 0.0_pReal +forall (s = 1:ns, t = 1_pInt:4_pInt) + plasticState(p)%deltaState(iRhoU(s,t,instance),o)= deltaRho(s,t) + plasticState(p)%deltaState(iRhoB(s,t,instance),o) = deltaRho(s,t+4_pInt) +endforall +forall (s = 1:ns, c = 1_pInt:2_pInt) & + plasticState(p)%deltaState(iRhoD(s,c,instance),o) = deltaRho(s,c+8_pInt) +#else deltaState%p = 0.0_pReal forall (s = 1:ns, t = 1_pInt:4_pInt) deltaState%p(iRhoU(s,t,instance)) = deltaRho(s,t) @@ -3273,6 +2511,7 @@ forall (s = 1:ns, t = 1_pInt:4_pInt) endforall forall (s = 1:ns, c = 1_pInt:2_pInt) & deltaState%p(iRhoD(s,c,instance)) = deltaRho(s,c+8_pInt) +#endif #ifndef _OPENMP @@ -3287,635 +2526,15 @@ forall (s = 1:ns, c = 1_pInt:2_pInt) & end subroutine constitutive_nonlocal_deltaState -#endif - +!--------------------------------------------------------------------------------------------------- +!> @brief calculates the rate of change of microstructure +!--------------------------------------------------------------------------------------------------- #ifdef NEWSTATE -!-------------------------------------------------------------------------------------------------- -!> @brief calculates the rate of change of microstructure -!-------------------------------------------------------------------------------------------------- -function constitutive_nonlocal_dotState(Tstar_v, Fe, Fp, Temperature, mapping, timestep, subfrac, ipc,ip,el) - -use prec, only: DAMASK_NaN -use numerics, only: numerics_integrationMode, & - numerics_timeSyncing -use IO, only: IO_error -use debug, only: debug_level, & - debug_constitutive, & - debug_levelBasic, & - debug_levelExtensive, & - debug_levelSelective, & - debug_g, & - debug_i, & - debug_e -use math, only: math_norm3, & - math_mul6x6, & - math_mul3x3, & - math_mul33x3, & - math_mul33x33, & - math_inv33, & - math_det33, & - math_transpose33, & - pi -use mesh, only: mesh_NcpElems, & - mesh_maxNips, & - mesh_element, & - mesh_ipNeighborhood, & - mesh_ipVolume, & - mesh_ipArea, & - mesh_ipAreaNormal, & - FE_NipNeighbors, & - FE_geomtype, & - FE_celltype -use material, only: homogenization_maxNgrains, & - material_phase, & - phase_plasticityInstance, & - phase_localPlasticity, & - phase_plasticity ,& - PLASTICITY_NONLOCAL_ID, & - plasticState -use lattice, only: lattice_Sslip_v, & - lattice_sd, & - lattice_st ,& - lattice_mu, & - lattice_nu, & - lattice_structure, & - LATTICE_bcc_ID, & - LATTICE_fcc_ID - -implicit none - -!*** input variables -integer(pInt), intent(in) :: ipc, & !< current grain number - ip, & !< current integration point - el !< current element number -real(pReal), intent(in) :: Temperature, & !< temperature - timestep !< substepped crystallite time increment -real(pReal), dimension(6), intent(in) :: Tstar_v !< current 2nd Piola-Kirchhoff stress in Mandel notation -real(pReal), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & - subfrac !< fraction of timestep at the beginning of the substepped crystallite time increment -real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & - Fe, & !< elastic deformation gradient - Fp !< plastic deformation gradient -integer(pInt),dimension(2,homogenization_maxngrains,mesh_maxNips,mesh_ncpelems), intent(in) :: & - mapping -real(pReal), allocatable, dimension(:) :: & - constitutive_nonlocal_dotState !< evolution of state variables / microstructure - -!*** local variables -integer(pInt) :: phase, & - instance, & !< current instance of this plasticity - neighbor_instance, & !< instance of my neighbor's plasticity - ns, & !< short notation for the total number of active slip systems - c, & !< character of dislocation - n, & !< index of my current neighbor - neighbor_el, & !< element number of my neighbor - neighbor_ip, & !< integration point of my neighbor - neighbor_n, & !< neighbor index pointing to me when looking from my neighbor - opposite_neighbor, & !< index of my opposite neighbor - opposite_ip, & !< ip of my opposite neighbor - opposite_el, & !< element index of my opposite neighbor - opposite_n, & !< neighbor index pointing to me when looking from my opposite neighbor - t, & !< type of dislocation - topp, & !< type of dislocation with opposite sign to t - s, & !< index of my current slip system - sLattice !< index of my current slip system according to lattice order -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),10) :: & - rhoDot, & !< density evolution - rhoDotMultiplication, & !< density evolution by multiplication - rhoDotFlux, & !< density evolution by flux - rhoDotSingle2DipoleGlide, & !< density evolution by dipole formation (by glide) - rhoDotAthermalAnnihilation, & !< density evolution by athermal annihilation - rhoDotThermalAnnihilation !< density evolution by thermal annihilation -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),8) :: & - rhoSgl, & !< current single dislocation densities (positive/negative screw and edge without dipoles) - rhoSglOriginal, & - neighbor_rhoSgl, & !< current single dislocation densities of neighboring ip (positive/negative screw and edge without dipoles) - rhoSgl0, & !< single dislocation densities at start of cryst inc (positive/negative screw and edge without dipoles) - my_rhoSgl !< single dislocation densities of central ip (positive/negative screw and edge without dipoles) -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),4) :: & - v, & !< current dislocation glide velocity - v0, & !< dislocation glide velocity at start of cryst inc - my_v, & !< dislocation glide velocity of central ip - neighbor_v, & !< dislocation glide velocity of enighboring ip - gdot !< shear rates -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & - rhoForest, & !< forest dislocation density - tauThreshold, & !< threshold shear stress - tau, & !< current resolved shear stress - tauBack, & !< current back stress from pileups on same slip system - vClimb, & !< climb velocity of edge dipoles - nSources -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),2) :: & - rhoDip, & !< current dipole dislocation densities (screw and edge dipoles) - rhoDipOriginal, & - dLower, & !< minimum stable dipole distance for edges and screws - dUpper !< current maximum stable dipole distance for edges and screws -real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),4) :: & - m !< direction of dislocation motion -real(pReal), dimension(3,3) :: my_F, & !< my total deformation gradient - neighbor_F, & !< total deformation gradient of my neighbor - my_Fe, & !< my elastic deformation gradient - neighbor_Fe, & !< elastic deformation gradient of my neighbor - Favg !< average total deformation gradient of me and my neighbor -real(pReal), dimension(3) :: normal_neighbor2me, & !< interface normal pointing from my neighbor to me in neighbor's lattice configuration - normal_neighbor2me_defConf, & !< interface normal pointing from my neighbor to me in shared deformed configuration - normal_me2neighbor, & !< interface normal pointing from me to my neighbor in my lattice configuration - normal_me2neighbor_defConf !< interface normal pointing from me to my neighbor in shared deformed configuration -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 - selfDiffusion, & !< self diffusion - rnd, & - meshlength -logical considerEnteringFlux, & - considerLeavingFlux - -#ifndef _OPENMP - if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt & - .and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)& - .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then - write(6,*) - write(6,'(a,i8,1x,i2,1x,i1)') '<< CONST >> nonlocal_dotState at el ip ipc ',el,ip,ipc - write(6,*) - endif -#endif - -phase = material_phase(ipc,ip,el) -instance = phase_plasticityInstance(phase) -ns = totalNslip(instance) - -tau = 0.0_pReal -gdot = 0.0_pReal - -allocate(constitutive_nonlocal_dotState(plasticState(mapping(2,ipc,ip,el))%sizeDotState),source=0.0_pReal) -!*** shortcut to state variables - - -forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) - rhoSgl(s,t) = max(plasticState(mapping(2,ipc,ip,el))%state(iRhoU(s,t,instance),mapping(1,ipc,ip,el)), 0.0_pReal) ! ensure positive single mobile densities - rhoSgl(s,t+4_pInt) = plasticState(mapping(2,ipc,ip,el))%state(iRhoB(s,t,instance),mapping(1,ipc,ip,el)) - v(s,t) = plasticState(mapping(2,ipc,ip,el))%state(iV(s,t,instance),mapping(1,ipc,ip,el)) -endforall -forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) - rhoDip(s,c) = max(plasticState(mapping(2,ipc,ip,el))%state(iRhoD(s,c,instance),mapping(1,ipc,ip,el)), 0.0_pReal) ! ensure positive dipole densities - -endforall -rhoForest = plasticState(mapping(2,ipc,ip,el))%state(iRhoF(1:ns,instance),mapping(1,ipc,ip,el)) -tauThreshold = plasticState(mapping(2,ipc,ip,el))%state(iTauF(1:ns,instance),mapping(1,ipc,ip,el)) -tauBack = plasticState(mapping(2,ipc,ip,el))%state(iTauB(1:ns,instance),mapping(1,ipc,ip,el)) - -rhoSglOriginal = rhoSgl -rhoDipOriginal = rhoDip -where (abs(rhoSgl) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(instance) & - .or. abs(rhoSgl) < significantRho(instance)) & - rhoSgl = 0.0_pReal -where (abs(rhoDip) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(instance) & - .or. abs(rhoDip) < significantRho(instance)) & - rhoDip = 0.0_pReal - -if (numerics_timeSyncing) then - forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) - - rhoSgl0(s,t) = max(plasticState(mapping(2,ipc,ip,el))%state0(iRhoU(s,t,instance), & - mapping(1,ipc,ip,el)), 0.0_pReal) - rhoSgl0(s,t+4_pInt) = plasticState(mapping(2,ipc,ip,el))%state0(iRhoB(s,t,instance),mapping(1,ipc,ip,el)) - v0(s,t) = plasticState(mapping(2,ipc,ip,el))%state0(iV(s,t,instance),mapping(1,ipc,ip,el)) - endforall - where (abs(rhoSgl0) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(instance) & - .or. abs(rhoSgl0) < significantRho(instance)) & - rhoSgl0 = 0.0_pReal -endif - - - -!*** sanity check for timestep - -if (timestep <= 0.0_pReal) then ! if illegal timestep... - constitutive_nonlocal_dotState = 0.0_pReal ! ...return without doing anything (-> zero dotState) - return -endif - - - -!**************************************************************************** -!*** Calculate shear rate - -forall (t = 1_pInt:4_pInt) & - gdot(1_pInt:ns,t) = rhoSgl(1_pInt:ns,t) * burgers(1:ns,instance) * v(1:ns,t) - -#ifndef _OPENMP - if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt & - .and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)& - .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then - write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> rho / 1/m^2', rhoSgl, rhoDip - write(6,'(a,/,4(12x,12(e12.5,1x),/))') '<< CONST >> gdot / 1/s',gdot - endif -#endif - - - -!**************************************************************************** -!*** calculate limits for stable dipole height - -do s = 1_pInt,ns ! loop over slip systems - sLattice = slipSystemLattice(s,instance) - tau(s) = math_mul6x6(Tstar_v, lattice_Sslip_v(1:6,1,sLattice,phase)) + tauBack(s) - if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal -enddo - -dLower = minDipoleHeight(1:ns,1:2,instance) -dUpper(1:ns,1) = lattice_mu(phase) * burgers(1:ns,instance) & - / (8.0_pReal * pi * (1.0_pReal - lattice_nu(phase)) * abs(tau)) -dUpper(1:ns,2) = lattice_mu(phase) * burgers(1:ns,instance) & - / (4.0_pReal * pi * abs(tau)) -forall (c = 1_pInt:2_pInt) & - dUpper(1:ns,c) = min(1.0_pReal / sqrt(rhoSgl(1:ns,2*c-1) + rhoSgl(1:ns,2*c) & - + abs(rhoSgl(1:ns,2*c+3)) + abs(rhoSgl(1:ns,2*c+4)) + rhoDip(1:ns,c)), & - dUpper(1:ns,c)) -dUpper = max(dUpper,dLower) - - - -!**************************************************************************** -!*** calculate dislocation multiplication - -rhoDotMultiplication = 0.0_pReal -if (lattice_structure(phase) == LATTICE_bcc_ID) then ! BCC - 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,instance) & ! assuming double-cross-slip of screws to be decisive for multiplication - * sqrt(rhoForest(s)) / lambda0(s,instance) ! & ! 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 - rhoDotMultiplication(s,3:4) = sum(abs(gdot(s,3:4))) / burgers(s,instance) & ! assuming double-cross-slip of screws to be decisive for multiplication - * sqrt(rhoForest(s)) / lambda0(s,instance) ! & ! 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 - endforall - -else ! ALL OTHER STRUCTURES - if (probabilisticMultiplication(instance)) then - meshlength = mesh_ipVolume(ip,el)**0.333_pReal - where(sum(rhoSgl(1:ns,1:4),2) > 0.0_pReal) - nSources = (sum(rhoSgl(1:ns,1:2),2) * fEdgeMultiplication(instance) + sum(rhoSgl(1:ns,3:4),2)) & - / sum(rhoSgl(1:ns,1:4),2) * meshlength / lambda0(1:ns,instance) * sqrt(rhoForest(1:ns)) - elsewhere - nSources = meshlength / lambda0(1:ns,instance) * sqrt(rhoForest(1:ns)) - endwhere - do s = 1_pInt,ns - if (nSources(s) < 1.0_pReal) then - if (sourceProbability(s,ipc,ip,el) > 1.0_pReal) then - call random_number(rnd) - sourceProbability(s,ipc,ip,el) = rnd - !$OMP FLUSH(sourceProbability) - endif - if (sourceProbability(s,ipc,ip,el) > 1.0_pReal - nSources(s)) then - rhoDotMultiplication(s,1:4) = sum(rhoSglOriginal(s,1:4) * abs(v(s,1:4))) / meshlength - endif - else - sourceProbability(s,ipc,ip,el) = 2.0_pReal - rhoDotMultiplication(s,1:4) = & - (sum(abs(gdot(s,1:2))) * fEdgeMultiplication(instance) + sum(abs(gdot(s,3:4)))) & - / burgers(s,instance) * sqrt(rhoForest(s)) / lambda0(s,instance) - endif - enddo -#ifndef _OPENMP - if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt & - .and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)& - .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then - write(6,'(a,/,4(12x,12(f12.5,1x),/))') '<< CONST >> sources', nSources - write(6,*) - endif -#endif - else - rhoDotMultiplication(1:ns,1:4) = spread( & - (sum(abs(gdot(1:ns,1:2)),2) * fEdgeMultiplication(instance) + sum(abs(gdot(1:ns,3:4)),2)) & - * sqrt(rhoForest(1:ns)) / lambda0(1:ns,instance) / burgers(1:ns,instance), 2, 4) - endif -endif - - - -!**************************************************************************** -!*** calculate dislocation fluxes (only for nonlocal plasticity) - -rhoDotFlux = 0.0_pReal - -if (.not. phase_localPlasticity(material_phase(ipc,ip,el))) then ! only for nonlocal plasticity - - - !*** check CFL (Courant-Friedrichs-Lewy) condition for flux - - if (any( abs(gdot) > 0.0_pReal & ! any active slip system ... - .and. CFLfactor(instance) * abs(v) * timestep & - > mesh_ipVolume(ip,el) / maxval(mesh_ipArea(:,ip,el)))) then ! ...with velocity above critical value (we use the reference volume and area for simplicity here) -#ifndef _OPENMP - if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt) then - write(6,'(a,i5,a,i2)') '<< CONST >> CFL condition not fullfilled at el ',el,' ip ',ip - write(6,'(a,e10.3,a,e10.3)') '<< CONST >> velocity is at ', & - maxval(abs(v), abs(gdot) > 0.0_pReal & - .and. CFLfactor(instance) * abs(v) * timestep & - > mesh_ipVolume(ip,el) / maxval(mesh_ipArea(:,ip,el))), & - ' at a timestep of ',timestep - write(6,'(a)') '<< CONST >> enforcing cutback !!!' - endif -#endif - constitutive_nonlocal_dotState = DAMASK_NaN ! -> return NaN and, hence, enforce cutback - return - endif - - - !*** be aware of the definition of lattice_st = lattice_sd x lattice_sn !!! - !*** opposite sign to our p vector in the (s,p,n) triplet !!! - - m(1:3,1:ns,1) = lattice_sd(1:3, slipSystemLattice(1:ns,instance), phase) - m(1:3,1:ns,2) = -lattice_sd(1:3, slipSystemLattice(1:ns,instance), phase) - m(1:3,1:ns,3) = -lattice_st(1:3, slipSystemLattice(1:ns,instance), phase) - m(1:3,1:ns,4) = lattice_st(1:3, slipSystemLattice(1:ns,instance), phase) - - my_Fe = Fe(1:3,1:3,ipc,ip,el) - my_F = math_mul33x33(my_Fe, Fp(1:3,1:3,ipc,ip,el)) - - do n = 1_pInt,FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,el)))) ! loop through my neighbors - neighbor_el = mesh_ipNeighborhood(1,n,ip,el) - neighbor_ip = mesh_ipNeighborhood(2,n,ip,el) - neighbor_n = mesh_ipNeighborhood(3,n,ip,el) - - opposite_neighbor = n + mod(n,2_pInt) - mod(n+1_pInt,2_pInt) - opposite_el = mesh_ipNeighborhood(1,opposite_neighbor,ip,el) - opposite_ip = mesh_ipNeighborhood(2,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 - neighbor_instance = phase_plasticityInstance(material_phase(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)) - Favg = 0.5_pReal * (my_F + neighbor_F) - else ! if no neighbor, take my value as average - Favg = my_F - endif - - - !* FLUX FROM MY NEIGHBOR TO ME - !* This is only considered, if I have a neighbor of nonlocal plasticity (also nonlocal constitutive law with local properties) 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. - neighbor_v = 0.0_pReal ! needed for check of sign change in flux density below - neighbor_rhoSgl = 0.0_pReal - if (neighbor_n > 0_pInt) then - if (phase_plasticity(material_phase(1,neighbor_ip,neighbor_el)) == PLASTICITY_NONLOCAL_ID & - .and. any(compatibility(:,:,:,n,ip,el) > 0.0_pReal)) & - considerEnteringFlux = .true. - endif - - 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 - forall (s = 1:ns, t = 1_pInt:4_pInt) - neighbor_v(s,t) = plasticState(mapping(2,ipc,neighbor_ip,neighbor_el))% & - state0(iV(s,t,neighbor_instance),mapping(1,ipc,ip,el)) - neighbor_rhoSgl(s,t) = max(plasticState(mapping(2,ipc,neighbor_ip,neighbor_el))% & - state0(iRhoU(s,t,neighbor_instance),mapping(1,ipc,ip,el)), 0.0_pReal) - endforall - else - forall (s = 1:ns, t = 1_pInt:4_pInt) - neighbor_v(s,t) = plasticState(mapping(2,ipc,neighbor_ip,neighbor_el))% & - state(iV(s,t,neighbor_instance),mapping(1,ipc,ip,el)) - neighbor_rhoSgl(s,t) = max(plasticState(mapping(2,ipc,neighbor_ip,neighbor_el))% & - state(iRhoU(s,t,neighbor_instance),mapping(1,ipc,ip,el)), 0.0_pReal) - endforall - endif - where (neighbor_rhoSgl * mesh_ipVolume(neighbor_ip,neighbor_el) ** 0.667_pReal < significantN(instance) & - .or. neighbor_rhoSgl < significantRho(instance)) & - neighbor_rhoSgl = 0.0_pReal - normal_neighbor2me_defConf = math_det33(Favg) * math_mul33x3(math_inv33(transpose(Favg)), & - mesh_ipAreaNormal(1:3,neighbor_n,neighbor_ip,neighbor_el)) ! calculate the normal of the interface in (average) deformed configuration (now pointing from my neighbor to me!!!) - normal_neighbor2me = math_mul33x3(transpose(neighbor_Fe), normal_neighbor2me_defConf) & - / math_det33(neighbor_Fe) ! interface normal in the lattice configuration of my neighbor - area = mesh_ipArea(neighbor_n,neighbor_ip,neighbor_el) * math_norm3(normal_neighbor2me) - normal_neighbor2me = normal_neighbor2me / math_norm3(normal_neighbor2me) ! normalize the surface normal to unit length - do s = 1_pInt,ns - do t = 1_pInt,4_pInt - c = (t + 1_pInt) / 2 - topp = t + mod(t,2_pInt) - mod(t+1_pInt,2_pInt) - if (neighbor_v(s,t) * math_mul3x3(m(1:3,s,t), normal_neighbor2me) > 0.0_pReal & ! flux from my neighbor to me == entering flux for me - .and. v(s,t) * neighbor_v(s,t) >= 0.0_pReal ) then ! ... only if no sign change in flux density - lineLength = neighbor_rhoSgl(s,t) * neighbor_v(s,t) & - * math_mul3x3(m(1:3,s,t), normal_neighbor2me) * area ! positive line length that wants to enter through this interface - where (compatibility(c,1_pInt:ns,s,n,ip,el) > 0.0_pReal) & ! positive compatibility... - rhoDotFlux(1_pInt:ns,t) = rhoDotFlux(1_pInt:ns,t) & - + lineLength / mesh_ipVolume(ip,el) & ! ... transferring to equally signed mobile dislocation type - * compatibility(c,1_pInt:ns,s,n,ip,el) ** 2.0_pReal - where (compatibility(c,1_pInt:ns,s,n,ip,el) < 0.0_pReal) & ! ..negative compatibility... - rhoDotFlux(1_pInt:ns,topp) = rhoDotFlux(1_pInt:ns,topp) & - + lineLength / mesh_ipVolume(ip,el) & ! ... transferring to opposite signed mobile dislocation type - * compatibility(c,1_pInt:ns,s,n,ip,el) ** 2.0_pReal - endif - enddo - enddo - endif - - - !* FLUX FROM ME TO MY NEIGHBOR - !* This is not considered, if my opposite neighbor has a different constitutive law than nonlocal (still considered for nonlocal law with lcal properties). - !* 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. - - considerLeavingFlux = .true. - if (opposite_n > 0_pInt) then - if (phase_plasticity(material_phase(1,opposite_ip,opposite_el)) /= PLASTICITY_NONLOCAL_ID) & - considerLeavingFlux = .false. - endif - - if (considerLeavingFlux) then - - !* timeSyncing mode: If the central ip has zero subfraction, always use "state0". This is needed in case of - !* a synchronization step for the central ip, because then "state" contains the values at the end of the - !* previously converged full time step. Also, if either me or my neighbor has zero subfraction, we have to - !* use "state0" to make sure that fluxes on both sides of the (potential) timestep are equal. - my_rhoSgl = rhoSgl - my_v = v - if(numerics_timeSyncing) then - if (subfrac(ipc,ip,el) == 0.0_pReal) then - my_rhoSgl = rhoSgl0 - my_v = v0 - elseif (neighbor_n > 0_pInt) then - if (subfrac(ipc,neighbor_ip,neighbor_el) == 0.0_pReal) then - my_rhoSgl = rhoSgl0 - my_v = v0 - endif - endif - endif - - normal_me2neighbor_defConf = math_det33(Favg) & - * math_mul33x3(math_inv33(math_transpose33(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(math_transpose33(my_Fe), normal_me2neighbor_defConf) & - / math_det33(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_pInt,ns - do t = 1_pInt,4_pInt - c = (t + 1_pInt) / 2_pInt - if (my_v(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) - if (my_v(s,t) * neighbor_v(s,t) >= 0.0_pReal) then ! no sign change in flux density - transmissivity = sum(compatibility(c,1_pInt:ns,s,n,ip,el)**2.0_pReal) ! overall transmissivity from this slip system to my neighbor - else ! sign change in flux density means sign change in stress which does not allow for dislocations to arive at the neighbor - transmissivity = 0.0_pReal - endif - lineLength = my_rhoSgl(s,t) * my_v(s,t) & - * math_mul3x3(m(1:3,s,t), normal_me2neighbor) * area ! positive line length of mobiles that wants to leave through this interface - rhoDotFlux(s,t) = rhoDotFlux(s,t) - lineLength / mesh_ipVolume(ip,el) ! subtract dislocation flux from current type - rhoDotFlux(s,t+4_pInt) = rhoDotFlux(s,t+4_pInt) & - + lineLength / mesh_ipVolume(ip,el) * (1.0_pReal - transmissivity) & - * sign(1.0_pReal, my_v(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 - - enddo ! neighbor loop -endif - - - -!**************************************************************************** -!*** calculate dipole formation and annihilation - -!*** formation by glide - -do c = 1_pInt,2_pInt - rhoDotSingle2DipoleGlide(1:ns,2*c-1) = -2.0_pReal * dUpper(1:ns,c) / burgers(1:ns,instance) & - * (rhoSgl(1:ns,2*c-1) * abs(gdot(1:ns,2*c)) & ! negative mobile --> positive mobile - + rhoSgl(1:ns,2*c) * abs(gdot(1:ns,2*c-1)) & ! positive mobile --> negative mobile - + abs(rhoSgl(1:ns,2*c+4)) * abs(gdot(1:ns,2*c-1))) ! positive mobile --> negative immobile - - rhoDotSingle2DipoleGlide(1:ns,2*c) = -2.0_pReal * dUpper(1:ns,c) / burgers(1:ns,instance) & - * (rhoSgl(1:ns,2*c-1) * abs(gdot(1:ns,2*c)) & ! negative mobile --> positive mobile - + rhoSgl(1:ns,2*c) * abs(gdot(1:ns,2*c-1)) & ! positive mobile --> negative mobile - + abs(rhoSgl(1:ns,2*c+3)) * abs(gdot(1:ns,2*c))) ! negative mobile --> positive immobile - - rhoDotSingle2DipoleGlide(1:ns,2*c+3) = -2.0_pReal * dUpper(1:ns,c) / burgers(1:ns,instance) & - * rhoSgl(1:ns,2*c+3) * abs(gdot(1:ns,2*c)) ! negative mobile --> positive immobile - - rhoDotSingle2DipoleGlide(1:ns,2*c+4) = -2.0_pReal * dUpper(1:ns,c) / burgers(1:ns,instance) & - * rhoSgl(1:ns,2*c+4) * abs(gdot(1:ns,2*c-1)) ! positive mobile --> negative immobile - - rhoDotSingle2DipoleGlide(1:ns,c+8) = - rhoDotSingle2DipoleGlide(1:ns,2*c-1) & - - rhoDotSingle2DipoleGlide(1:ns,2*c) & - + abs(rhoDotSingle2DipoleGlide(1:ns,2*c+3)) & - + abs(rhoDotSingle2DipoleGlide(1:ns,2*c+4)) -enddo - - -!*** athermal annihilation - -rhoDotAthermalAnnihilation = 0.0_pReal - -forall (c=1_pInt:2_pInt) & - rhoDotAthermalAnnihilation(1:ns,c+8_pInt) = -2.0_pReal * dLower(1:ns,c) / burgers(1:ns,instance) & - * ( 2.0_pReal * (rhoSgl(1:ns,2*c-1) * abs(gdot(1:ns,2*c)) + rhoSgl(1:ns,2*c) * abs(gdot(1:ns,2*c-1))) & ! was single hitting single - + 2.0_pReal * (abs(rhoSgl(1:ns,2*c+3)) * abs(gdot(1:ns,2*c)) + abs(rhoSgl(1:ns,2*c+4)) * abs(gdot(1:ns,2*c-1))) & ! was single hitting immobile single or was immobile single hit by single - + rhoDip(1:ns,c) * (abs(gdot(1:ns,2*c-1)) + abs(gdot(1:ns,2*c)))) ! single knocks dipole constituent -! annihilated screw dipoles leave edge jogs behind on the colinear system -if (lattice_structure(phase) == LATTICE_fcc_ID) then ! only fcc - forall (s = 1:ns, colinearSystem(s,instance) > 0_pInt) & - rhoDotAthermalAnnihilation(colinearSystem(s,instance),1:2) = - rhoDotAthermalAnnihilation(s,10) & - * 0.25_pReal * sqrt(rhoForest(s)) * (dUpper(s,2) + dLower(s,2)) * edgeJogFactor(instance) -endif - - -!*** thermally activated annihilation of edge dipoles by climb - -rhoDotThermalAnnihilation = 0.0_pReal -selfDiffusion = Dsd0(instance) * exp(-selfDiffusionEnergy(instance) / (KB * Temperature)) -vClimb = atomicVolume(instance) * selfDiffusion / ( KB * Temperature ) & - * lattice_mu(phase) / ( 2.0_pReal * PI * (1.0_pReal-lattice_nu(phase)) ) & - * 2.0_pReal / ( dUpper(1:ns,1) + dLower(1:ns,1) ) -forall (s = 1_pInt:ns, dUpper(s,1) > dLower(s,1)) & - rhoDotThermalAnnihilation(s,9) = max(- 4.0_pReal * rhoDip(s,1) * vClimb(s) / (dUpper(s,1) - dLower(s,1)), & - - rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) & - - rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have - - - -!**************************************************************************** -!*** assign the rates of dislocation densities to my dotState -!*** if evolution rates lead to negative densities, a cutback is enforced - -rhoDot = 0.0_pReal -rhoDot = rhoDotFlux & - + rhoDotMultiplication & - + rhoDotSingle2DipoleGlide & - + rhoDotAthermalAnnihilation & - + rhoDotThermalAnnihilation - -if (numerics_integrationMode == 1_pInt) then ! save rates for output if in central integration mode - rhoDotFluxOutput(1:ns,1:8,ipc,ip,el) = rhoDotFlux(1:ns,1:8) - rhoDotMultiplicationOutput(1:ns,1:2,ipc,ip,el) = rhoDotMultiplication(1:ns,[1,3]) - rhoDotSingle2DipoleGlideOutput(1:ns,1:2,ipc,ip,el) = rhoDotSingle2DipoleGlide(1:ns,9:10) - rhoDotAthermalAnnihilationOutput(1:ns,1:2,ipc,ip,el) = rhoDotAthermalAnnihilation(1:ns,9:10) - rhoDotThermalAnnihilationOutput(1:ns,1:2,ipc,ip,el) = rhoDotThermalAnnihilation(1:ns,9:10) - rhoDotEdgeJogsOutput(1:ns,ipc,ip,el) = 2.0_pReal * rhoDotThermalAnnihilation(1:ns,1) -endif - - -#ifndef _OPENMP - if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt & - .and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)& - .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then - write(6,'(a,/,4(12x,12(e12.5,1x),/))') '<< CONST >> dislocation multiplication', & - rhoDotMultiplication(1:ns,1:4) * timestep - write(6,'(a,/,8(12x,12(e12.5,1x),/))') '<< CONST >> dislocation flux', & - rhoDotFlux(1:ns,1:8) * timestep - write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> dipole formation by glide', & - rhoDotSingle2DipoleGlide * timestep - write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> athermal dipole annihilation', & - rhoDotAthermalAnnihilation * timestep - write(6,'(a,/,2(12x,12(e12.5,1x),/))') '<< CONST >> thermally activated dipole annihilation', & - rhoDotThermalAnnihilation(1:ns,9:10) * timestep - write(6,'(a,/,10(12x,12(e12.5,1x),/))') '<< CONST >> total density change', & - rhoDot * timestep - write(6,'(a,/,10(12x,12(f12.5,1x),/))') '<< CONST >> relative density change', & - rhoDot(1:ns,1:8) * timestep / (abs(rhoSglOriginal)+1.0e-10), & - rhoDot(1:ns,9:10) * timestep / (rhoDipOriginal+1.0e-10) - write(6,*) - endif -#endif - - -if ( any(rhoSglOriginal(1:ns,1:4) + rhoDot(1:ns,1:4) * timestep < -aTolRho(instance)) & - .or. any(rhoDipOriginal(1:ns,1:2) + rhoDot(1:ns,9:10) * timestep < -aTolRho(instance))) then -#ifndef _OPENMP - if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt) then - write(6,'(a,i5,a,i2)') '<< CONST >> evolution rate leads to negative density at el ',el,' ip ',ip - write(6,'(a)') '<< CONST >> enforcing cutback !!!' - endif -#endif - constitutive_nonlocal_dotState = DAMASK_NaN - return -else - forall (s = 1:ns, t = 1_pInt:4_pInt) - constitutive_nonlocal_dotState(iRhoU(s,t,instance)) = rhoDot(s,t) - constitutive_nonlocal_dotState(iRhoB(s,t,instance)) = rhoDot(s,t+4_pInt) - endforall - forall (s = 1:ns, c = 1_pInt:2_pInt) & - constitutive_nonlocal_dotState(iRhoD(s,c,instance)) = rhoDot(s,c+8_pInt) - forall (s = 1:ns) & - constitutive_nonlocal_dotState(iGamma(s,instance)) = sum(gdot(s,1:4)) -endif - -end function constitutive_nonlocal_dotState - - - +function constitutive_nonlocal_dotState(Tstar_v, Fe, Fp, Temperature, timestep,subfrac, ipc,ip,el) #else -!-------------------------------------------------------------------------------------------------- -!> @brief calculates the rate of change of microstructure -!-------------------------------------------------------------------------------------------------- -function constitutive_nonlocal_dotState(Tstar_v, Fe, Fp, Temperature, state, state0, timestep, subfrac, ipc,ip,el) - +function constitutive_nonlocal_dotState(Tstar_v, Fe, Fp, Temperature, state, state0, timestep, & + subfrac, ipc,ip,el) +#endif use prec, only: DAMASK_NaN use numerics, only: numerics_integrationMode, & numerics_timeSyncing @@ -3951,6 +2570,10 @@ use material, only: homogenization_maxNgrains, & material_phase, & phase_plasticityInstance, & phase_localPlasticity, & +#ifdef NEWSTATE + plasticState, & + mappingConstitutive, & +#endif phase_plasticity ,& PLASTICITY_NONLOCAL_ID use lattice, only: lattice_Sslip_v, & @@ -3965,93 +2588,112 @@ use lattice, only: lattice_Sslip_v, & implicit none !*** input variables -integer(pInt), intent(in) :: ipc, & !< current grain number - ip, & !< current integration point - el !< current element number -real(pReal), intent(in) :: Temperature, & !< temperature - timestep !< substepped crystallite time increment -real(pReal), dimension(6), intent(in) :: Tstar_v !< current 2nd Piola-Kirchhoff stress in Mandel notation +integer(pInt), intent(in) :: ipc, & !< current grain number + ip, & !< current integration point + el !< current element number +real(pReal), intent(in) :: Temperature, & !< temperature + timestep !< substepped crystallite time increment +real(pReal), dimension(6), intent(in) :: Tstar_v !< current 2nd Piola-Kirchhoff stress in Mandel notation real(pReal), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & - subfrac !< fraction of timestep at the beginning of the substepped crystallite time increment + subfrac !< fraction of timestep at the beginning of the substepped crystallite time increment real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & - Fe, & !< elastic deformation gradient - Fp !< plastic deformation gradient + Fe, & !< elastic deformation gradient + Fp !< plastic deformation gradient +#ifndef NEWSTATE type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & - state, & !< current microstructural state - state0 !< microstructural state at beginning of crystallite increment + state, & !< current microstructural state + state0 !< microstructural state at beginning of crystallite increment +#endif !*** output variables +#ifdef NEWSTATE +real(pReal), dimension(:), allocatable :: constitutive_nonlocal_dotState !< evolution of state variables / microstructure +#else real(pReal), dimension(constitutive_nonlocal_sizeDotState(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & - constitutive_nonlocal_dotState !< evolution of state variables / microstructure + constitutive_nonlocal_dotState !< evolution of state variables / microstructure +#endif !*** local variables integer(pInt) :: phase, & - instance, & !< current instance of this plasticity - neighbor_instance, & !< instance of my neighbor's plasticity - ns, & !< short notation for the total number of active slip systems - c, & !< character of dislocation - n, & !< index of my current neighbor - neighbor_el, & !< element number of my neighbor - neighbor_ip, & !< integration point of my neighbor - neighbor_n, & !< neighbor index pointing to me when looking from my neighbor - opposite_neighbor, & !< index of my opposite neighbor - opposite_ip, & !< ip of my opposite neighbor - opposite_el, & !< element index of my opposite neighbor - opposite_n, & !< neighbor index pointing to me when looking from my opposite neighbor - t, & !< type of dislocation - topp, & !< type of dislocation with opposite sign to t - s, & !< index of my current slip system - sLattice !< index of my current slip system according to lattice order + instance, & !< current instance of this plasticity + neighbor_instance, & !< instance of my neighbor's plasticity + ns, & !< short notation for the total number of active slip systems + c, & !< character of dislocation + n, & !< index of my current neighbor + neighbor_el, & !< element number of my neighbor + neighbor_ip, & !< integration point of my neighbor + neighbor_n, & !< neighbor index pointing to me when looking from my neighbor + opposite_neighbor, & !< index of my opposite neighbor + opposite_ip, & !< ip of my opposite neighbor + opposite_el, & !< element index of my opposite neighbor + opposite_n, & !< neighbor index pointing to me when looking from my opposite neighbor + t, & !< type of dislocation +#ifdef NEWSTATE + o,& !< offset shortcut + no,& !< neighbour offset shortcut + p,& !< phase shortcut + np,& !< neighbour phase shortcut +#endif + topp, & !< type of dislocation with opposite sign to t + s, & !< index of my current slip system + sLattice !< index of my current slip system according to lattice order real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),10) :: & - rhoDot, & !< density evolution - rhoDotMultiplication, & !< density evolution by multiplication - rhoDotFlux, & !< density evolution by flux - rhoDotSingle2DipoleGlide, & !< density evolution by dipole formation (by glide) - rhoDotAthermalAnnihilation, & !< density evolution by athermal annihilation - rhoDotThermalAnnihilation !< density evolution by thermal annihilation + rhoDot, & !< density evolution + rhoDotMultiplication, & !< density evolution by multiplication + rhoDotFlux, & !< density evolution by flux + rhoDotSingle2DipoleGlide, & !< density evolution by dipole formation (by glide) + rhoDotAthermalAnnihilation, & !< density evolution by athermal annihilation + rhoDotThermalAnnihilation !< density evolution by thermal annihilation real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,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) rhoSglOriginal, & - neighbor_rhoSgl, & !< current single dislocation densities of neighboring ip (positive/negative screw and edge without dipoles) - rhoSgl0, & !< single dislocation densities at start of cryst inc (positive/negative screw and edge without dipoles) - my_rhoSgl !< single dislocation densities of central ip (positive/negative screw and edge without dipoles) + neighbor_rhoSgl, & !< current single dislocation densities of neighboring ip (positive/negative screw and edge without dipoles) + rhoSgl0, & !< single dislocation densities at start of cryst inc (positive/negative screw and edge without dipoles) + my_rhoSgl !< single dislocation densities of central ip (positive/negative screw and edge without dipoles) real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),4) :: & - v, & !< current dislocation glide velocity - v0, & !< dislocation glide velocity at start of cryst inc - my_v, & !< dislocation glide velocity of central ip - neighbor_v, & !< dislocation glide velocity of enighboring ip - gdot !< shear rates + v, & !< current dislocation glide velocity + v0, & !< dislocation glide velocity at start of cryst inc + my_v, & !< dislocation glide velocity of central ip + neighbor_v, & !< dislocation glide velocity of enighboring ip + gdot !< shear rates real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: & - rhoForest, & !< forest dislocation density - tauThreshold, & !< threshold shear stress - tau, & !< current resolved shear stress - tauBack, & !< current back stress from pileups on same slip system - vClimb, & !< climb velocity of edge dipoles + rhoForest, & !< forest dislocation density + tauThreshold, & !< threshold shear stress + tau, & !< current resolved shear stress + tauBack, & !< current back stress from pileups on same slip system + vClimb, & !< climb velocity of edge dipoles nSources real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),2) :: & - rhoDip, & !< current dipole dislocation densities (screw and edge dipoles) + rhoDip, & !< current dipole dislocation densities (screw and edge dipoles) rhoDipOriginal, & - dLower, & !< minimum stable dipole distance for edges and screws - dUpper !< current maximum 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 real(pReal), dimension(3,totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),4) :: & - m !< direction of dislocation motion -real(pReal), dimension(3,3) :: my_F, & !< my total deformation gradient - neighbor_F, & !< total deformation gradient of my neighbor - my_Fe, & !< my elastic deformation gradient - neighbor_Fe, & !< elastic deformation gradient of my neighbor - Favg !< average total deformation gradient of me and my neighbor -real(pReal), dimension(3) :: normal_neighbor2me, & !< interface normal pointing from my neighbor to me in neighbor's lattice configuration - normal_neighbor2me_defConf, & !< interface normal pointing from my neighbor to me in shared deformed configuration - normal_me2neighbor, & !< interface normal pointing from me to my neighbor in my lattice configuration - normal_me2neighbor_defConf !< interface normal pointing from me to my neighbor in shared deformed configuration -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 - selfDiffusion, & !< self diffusion + m !< direction of dislocation motion +real(pReal), dimension(3,3) :: my_F, & !< my total deformation gradient + neighbor_F, & !< total deformation gradient of my neighbor + my_Fe, & !< my elastic deformation gradient + neighbor_Fe, & !< elastic deformation gradient of my neighbor + Favg !< average total deformation gradient of me and my neighbor +real(pReal), dimension(3) :: normal_neighbor2me, & !< interface normal pointing from my neighbor to me in neighbor's lattice configuration + normal_neighbor2me_defConf, & !< interface normal pointing from my neighbor to me in shared deformed configuration + normal_me2neighbor, & !< interface normal pointing from me to my neighbor in my lattice configuration + normal_me2neighbor_defConf !< interface normal pointing from me to my neighbor in shared deformed configuration +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 + selfDiffusion, & !< self diffusion rnd, & meshlength logical considerEnteringFlux, & considerLeavingFlux + +#ifdef NEWSTATE + p = mappingConstitutive(2,1,ip,el) + o = mappingConstitutive(1,1,ip,el) +!allocate function var +allocate(constitutive_nonlocal_dotState(plasticState(p)%sizeDotState), source = 0.0_pReal) +#endif #ifndef _OPENMP if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt & @@ -4073,18 +2715,31 @@ gdot = 0.0_pReal !*** shortcut to state variables - +#ifdef NEWSTATE 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(plasticState(p)%state((iRhoU(s,t,instance)),o), 0.0_pReal) ! ensure positive single mobile densities + rhoSgl(s,t+4_pInt) = plasticState(p)%state((iRhoB(s,t,instance)),o) + v(s,t) = plasticState(p)%state((iV(s,t,instance)),o) +endforall +forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) + rhoDip(s,c) = max(plasticState(p)%state((iRhoD(s,c,instance)),o), 0.0_pReal) ! ensure positive dipole densities +endforall +rhoForest = plasticState(p)%state((iRhoF(1:ns,instance)),o) +tauThreshold = plasticState(p)%state((iTauF(1:ns,instance)),o) +tauBack = plasticState(p)%state((iTauB(1:ns,instance)),o) +#else +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+4_pInt) = state(ipc,ip,el)%p(iRhoB(s,t,instance)) v(s,t) = state(ipc,ip,el)%p(iV(s,t,instance)) endforall forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) - rhoDip(s,c) = max(state(ipc,ip,el)%p(iRhoD(s,c,instance)), 0.0_pReal) ! ensure positive dipole densities + rhoDip(s,c) = max(state(ipc,ip,el)%p(iRhoD(s,c,instance)), 0.0_pReal) ! ensure positive dipole densities endforall rhoForest = state(ipc,ip,el)%p(iRhoF(1:ns,instance)) tauThreshold = state(ipc,ip,el)%p(iTauF(1:ns,instance)) tauBack = state(ipc,ip,el)%p(iTauB(1:ns,instance)) +#endif rhoSglOriginal = rhoSgl rhoDipOriginal = rhoDip @@ -4096,11 +2751,19 @@ where (abs(rhoDip) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(instance rhoDip = 0.0_pReal if (numerics_timeSyncing) then +#ifdef NEWSTATE + forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) + rhoSgl0(s,t) = max(plasticState(p)%state0(iRhoU(s,t,instance),o), 0.0_pReal) + rhoSgl0(s,t+4_pInt) = plasticState(p)%state0(iRhoB(s,t,instance),o) + v0(s,t) = plasticState(p)%state0(iV(s,t,instance),o) + endforall +#else forall (s = 1_pInt:ns, t = 1_pInt:4_pInt) rhoSgl0(s,t) = max(state0(ipc,ip,el)%p(iRhoU(s,t,instance)), 0.0_pReal) rhoSgl0(s,t+4_pInt) = state0(ipc,ip,el)%p(iRhoB(s,t,instance)) v0(s,t) = state0(ipc,ip,el)%p(iV(s,t,instance)) endforall +#endif where (abs(rhoSgl0) * mesh_ipVolume(ip,el) ** 0.667_pReal < significantN(instance) & .or. abs(rhoSgl0) < significantRho(instance)) & rhoSgl0 = 0.0_pReal @@ -4110,8 +2773,8 @@ endif !*** sanity check for timestep -if (timestep <= 0.0_pReal) then ! if illegal timestep... - constitutive_nonlocal_dotState = 0.0_pReal ! ...return without doing anything (-> zero dotState) +if (timestep <= 0.0_pReal) then ! if illegal timestep... + constitutive_nonlocal_dotState = 0.0_pReal ! ...return without doing anything (-> zero dotState) return endif @@ -4175,7 +2838,7 @@ else meshlength = mesh_ipVolume(ip,el)**0.333_pReal where(sum(rhoSgl(1:ns,1:4),2) > 0.0_pReal) nSources = (sum(rhoSgl(1:ns,1:2),2) * fEdgeMultiplication(instance) + sum(rhoSgl(1:ns,3:4),2)) & - / sum(rhoSgl(1:ns,1:4),2) * meshlength / lambda0(1:ns,instance) * sqrt(rhoForest(1:ns)) + / sum(rhoSgl(1:ns,1:4),2) * meshlength / lambda0(1:ns,instance)*sqrt(rhoForest(1:ns)) elsewhere nSources = meshlength / lambda0(1:ns,instance) * sqrt(rhoForest(1:ns)) endwhere @@ -4199,7 +2862,7 @@ else #ifndef _OPENMP if (iand(debug_level(debug_constitutive),debug_levelExtensive) /= 0_pInt & .and. ((debug_e == el .and. debug_i == ip .and. debug_g == ipc)& - .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then + .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt )) then write(6,'(a,/,4(12x,12(f12.5,1x),/))') '<< CONST >> sources', nSources write(6,*) endif @@ -4257,13 +2920,18 @@ if (.not. phase_localPlasticity(material_phase(ipc,ip,el))) then neighbor_el = mesh_ipNeighborhood(1,n,ip,el) neighbor_ip = mesh_ipNeighborhood(2,n,ip,el) neighbor_n = mesh_ipNeighborhood(3,n,ip,el) - + +#ifdef NEWSTATE + np = mappingConstitutive(2,ipc,neighbor_ip,neighbor_el) + no = mappingConstitutive(1,ipc,neighbor_ip,neighbor_el) +#endif + opposite_neighbor = n + mod(n,2_pInt) - mod(n+1_pInt,2_pInt) opposite_el = mesh_ipNeighborhood(1,opposite_neighbor,ip,el) opposite_ip = mesh_ipNeighborhood(2,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_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)) @@ -4274,9 +2942,13 @@ if (.not. phase_localPlasticity(material_phase(ipc,ip,el))) then !* FLUX FROM MY NEIGHBOR TO ME - !* This is only considered, if I have a neighbor of nonlocal plasticity (also nonlocal constitutive law with local properties) 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 + !* This is only considered, if I have a neighbor of nonlocal plasticity + !* (also nonlocal constitutive law with local properties) 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. neighbor_v = 0.0_pReal ! needed for check of sign change in flux density below @@ -4288,15 +2960,31 @@ if (.not. phase_localPlasticity(material_phase(ipc,ip,el))) then endif 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 +#ifdef NEWSTATE + forall (s = 1:ns, t = 1_pInt:4_pInt) + neighbor_v(s,t) = plasticState(np)%state0(iV(s,t,neighbor_instance),no) + neighbor_rhoSgl(s,t) = max(plasticState(np)%state0(iRhoU(s,t,neighbor_instance),no),0.0_pReal) + endforall +#else 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_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_instance)), 0.0_pReal) endforall +#endif else forall (s = 1:ns, t = 1_pInt:4_pInt) +#ifdef NEWSTATE + neighbor_v(s,t) = plasticState(np)%state(iV(s,t, neighbor_instance),no) + neighbor_rhoSgl(s,t) = max(plasticState(np)%state(iRhoU(s,t,neighbor_instance), no), & + 0.0_pReal) +#else neighbor_v(s,t) = state(ipc,neighbor_ip,neighbor_el)%p(iV(s,t,neighbor_instance)) - 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_instance)), & + 0.0_pReal) +#endif endforall endif where (neighbor_rhoSgl * mesh_ipVolume(neighbor_ip,neighbor_el) ** 0.667_pReal < significantN(instance) & @@ -4524,8 +3212,6 @@ endif end function constitutive_nonlocal_dotState -#endif - !********************************************************************* !* COMPATIBILITY UPDATE * !* Compatibility is defined as normalized product of signed cosine * @@ -4697,370 +3383,15 @@ enddo ! neighbor cycle compatibility(1:2,1:ns,1:ns,1:Nneighbors,i,e) = my_compatibility end subroutine constitutive_nonlocal_updateCompatibility + +!********************************************************************* +!* calculates quantities characterizing the microstructure * +!********************************************************************* #ifdef NEWSTATE -!********************************************************************* -!* calculates quantities characterizing the microstructure * -!********************************************************************* -function constitutive_nonlocal_dislocationstress(mapp, Fe, ipc, ip, el) - -use math, only: math_mul33x33, & - math_mul33x3, & - math_invert33, & - math_transpose33, & - pi -use mesh, only: mesh_NcpElems, & - mesh_maxNips, & - mesh_element, & - mesh_node0, & - mesh_cellCenterCoordinates, & - mesh_ipVolume, & - mesh_periodicSurface, & - FE_Nips, & - FE_geomtype -use material, only: homogenization_maxNgrains, & - material_phase, & - phase_localPlasticity, & - phase_plasticityInstance, & - plasticState -use lattice, only: lattice_mu, & - lattice_nu - -implicit none - -!*** input variables -integer(pInt), intent(in) :: ipc, & !< current grain ID - ip, & !< current integration point - el !< current element -real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & - Fe !< elastic deformation gradient -!type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & -! state !< microstructural state -integer(pInt),dimension(2,homogenization_maxngrains,mesh_maxNips,mesh_ncpelems), intent(in) :: & - mapp -!*** output variables -real(pReal), dimension(3,3) :: constitutive_nonlocal_dislocationstress - -!*** local variables -integer(pInt) neighbor_el, & !< element number of neighbor material point - neighbor_ip, & !< integration point of neighbor material point - instance, & !< my instance of this plasticity - neighbor_instance, & !< instance of this plasticity of neighbor material point - phase, & - neighbor_phase, & - ns, & !< total number of active slip systems at my material point - neighbor_ns, & !< total number of active slip systems at neighbor material point - c, & !< index of dilsocation character (edge, screw) - s, & !< slip system index - t, & !< index of dilsocation type (e+, e-, s+, s-, used e+, used e-, used s+, used s-) - dir, & - deltaX, deltaY, deltaZ, & - side, & - j -integer(pInt), dimension(2,3) :: periodicImages -real(pReal) x, y, z, & !< coordinates of connection vector in neighbor lattice frame - xsquare, ysquare, zsquare, & !< squares of respective coordinates - distance, & !< length of connection vector - segmentLength, & !< segment length of dislocations - lambda, & - R, Rsquare, Rcube, & - denominator, & - flipSign, & - neighbor_ipVolumeSideLength, & - detFe -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_neighborSlip, & !< connection vector between me and my neighbor in the slip system frame of my neighbor - maxCoord, minCoord, & - meshSize, & - coords, & !< x,y,z coordinates of cell center of ip volume - neighbor_coords !< x,y,z coordinates of cell center of neighbor ip volume -real(pReal), dimension(3,3) :: sigma, & !< dislocation stress for one slip system in neighbor material point's slip system frame - Tdislo_neighborLattice, & !< dislocation stress as 2nd Piola-Kirchhoff stress at neighbor material point - invFe, & !< inverse of my elastic deformation gradient - neighbor_invFe, & - neighborLattice2myLattice !< mapping from neighbor MPs lattice configuration to my lattice configuration -real(pReal), dimension(2,2,maxval(totalNslip)) :: & - neighbor_rhoExcess !< excess density at neighbor material point (edge/screw,mobile/dead,slipsystem) -real(pReal), dimension(2,maxval(totalNslip)) :: & - rhoExcessDead -real(pReal), dimension(totalNslip(phase_plasticityInstance(material_phase(ipc,ip,el))),8) :: & - rhoSgl ! single dislocation density (edge+, edge-, screw+, screw-, used edge+, used edge-, used screw+, used screw-) -logical inversionError - -phase = material_phase(ipc,ip,el) -instance = phase_plasticityInstance(phase) -ns = totalNslip(instance) - -!*** get basic states - -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+4_pInt) = state(ipc,ip,el)%p(iRhoB(s,t,instance)) - - rhoSgl(s,t) = max(plasticState(mapp(2,ipc,ip,el))%state(iRhoU(s,t,instance),mapp(1,ipc,ip,el)), & - 0.0_pReal) ! ensure positive single mobile densities - rhoSgl(s,t+4_pInt) = plasticState(mapp(2,ipc,ip,el))%state(iRhoB(s,t,instance),mapp(1,ipc,ip,el)) -endforall - - - -!*** calculate the dislocation stress of the neighboring excess dislocation densities -!*** zero for material points of local plasticity - -constitutive_nonlocal_dislocationstress = 0.0_pReal - -if (.not. phase_localPlasticity(phase)) then - call math_invert33(Fe(1:3,1:3,ipc,ip,el), invFe, detFe, inversionError) - - !* in case of periodic surfaces we have to find out how many periodic images in each direction we need - - do dir = 1_pInt,3_pInt - maxCoord(dir) = maxval(mesh_node0(dir,:)) - minCoord(dir) = minval(mesh_node0(dir,:)) - enddo - meshSize = maxCoord - minCoord - coords = mesh_cellCenterCoordinates(ip,el) - periodicImages = 0_pInt - do dir = 1_pInt,3_pInt - if (mesh_periodicSurface(dir)) then - periodicImages(1,dir) = floor((coords(dir) - cutoffRadius(instance) - minCoord(dir)) / meshSize(dir), pInt) - periodicImages(2,dir) = ceiling((coords(dir) + cutoffRadius(instance) - maxCoord(dir)) / meshSize(dir), pInt) - endif - enddo - - - !* loop through all material points (also through their periodic images if present), - !* but only consider nonlocal neighbors within a certain cutoff radius R - - do neighbor_el = 1_pInt,mesh_NcpElems -ipLoop: do neighbor_ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,neighbor_el))) - neighbor_phase = material_phase(ipc,neighbor_ip,neighbor_el) - if (phase_localPlasticity(neighbor_phase)) then - cycle - endif - neighbor_instance = phase_plasticityInstance(neighbor_phase) - neighbor_ns = totalNslip(neighbor_instance) - 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 - forall (s = 1_pInt:neighbor_ns, c = 1_pInt:2_pInt) - neighbor_rhoExcess(c,1,s) = plasticState(mapp(2,ipc,neighbor_ip,neighbor_el))% & - state(iRhoU(s,2*c-1,neighbor_instance),mapp(1,ipc,neighbor_ip,neighbor_el)) & ! positive mobiles - - plasticState(mapp(2,ipc,neighbor_ip,neighbor_el))% & - state(iRhoU(s,2*c,neighbor_instance),mapp(1,ipc,neighbor_ip,neighbor_el)) ! negative mobiles - neighbor_rhoExcess(c,2,s) = abs(plasticState(mapp(2,ipc,neighbor_ip,neighbor_el))% & - state(iRhoB(s,2*c-1,neighbor_instance),mapp(1,ipc,neighbor_ip,neighbor_el))) & ! positive deads - - abs(plasticState(mapp(2,ipc,neighbor_ip,neighbor_el))% & - state(iRhoB(s,2*c,neighbor_instance),mapp(1,ipc,neighbor_ip,neighbor_el))) ! negative deads - endforall - Tdislo_neighborLattice = 0.0_pReal - do deltaX = periodicImages(1,1),periodicImages(2,1) - do deltaY = periodicImages(1,2),periodicImages(2,2) - do deltaZ = periodicImages(1,3),periodicImages(2,3) - - - !* regular case - - if (neighbor_el /= el .or. neighbor_ip /= ip & - .or. deltaX /= 0_pInt .or. deltaY /= 0_pInt .or. deltaZ /= 0_pInt) then - - neighbor_coords = mesh_cellCenterCoordinates(neighbor_ip,neighbor_el) & - + [real(deltaX,pReal), real(deltaY,pReal), real(deltaZ,pReal)] * meshSize - connection = neighbor_coords - coords - distance = sqrt(sum(connection * connection)) - if (distance > cutoffRadius(instance)) then - cycle - endif - - - !* the segment length is the minimum of the third root of the control volume and the ip distance - !* this ensures, that the central MP never sits on a neighbor dislocation segment - - connection_neighborLattice = math_mul33x3(neighbor_invFe, connection) - segmentLength = min(neighbor_ipVolumeSideLength, distance) - - - !* loop through all slip systems of the neighbor material point - !* and add up the stress contributions from egde and screw excess on these slip systems (if significant) - - do s = 1_pInt,neighbor_ns - if (all(abs(neighbor_rhoExcess(:,:,s)) < significantRho(instance))) cycle ! not significant - - - !* 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_neighborLattice) - - - !* edge contribution to stress - sigma = 0.0_pReal - - x = connection_neighborSlip(1) - y = connection_neighborSlip(2) - z = connection_neighborSlip(3) - xsquare = x * x - ysquare = y * y - zsquare = z * z - - do j = 1_pInt,2_pInt - if (abs(neighbor_rhoExcess(1,j,s)) < significantRho(instance)) then - cycle - elseif (j > 1_pInt) then -! 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,2,neighbor_instance))) - x = connection_neighborSlip(1) & - + sign(0.5_pReal * segmentLength, & - plasticState(mapp(2,ipc,neighbor_ip,neighbor_el))% & - state(iRhoB(s,1,neighbor_instance),mapp(1,ipc,neighbor_ip,neighbor_el)) & - - plasticState(mapp(2,ipc,neighbor_ip,neighbor_el))% & - state(iRhoB(s,2,neighbor_instance),mapp(1,ipc,neighbor_ip,neighbor_el))) - xsquare = x * x - endif - - flipSign = sign(1.0_pReal, -y) - do side = 1_pInt,-1_pInt,-2_pInt - lambda = real(side,pReal) * 0.5_pReal * segmentLength - y - R = sqrt(xsquare + zsquare + lambda * lambda) - Rsquare = R * R - Rcube = Rsquare * R - denominator = R * (R + flipSign * lambda) - if (denominator == 0.0_pReal) exit ipLoop - - sigma(1,1) = sigma(1,1) - real(side,pReal) & - * flipSign * z / denominator & - * (1.0_pReal + xsquare / Rsquare + xsquare / denominator) & - * neighbor_rhoExcess(1,j,s) - sigma(2,2) = sigma(2,2) - real(side,pReal) & - * (flipSign * 2.0_pReal * lattice_nu(phase) * z / denominator + z * lambda / Rcube) & - * neighbor_rhoExcess(1,j,s) - sigma(3,3) = sigma(3,3) + real(side,pReal) & - * flipSign * z / denominator & - * (1.0_pReal - zsquare / Rsquare - zsquare / denominator) & - * neighbor_rhoExcess(1,j,s) - sigma(1,2) = sigma(1,2) + real(side,pReal) & - * x * z / Rcube * neighbor_rhoExcess(1,j,s) - sigma(1,3) = sigma(1,3) + real(side,pReal) & - * flipSign * x / denominator & - * (1.0_pReal - zsquare / Rsquare - zsquare / denominator) & - * neighbor_rhoExcess(1,j,s) - sigma(2,3) = sigma(2,3) - real(side,pReal) & - * (lattice_nu(phase) / R - zsquare / Rcube) * neighbor_rhoExcess(1,j,s) - enddo - enddo - - !* screw contribution to stress - - 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 - if (abs(neighbor_rhoExcess(2,j,s)) < significantRho(instance)) then - cycle - elseif (j > 1_pInt) then - y = connection_neighborSlip(2) & - + sign(0.5_pReal * segmentLength, & - plasticState(mapp(2,ipc,neighbor_ip,neighbor_el))% & - state(iRhoB(s,3,neighbor_instance),mapp(1,ipc,neighbor_ip,neighbor_el)) & - - plasticState(mapp(2,ipc,neighbor_ip,neighbor_el))% & - state(iRhoB(s,4,neighbor_instance),mapp(1,ipc,neighbor_ip,neighbor_el))) - ysquare = y * y - endif - - flipSign = sign(1.0_pReal, x) - do side = 1_pInt,-1_pInt,-2_pInt - lambda = x + real(side,pReal) * 0.5_pReal * segmentLength - R = sqrt(ysquare + zsquare + lambda * lambda) - Rsquare = R * R - Rcube = Rsquare * R - denominator = R * (R + flipSign * lambda) - if (denominator == 0.0_pReal) exit ipLoop - - sigma(1,2) = sigma(1,2) - real(side,pReal) * flipSign * z & - * (1.0_pReal - lattice_nu(phase)) / denominator & - * neighbor_rhoExcess(2,j,s) - sigma(1,3) = sigma(1,3) + real(side,pReal) * flipSign * y & - * (1.0_pReal - lattice_nu(phase)) / denominator & - * neighbor_rhoExcess(2,j,s) - enddo - enddo - - if (all(abs(sigma) < 1.0e-10_pReal)) cycle ! SIGMA IS NOT A REAL STRESS, THATS WHY WE NEED A REALLY SMALL VALUE HERE - - !* copy symmetric parts - - sigma(2,1) = sigma(1,2) - sigma(3,1) = sigma(1,3) - sigma(3,2) = sigma(2,3) - - - !* scale stresses and map them into the neighbor material point's lattice configuration - - sigma = sigma * lattice_mu(neighbor_phase) * burgers(s,neighbor_instance) & - / (4.0_pReal * pi * (1.0_pReal - lattice_nu(neighbor_phase))) & - * mesh_ipVolume(neighbor_ip,neighbor_el) / segmentLength ! reference volume is used here (according to the segment length calculation) - Tdislo_neighborLattice = Tdislo_neighborLattice & - + math_mul33x33(math_transpose33(lattice2slip(1:3,1:3,s,neighbor_instance)), & - math_mul33x33(sigma, lattice2slip(1:3,1:3,s,neighbor_instance))) - - enddo ! slip system loop - - - !* special case of central ip volume - !* only consider dead dislocations - !* we assume that they all sit at a distance equal to half the third root of V - !* in direction of the according slip direction - - else - - forall (s = 1_pInt:ns, c = 1_pInt:2_pInt) & - rhoExcessDead(c,s) = plasticState(mapp(2,ipc,ip,el))% & - state(iRhoB(s,2*c-1,instance),mapp(1,ipc,ip,el)) & ! 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) - + plasticState(mapp(2,ipc,neighbor_ip,neighbor_el))% & - state(iRhoB(s,2*c,instance),mapp(1,ipc,ip,el)) ! 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 - if (all(abs(rhoExcessDead(:,s)) < significantRho(instance))) cycle ! not significant - sigma = 0.0_pReal ! all components except for sigma13 are zero - sigma(1,3) = - (rhoExcessDead(1,s) + rhoExcessDead(2,s) * (1.0_pReal - lattice_nu(phase))) & - * neighbor_ipVolumeSideLength * lattice_mu(phase) * burgers(s,instance) & - / (sqrt(2.0_pReal) * pi * (1.0_pReal - lattice_nu(phase))) - sigma(3,1) = sigma(1,3) - - Tdislo_neighborLattice = Tdislo_neighborLattice & - + math_mul33x33(math_transpose33(lattice2slip(1:3,1:3,s,instance)), & - math_mul33x33(sigma, lattice2slip(1:3,1:3,s,instance))) - - enddo ! slip system loop - - endif - - enddo ! deltaZ loop - enddo ! deltaY loop - enddo ! deltaX loop - - - !* map the stress from the neighbor MP's lattice configuration into the deformed configuration - !* and back into my lattice configuration - - neighborLattice2myLattice = math_mul33x33(invFe, Fe(1:3,1:3,1,neighbor_ip,neighbor_el)) - constitutive_nonlocal_dislocationstress = constitutive_nonlocal_dislocationstress & - + math_mul33x33(neighborLattice2myLattice, & - math_mul33x33(Tdislo_neighborLattice, & - math_transpose33(neighborLattice2myLattice))) - - enddo ipLoop - enddo ! element loop - -endif - -end function constitutive_nonlocal_dislocationstress - +pure function constitutive_nonlocal_dislocationstress(Fe, ipc, ip, el) #else -!********************************************************************* -!* calculates quantities characterizing the microstructure * -!********************************************************************* pure function constitutive_nonlocal_dislocationstress(state, Fe, ipc, ip, el) - +#endif use math, only: math_mul33x33, & math_mul33x3, & math_invert33, & @@ -5077,6 +3408,10 @@ use mesh, only: mesh_NcpElems, & FE_geomtype use material, only: homogenization_maxNgrains, & material_phase, & +#ifdef NEWSTATE + plasticState, & + mappingConstitutive,& +#endif phase_localPlasticity, & phase_plasticityInstance use lattice, only: lattice_mu, & @@ -5090,8 +3425,10 @@ integer(pInt), intent(in) :: ipc, & el !< current element real(pReal), dimension(3,3,homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & Fe !< elastic deformation gradient +#ifndef NEWSTATE type(p_vec), dimension(homogenization_maxNgrains,mesh_maxNips,mesh_NcpElems), intent(in) :: & state !< microstructural state +#endif !*** output variables real(pReal), dimension(3,3) :: constitutive_nonlocal_dislocationstress @@ -5107,6 +3444,12 @@ integer(pInt) neighbor_el, & neighbor_ns, & !< total number of active slip systems at neighbor material point c, & !< index of dilsocation character (edge, screw) s, & !< slip system index +#ifdef NEWSTATE + o,& !< offset shortcut + no,& !< neighbour offset shortcut + p,& !< phase shortcut + np,& !< neighbour phase shortcut +#endif t, & !< index of dilsocation type (e+, e-, s+, s-, used e+, used e-, used s+, used s-) dir, & deltaX, deltaY, deltaZ, & @@ -5146,14 +3489,22 @@ logical inversionError phase = material_phase(ipc,ip,el) instance = phase_plasticityInstance(phase) ns = totalNslip(instance) - +#ifdef NEWSTATE + p = mappingConstitutive(2,1,ip,el) + o = mappingConstitutive(1,1,ip,el) +#endif !*** get basic states 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 +#ifdef NEWSTATE + rhoSgl(s,t) = max(plasticState(p)%state(iRhoU(s,t,instance),o), 0.0_pReal) ! ensure positive single mobile densities + rhoSgl(s,t+4_pInt) = plasticState(p)%state(iRhoB(s,t,instance),o) +#else + rhoSgl(s,t) = max(state(ipc,ip,el)%p(iRhoU(s,t,instance)), 0.0_pReal) ! ensure positive single mobile densities rhoSgl(s,t+4_pInt) = state(ipc,ip,el)%p(iRhoB(s,t,instance)) +#endif endforall @@ -5189,6 +3540,10 @@ if (.not. phase_localPlasticity(phase)) then do neighbor_el = 1_pInt,mesh_NcpElems ipLoop: do neighbor_ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,neighbor_el))) neighbor_phase = material_phase(ipc,neighbor_ip,neighbor_el) +#ifdef NEWSTATE + np = mappingConstitutive(2,1,neighbor_ip,neighbor_el) + no = mappingConstitutive(1,1,neighbor_ip,neighbor_el) +#endif if (phase_localPlasticity(neighbor_phase)) then cycle endif @@ -5196,11 +3551,19 @@ ipLoop: do neighbor_ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,neighbor_el)) neighbor_ns = totalNslip(neighbor_instance) 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 + 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 - - state(ipc,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c,neighbor_instance)) ! negative mobiles +#ifdef NEWSTATE + neighbor_rhoExcess(c,1,s) = plasticState(np)%state(iRhoU(s,2*c-1,neighbor_instance),no) & ! positive mobiles + - plasticState(np)%state(iRhoU(s,2*c,neighbor_instance),no) ! negative mobiles + neighbor_rhoExcess(c,2,s) = abs(plasticState(np)%state(iRhoB(s,2*c-1,neighbor_instance),no)) & ! positive deads + - abs(plasticState(np)%state(iRhoB(s,2*c,neighbor_instance),no)) ! negative deads +#else + neighbor_rhoExcess(c,1,s) = state(ipc,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c-1,neighbor_instance)) & ! positive mobiles + - state(ipc,neighbor_ip,neighbor_el)%p(iRhoU(s,2*c,neighbor_instance)) ! negative mobiles neighbor_rhoExcess(c,2,s) = abs(state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,2*c-1,neighbor_instance))) & ! positive deads - abs(state(ipc,neighbor_ip,neighbor_el)%p(iRhoB(s,2*c,neighbor_instance))) ! negative deads +#endif endforall Tdislo_neighborLattice = 0.0_pReal do deltaX = periodicImages(1,1),periodicImages(2,1) @@ -5256,10 +3619,17 @@ ipLoop: do neighbor_ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,neighbor_el)) if (abs(neighbor_rhoExcess(1,j,s)) < significantRho(instance)) then cycle elseif (j > 1_pInt) then +#ifdef NEWSTATE + x = connection_neighborSlip(1) & + + sign(0.5_pReal * segmentLength, & + plasticState(np)%state(iRhoB(s,1,neighbor_instance),no) & + - plasticState(np)%state(iRhoB(s,2,neighbor_instance),no)) +#else 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,2,neighbor_instance))) +#endif xsquare = x * x endif @@ -5301,10 +3671,18 @@ ipLoop: do neighbor_ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,neighbor_el)) if (abs(neighbor_rhoExcess(2,j,s)) < significantRho(instance)) then cycle elseif (j > 1_pInt) then +#ifdef NEWSTATE + y = connection_neighborSlip(2) & + + sign(0.5_pReal * segmentLength, & + plasticState(np)%state(iRhoB(s,3,neighbor_instance),no) & + - plasticState(np)%state(iRhoB(s,4,neighbor_instance),no)) +#else 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,4,neighbor_instance))) +#endif + ysquare = y * y endif @@ -5355,9 +3733,21 @@ ipLoop: do neighbor_ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,neighbor_el)) else 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) - + 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) - +#ifdef NEWSTATE + rhoExcessDead(c,s) = plasticState(p)%state(iRhoB(s,2*c-1,instance),o) & ! 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) + + plasticState(p)%state(iRhoB(s,2*c,instance),o) ! 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) +#else + 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) + + 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) +#endif do s = 1_pInt,ns if (all(abs(rhoExcessDead(:,s)) < significantRho(instance))) cycle ! not significant sigma = 0.0_pReal ! all components except for sigma13 are zero @@ -5394,7 +3784,10 @@ ipLoop: do neighbor_ip = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,neighbor_el)) endif end function constitutive_nonlocal_dislocationstress -#endif + +!!!!!!!!!! +!!!!!!!!!!! +!!!!!!!!!!!! #ifdef NEWSTATE !-------------------------------------------------------------------------------------------------- @@ -5893,7 +4286,7 @@ outputsLoop: do o = 1_pInt,phase_Noutput(material_phase(ipc,ip,el)) cs = cs + ns case(dislocationstress_ID) - sigma = constitutive_nonlocal_dislocationstress(mapping, Fe, ipc, ip, el) + sigma = constitutive_nonlocal_dislocationstress(Fe, ipc, ip, el) constitutive_nonlocal_postResults(cs+1_pInt) = sigma(1,1) constitutive_nonlocal_postResults(cs+2_pInt) = sigma(2,2) constitutive_nonlocal_postResults(cs+3_pInt) = sigma(3,3) @@ -5914,6 +4307,9 @@ enddo outputsLoop end function constitutive_nonlocal_postResults #else + +!!!!!!!!! +!!!!!!!!! !-------------------------------------------------------------------------------------------------- !> @brief return array of constitutive results !-------------------------------------------------------------------------------------------------- diff --git a/code/constitutive_phenopowerlaw.f90 b/code/constitutive_phenopowerlaw.f90 index a8637a3f0..6b58a6424 100644 --- a/code/constitutive_phenopowerlaw.f90 +++ b/code/constitutive_phenopowerlaw.f90 @@ -264,53 +264,53 @@ allocate(constitutive_phenopowerlaw_sizePostResults(maxNinstance), case ('(output)') select case(IO_lc(IO_stringValue(line,positions,2_pInt))) case ('resistance_slip') - constitutive_phenopowerlaw_outputID(constitutive_phenopowerlaw_Noutput(instance),instance) = resistance_slip_ID constitutive_phenopowerlaw_Noutput(instance) = constitutive_phenopowerlaw_Noutput(instance) + 1_pInt + constitutive_phenopowerlaw_outputID(constitutive_phenopowerlaw_Noutput(instance),instance) = resistance_slip_ID constitutive_phenopowerlaw_output(constitutive_phenopowerlaw_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('accumulatedshear_slip') - constitutive_phenopowerlaw_outputID(constitutive_phenopowerlaw_Noutput(instance),instance) = accumulatedshear_slip_ID constitutive_phenopowerlaw_Noutput(instance) = constitutive_phenopowerlaw_Noutput(instance) + 1_pInt + constitutive_phenopowerlaw_outputID(constitutive_phenopowerlaw_Noutput(instance),instance) = accumulatedshear_slip_ID constitutive_phenopowerlaw_output(constitutive_phenopowerlaw_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('shearrate_slip') - constitutive_phenopowerlaw_outputID(constitutive_phenopowerlaw_Noutput(instance),instance) = shearrate_slip_ID constitutive_phenopowerlaw_Noutput(instance) = constitutive_phenopowerlaw_Noutput(instance) + 1_pInt + constitutive_phenopowerlaw_outputID(constitutive_phenopowerlaw_Noutput(instance),instance) = shearrate_slip_ID constitutive_phenopowerlaw_output(constitutive_phenopowerlaw_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('resolvedstress_slip') - constitutive_phenopowerlaw_outputID(constitutive_phenopowerlaw_Noutput(instance),instance) = resolvedstress_slip_ID constitutive_phenopowerlaw_Noutput(instance) = constitutive_phenopowerlaw_Noutput(instance) + 1_pInt + constitutive_phenopowerlaw_outputID(constitutive_phenopowerlaw_Noutput(instance),instance) = resolvedstress_slip_ID constitutive_phenopowerlaw_output(constitutive_phenopowerlaw_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('totalshear') - constitutive_phenopowerlaw_outputID(constitutive_phenopowerlaw_Noutput(instance),instance) = totalshear_ID constitutive_phenopowerlaw_Noutput(instance) = constitutive_phenopowerlaw_Noutput(instance) + 1_pInt + constitutive_phenopowerlaw_outputID(constitutive_phenopowerlaw_Noutput(instance),instance) = totalshear_ID constitutive_phenopowerlaw_output(constitutive_phenopowerlaw_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('resistance_twin') - constitutive_phenopowerlaw_outputID(constitutive_phenopowerlaw_Noutput(instance),instance) = resistance_twin_ID constitutive_phenopowerlaw_Noutput(instance) = constitutive_phenopowerlaw_Noutput(instance) + 1_pInt + constitutive_phenopowerlaw_outputID(constitutive_phenopowerlaw_Noutput(instance),instance) = resistance_twin_ID constitutive_phenopowerlaw_output(constitutive_phenopowerlaw_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('accumulatedshear_twin') - constitutive_phenopowerlaw_outputID(constitutive_phenopowerlaw_Noutput(instance),instance) = accumulatedshear_twin_ID constitutive_phenopowerlaw_Noutput(instance) = constitutive_phenopowerlaw_Noutput(instance) + 1_pInt + constitutive_phenopowerlaw_outputID(constitutive_phenopowerlaw_Noutput(instance),instance) = accumulatedshear_twin_ID constitutive_phenopowerlaw_output(constitutive_phenopowerlaw_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('shearrate_twin') - constitutive_phenopowerlaw_outputID(constitutive_phenopowerlaw_Noutput(instance),instance) = shearrate_twin_ID constitutive_phenopowerlaw_Noutput(instance) = constitutive_phenopowerlaw_Noutput(instance) + 1_pInt + constitutive_phenopowerlaw_outputID(constitutive_phenopowerlaw_Noutput(instance),instance) = shearrate_twin_ID constitutive_phenopowerlaw_output(constitutive_phenopowerlaw_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('resolvedstress_twin') - constitutive_phenopowerlaw_outputID(constitutive_phenopowerlaw_Noutput(instance),instance) = resolvedstress_twin_ID constitutive_phenopowerlaw_Noutput(instance) = constitutive_phenopowerlaw_Noutput(instance) + 1_pInt + constitutive_phenopowerlaw_outputID(constitutive_phenopowerlaw_Noutput(instance),instance) = resolvedstress_twin_ID constitutive_phenopowerlaw_output(constitutive_phenopowerlaw_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('totalvolfrac') - constitutive_phenopowerlaw_outputID(constitutive_phenopowerlaw_Noutput(instance),instance) = totalvolfrac_ID constitutive_phenopowerlaw_Noutput(instance) = constitutive_phenopowerlaw_Noutput(instance) + 1_pInt + constitutive_phenopowerlaw_outputID(constitutive_phenopowerlaw_Noutput(instance),instance) = totalvolfrac_ID constitutive_phenopowerlaw_output(constitutive_phenopowerlaw_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case default diff --git a/code/constitutive_titanmod.f90 b/code/constitutive_titanmod.f90 index 3adc47820..637158d78 100644 --- a/code/constitutive_titanmod.f90 +++ b/code/constitutive_titanmod.f90 @@ -363,148 +363,148 @@ subroutine constitutive_titanmod_init(fileUnit) case ('(output)') select case(IO_lc(IO_stringValue(line,positions,2_pInt))) case ('rhoedge') - constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = rhoedge_ID constitutive_titanmod_Noutput(instance) = constitutive_titanmod_Noutput(instance) + 1_pInt + constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = rhoedge_ID constitutive_titanmod_output(constitutive_titanmod_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('rhoscrew') - constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = rhoscrew_ID constitutive_titanmod_Noutput(instance) = constitutive_titanmod_Noutput(instance) + 1_pInt + constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = rhoscrew_ID constitutive_titanmod_output(constitutive_titanmod_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('segment_edge') - constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = segment_edge_ID constitutive_titanmod_Noutput(instance) = constitutive_titanmod_Noutput(instance) + 1_pInt + constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = segment_edge_ID constitutive_titanmod_output(constitutive_titanmod_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('segment_screw') - constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = segment_screw_ID constitutive_titanmod_Noutput(instance) = constitutive_titanmod_Noutput(instance) + 1_pInt + constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = segment_screw_ID constitutive_titanmod_output(constitutive_titanmod_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('resistance_edge') - constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = resistance_edge_ID constitutive_titanmod_Noutput(instance) = constitutive_titanmod_Noutput(instance) + 1_pInt + constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = resistance_edge_ID constitutive_titanmod_output(constitutive_titanmod_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('resistance_screw') - constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = resistance_screw_ID constitutive_titanmod_Noutput(instance) = constitutive_titanmod_Noutput(instance) + 1_pInt + constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = resistance_screw_ID constitutive_titanmod_output(constitutive_titanmod_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('velocity_edge') - constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = velocity_edge_ID constitutive_titanmod_Noutput(instance) = constitutive_titanmod_Noutput(instance) + 1_pInt + constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = velocity_edge_ID constitutive_titanmod_output(constitutive_titanmod_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('velocity_screw') - constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = velocity_screw_ID constitutive_titanmod_Noutput(instance) = constitutive_titanmod_Noutput(instance) + 1_pInt + constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = velocity_screw_ID constitutive_titanmod_output(constitutive_titanmod_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('tau_slip') - constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = tau_slip_ID constitutive_titanmod_Noutput(instance) = constitutive_titanmod_Noutput(instance) + 1_pInt + constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = tau_slip_ID constitutive_titanmod_output(constitutive_titanmod_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('gdot_slip_edge') - constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = gdot_slip_edge_ID constitutive_titanmod_Noutput(instance) = constitutive_titanmod_Noutput(instance) + 1_pInt + constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = gdot_slip_edge_ID constitutive_titanmod_output(constitutive_titanmod_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('gdot_slip_screw') - constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = gdot_slip_screw_ID constitutive_titanmod_Noutput(instance) = constitutive_titanmod_Noutput(instance) + 1_pInt + constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = gdot_slip_screw_ID constitutive_titanmod_output(constitutive_titanmod_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('gdot_slip') - constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = gdot_slip_ID constitutive_titanmod_Noutput(instance) = constitutive_titanmod_Noutput(instance) + 1_pInt + constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = gdot_slip_ID constitutive_titanmod_output(constitutive_titanmod_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('stressratio_edge_p') - constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = stressratio_edge_p_ID constitutive_titanmod_Noutput(instance) = constitutive_titanmod_Noutput(instance) + 1_pInt + constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = stressratio_edge_p_ID constitutive_titanmod_output(constitutive_titanmod_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('stressratio_screw_p') - constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = stressratio_screw_p_ID constitutive_titanmod_Noutput(instance) = constitutive_titanmod_Noutput(instance) + 1_pInt + constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = stressratio_screw_p_ID constitutive_titanmod_output(constitutive_titanmod_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('shear_system') - constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = shear_system_ID constitutive_titanmod_Noutput(instance) = constitutive_titanmod_Noutput(instance) + 1_pInt + constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = shear_system_ID constitutive_titanmod_output(constitutive_titanmod_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('twin_fraction') - constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = twin_fraction_ID constitutive_titanmod_Noutput(instance) = constitutive_titanmod_Noutput(instance) + 1_pInt + constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = twin_fraction_ID constitutive_titanmod_output(constitutive_titanmod_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('shear_basal') - constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = shear_basal_ID constitutive_titanmod_Noutput(instance) = constitutive_titanmod_Noutput(instance) + 1_pInt + constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = shear_basal_ID constitutive_titanmod_output(constitutive_titanmod_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('shear_prism') - constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = shear_prism_ID constitutive_titanmod_Noutput(instance) = constitutive_titanmod_Noutput(instance) + 1_pInt + constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = shear_prism_ID constitutive_titanmod_output(constitutive_titanmod_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('shear_pyra') - constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = shear_pyra_ID constitutive_titanmod_Noutput(instance) = constitutive_titanmod_Noutput(instance) + 1_pInt + constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = shear_pyra_ID constitutive_titanmod_output(constitutive_titanmod_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('shear_pyrca') - constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = shear_pyrca_ID constitutive_titanmod_Noutput(instance) = constitutive_titanmod_Noutput(instance) + 1_pInt + constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = shear_pyrca_ID constitutive_titanmod_output(constitutive_titanmod_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('rhoedge_basal') - constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = rhoedge_basal_ID constitutive_titanmod_Noutput(instance) = constitutive_titanmod_Noutput(instance) + 1_pInt + constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = rhoedge_basal_ID constitutive_titanmod_output(constitutive_titanmod_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('rhoedge_prism') - constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = rhoedge_prism_ID constitutive_titanmod_Noutput(instance) = constitutive_titanmod_Noutput(instance) + 1_pInt + constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = rhoedge_prism_ID constitutive_titanmod_output(constitutive_titanmod_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('rhoedge_pyra') - constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = rhoedge_pyra_ID constitutive_titanmod_Noutput(instance) = constitutive_titanmod_Noutput(instance) + 1_pInt + constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = rhoedge_pyra_ID constitutive_titanmod_output(constitutive_titanmod_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('rhoedge_pyrca') - constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = rhoedge_pyrca_ID constitutive_titanmod_Noutput(instance) = constitutive_titanmod_Noutput(instance) + 1_pInt + constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = rhoedge_pyrca_ID constitutive_titanmod_output(constitutive_titanmod_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('rhoscrew_basal') - constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = rhoscrew_basal_ID constitutive_titanmod_Noutput(instance) = constitutive_titanmod_Noutput(instance) + 1_pInt + constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = rhoscrew_basal_ID constitutive_titanmod_output(constitutive_titanmod_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('rhoscrew_prism') - constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = rhoscrew_prism_ID constitutive_titanmod_Noutput(instance) = constitutive_titanmod_Noutput(instance) + 1_pInt + constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = rhoscrew_prism_ID constitutive_titanmod_output(constitutive_titanmod_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('rhoscrew_pyra') - constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = rhoscrew_pyra_ID constitutive_titanmod_Noutput(instance) = constitutive_titanmod_Noutput(instance) + 1_pInt + constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = rhoscrew_pyra_ID constitutive_titanmod_output(constitutive_titanmod_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('rhoscrew_pyrca') - constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = rhoscrew_pyrca_ID constitutive_titanmod_Noutput(instance) = constitutive_titanmod_Noutput(instance) + 1_pInt + constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = rhoscrew_pyrca_ID constitutive_titanmod_output(constitutive_titanmod_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) case ('shear_total') - constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = shear_total_ID constitutive_titanmod_Noutput(instance) = constitutive_titanmod_Noutput(instance) + 1_pInt + constitutive_titanmod_outputID(constitutive_titanmod_Noutput(instance),instance) = shear_total_ID constitutive_titanmod_output(constitutive_titanmod_Noutput(instance),instance) = & IO_lc(IO_stringValue(line,positions,2_pInt)) end select