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

Merge remote-tracking branch 'origin/some-polishing' into development

This commit is contained in:
Martin Diehl 2019-01-26 08:48:41 +01:00
parent bfe5438be7 e9087f83fe
commit ffd29bdcdc
4 changed files with 168 additions and 127 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

157
src/lattice.f90
View File

@ -869,8 +869,8 @@ subroutine lattice_initializeStructure(myPhase,CoverA)
math_mul33x3, & math_mul33x3, &
math_trace33, & math_trace33, &
math_symmetric33, & math_symmetric33, &
math_Mandel33to6, & math_sym33to6, &
math_Mandel3333to66, & math_sym3333to66, &
math_Voigt66to3333, & math_Voigt66to3333, &
math_axisAngleToR, & math_axisAngleToR, &
INRAD, & INRAD, &
@ -908,7 +908,7 @@ subroutine lattice_initializeStructure(myPhase,CoverA)
+ 6.0_pReal*lattice_C66(1,2,myPhase) & + 6.0_pReal*lattice_C66(1,2,myPhase) &
+ 2.0_pReal*lattice_C66(4,4,myPhase))! C12iso/(C11iso+C12iso) with C11iso=(3*C11+2*C12+4*C44)/5 and C12iso=(C11+4*C12-2*C44)/5 + 2.0_pReal*lattice_C66(4,4,myPhase))! C12iso/(C11iso+C12iso) with C11iso=(3*C11+2*C12+4*C44)/5 and C12iso=(C11+4*C12-2*C44)/5
lattice_C3333(1:3,1:3,1:3,1:3,myPhase) = math_Voigt66to3333(lattice_C66(1:6,1:6,myPhase)) ! Literature data is Voigt lattice_C3333(1:3,1:3,1:3,1:3,myPhase) = math_Voigt66to3333(lattice_C66(1:6,1:6,myPhase)) ! Literature data is Voigt
lattice_C66(1:6,1:6,myPhase) = math_Mandel3333to66(lattice_C3333(1:3,1:3,1:3,1:3,myPhase)) ! DAMASK uses Mandel lattice_C66(1:6,1:6,myPhase) = math_sym3333to66(lattice_C3333(1:3,1:3,1:3,1:3,myPhase)) ! DAMASK uses Mandel-weighting
do i = 1_pInt, 6_pInt do i = 1_pInt, 6_pInt
if (abs(lattice_C66(i,i,myPhase))<tol_math_check) & if (abs(lattice_C66(i,i,myPhase))<tol_math_check) &
call IO_error(135_pInt,el=i,ip=myPhase,ext_msg='matrix diagonal "el"ement of phase "ip"') call IO_error(135_pInt,el=i,ip=myPhase,ext_msg='matrix diagonal "el"ement of phase "ip"')
@ -1056,14 +1056,14 @@ subroutine lattice_initializeStructure(myPhase,CoverA)
enddo enddo
do j = 1_pInt,1_pInt+2_pInt*lattice_NnonSchmid(myPhase) do j = 1_pInt,1_pInt+2_pInt*lattice_NnonSchmid(myPhase)
lattice_Sslip_v(1:6,j,i,myPhase) = & lattice_Sslip_v(1:6,j,i,myPhase) = &
math_Mandel33to6(math_symmetric33(lattice_Sslip(1:3,1:3,j,i,myPhase))) math_sym33to6(math_symmetric33(lattice_Sslip(1:3,1:3,j,i,myPhase)))
enddo enddo
enddo enddo
do i = 1_pInt,myNcleavage ! store slip system vectors and Schmid matrix for my structure do i = 1_pInt,myNcleavage ! store slip system vectors and Schmid matrix for my structure
do j = 1_pInt,3_pInt do j = 1_pInt,3_pInt
lattice_Scleavage_v(1:6,j,i,myPhase) = & lattice_Scleavage_v(1:6,j,i,myPhase) = &
math_Mandel33to6(math_symmetric33(lattice_Scleavage(1:3,1:3,j,i,myPhase))) math_sym33to6(math_symmetric33(lattice_Scleavage(1:3,1:3,j,i,myPhase)))
enddo enddo
enddo enddo
@ -1366,11 +1366,14 @@ function lattice_characteristicShear_Twin(Ntwin,structure,CoverA) result(charact
4 & 4 &
],pInt),[LATTICE_HEX_NTWIN]) ! indicator to formulas below ],pInt),[LATTICE_HEX_NTWIN]) ! indicator to formulas below
if (len_trim(structure) /= 3_pInt) &
call IO_error(137_pInt,ext_msg='lattice_characteristicShear_Twin: '//trim(structure))
a = 0_pInt a = 0_pInt
myFamilies: do mf = 1_pInt,size(Ntwin,1) myFamilies: do mf = 1_pInt,size(Ntwin,1)
mySystems: do ms = 1_pInt,Ntwin(mf) mySystems: do ms = 1_pInt,Ntwin(mf)
a = a + 1_pInt a = a + 1_pInt
select case(trim(structure)) select case(structure(1:3))
case('fcc','bcc') case('fcc','bcc')
characteristicShear(a) = 0.5_pReal*sqrt(2.0_pReal) characteristicShear(a) = 0.5_pReal*sqrt(2.0_pReal)
case('hex') case('hex')
@ -1397,7 +1400,7 @@ end function lattice_characteristicShear_Twin
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Rotated elasticity matrices for twinning in Mandel notation !> @brief Rotated elasticity matrices for twinning in 66-vector notation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_C66_twin(Ntwin,C66,structure,CoverA) function lattice_C66_twin(Ntwin,C66,structure,CoverA)
use IO, only: & use IO, only: &
@ -1405,8 +1408,8 @@ function lattice_C66_twin(Ntwin,C66,structure,CoverA)
use math, only: & use math, only: &
INRAD, & INRAD, &
math_axisAngleToR, & math_axisAngleToR, &
math_Mandel3333to66, & math_sym3333to66, &
math_Mandel66to3333, & math_66toSym3333, &
math_rotate_forward3333 math_rotate_forward3333
implicit none implicit none
@ -1420,15 +1423,18 @@ function lattice_C66_twin(Ntwin,C66,structure,CoverA)
real(pReal), dimension(3,3) :: R real(pReal), dimension(3,3) :: R
integer(pInt) :: i integer(pInt) :: i
if (len_trim(structure) /= 3_pInt) &
call IO_error(137_pInt,ext_msg='lattice_C66_twin: '//trim(structure))
select case(trim(structure)) select case(structure(1:3))
case('fcc') case('fcc')
coordinateSystem = buildCoordinateSystem(Ntwin,LATTICE_FCC_NSLIPSYSTEM,LATTICE_FCC_SYSTEMTWIN,& coordinateSystem = buildCoordinateSystem(Ntwin,LATTICE_FCC_NSLIPSYSTEM,LATTICE_FCC_SYSTEMTWIN,&
trim(structure),0.0_pReal) trim(structure),0.0_pReal)
case('bcc') case('bcc')
coordinateSystem = buildCoordinateSystem(Ntwin,LATTICE_BCC_NSLIPSYSTEM,LATTICE_BCC_SYSTEMTWIN,& coordinateSystem = buildCoordinateSystem(Ntwin,LATTICE_BCC_NSLIPSYSTEM,LATTICE_BCC_SYSTEMTWIN,&
trim(structure),0.0_pReal) trim(structure),0.0_pReal)
case('hex','hexagonal') !ToDo: "No alias policy": long or short? case('hex')
coordinateSystem = buildCoordinateSystem(Ntwin,LATTICE_HEX_NSLIPSYSTEM,LATTICE_HEX_SYSTEMTWIN,& coordinateSystem = buildCoordinateSystem(Ntwin,LATTICE_HEX_NSLIPSYSTEM,LATTICE_HEX_SYSTEMTWIN,&
'hex',cOverA) 'hex',cOverA)
case default case default
@ -1437,18 +1443,17 @@ function lattice_C66_twin(Ntwin,C66,structure,CoverA)
do i = 1, sum(Ntwin) do i = 1, sum(Ntwin)
R = math_axisAngleToR(coordinateSystem(1:3,2,i), 180.0_pReal * INRAD) ! ToDo: Why always 180 deg? R = math_axisAngleToR(coordinateSystem(1:3,2,i), 180.0_pReal * INRAD) ! ToDo: Why always 180 deg?
lattice_C66_twin(1:6,1:6,i) = math_Mandel3333to66(math_rotate_forward3333(math_Mandel66to3333(C66),R)) lattice_C66_twin(1:6,1:6,i) = math_sym3333to66(math_rotate_forward3333(math_66toSym3333(C66),R))
enddo enddo
end function lattice_C66_twin end function lattice_C66_twin
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Rotated elasticity matrices for transformation in Mandel notation !> @brief Rotated elasticity matrices for transformation in 66-vector notation
!> ToDo: Completely untested and incomplete and undocumented !> ToDo: Completely untested and incomplete and undocumented
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_C66_trans(Ntrans,C_parent66, & function lattice_C66_trans(Ntrans,C_parent66,structure_target, &
structure_target, & CoverA_trans,a_bcc,a_fcc)
CoverA_trans,a_bcc,a_fcc)
use prec, only: & use prec, only: &
tol_math_check tol_math_check
use IO, only: & use IO, only: &
@ -1465,21 +1470,25 @@ function lattice_C66_trans(Ntrans,C_parent66, &
math_crossproduct math_crossproduct
implicit none implicit none
integer(pInt), dimension(:), intent(in) :: Ntrans !< number of active twin systems per family integer(pInt), dimension(:), intent(in) :: Ntrans !< number of active twin systems per family
character(len=*), intent(in) :: & character(len=*), intent(in) :: &
structure_target !< lattice structure structure_target !< lattice structure
real(pReal), dimension(6,6), intent(in) :: C_parent66 real(pReal), dimension(6,6), intent(in) :: C_parent66
real(pReal), dimension(6,6) :: C_bar66, C_target_unrotated66 real(pReal), dimension(6,6) :: C_bar66, C_target_unrotated66
real(pReal), dimension(3,3,3,3) :: C_target_unrotated real(pReal), dimension(3,3,3,3) :: C_target_unrotated
real(pReal), dimension(6,6,sum(Ntrans)) :: lattice_C66_trans real(pReal), dimension(6,6,sum(Ntrans)) :: lattice_C66_trans
real(pReal), dimension(3,3,sum(Ntrans)) :: Q,S real(pReal), dimension(3,3,sum(Ntrans)) :: Q,S
real(pReal) :: a_bcc, a_fcc, CoverA_trans real(pReal) :: a_bcc, a_fcc, CoverA_trans
integer(pInt) :: i integer(pInt) :: i
if (len_trim(structure_target) /= 3_pInt) &
call IO_error(137_pInt,ext_msg='lattice_C66_trans (target): '//trim(structure_target))
!ToDo: add checks for CoverA_trans,a_fcc,a_bcc
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! elasticity matrix of the target phase in cube orientation ! elasticity matrix of the target phase in cube orientation
if (trim(structure_target) == 'hex') then if (structure_target(1:3) == 'hex') then
C_bar66(1,1) = (C_parent66(1,1) + C_parent66(1,2) + 2.0_pReal*C_parent66(4,4))/2.0_pReal C_bar66(1,1) = (C_parent66(1,1) + C_parent66(1,2) + 2.0_pReal*C_parent66(4,4))/2.0_pReal
C_bar66(1,2) = (C_parent66(1,1) + 5.0_pReal*C_parent66(1,2) - 2.0_pReal*C_parent66(4,4))/6.0_pReal C_bar66(1,2) = (C_parent66(1,1) + 5.0_pReal*C_parent66(1,2) - 2.0_pReal*C_parent66(4,4))/6.0_pReal
C_bar66(3,3) = (C_parent66(1,1) + 2.0_pReal*C_parent66(1,2) + 4.0_pReal*C_parent66(4,4))/3.0_pReal C_bar66(3,3) = (C_parent66(1,1) + 2.0_pReal*C_parent66(1,2) + 4.0_pReal*C_parent66(4,4))/3.0_pReal
@ -1494,10 +1503,10 @@ function lattice_C66_trans(Ntrans,C_parent66, &
C_target_unrotated66(3,3) = C_bar66(3,3) C_target_unrotated66(3,3) = C_bar66(3,3)
C_target_unrotated66(4,4) = C_bar66(4,4) - C_bar66(1,4)**2.0_pReal/(0.5_pReal*(C_bar66(1,1) - C_bar66(1,2))) C_target_unrotated66(4,4) = C_bar66(4,4) - C_bar66(1,4)**2.0_pReal/(0.5_pReal*(C_bar66(1,1) - C_bar66(1,2)))
C_target_unrotated66 = lattice_symmetrizeC66(LATTICE_HEX_ID,C_target_unrotated66) C_target_unrotated66 = lattice_symmetrizeC66(LATTICE_HEX_ID,C_target_unrotated66)
elseif (trim(structure_target) == 'bcc') then elseif (structure_target(1:3) == 'bcc') then
C_target_unrotated66 = C_parent66 C_target_unrotated66 = C_parent66
else else
write(6,*) "Mist" call IO_error(137_pInt,ext_msg='lattice_C66_trans (target): '//trim(structure_target))
endif endif
@ -1511,7 +1520,7 @@ function lattice_C66_trans(Ntrans,C_parent66, &
do i = 1, sum(Ntrans) do i = 1, sum(Ntrans)
lattice_C66_trans(1:6,1:6,i) = math_Mandel3333to66(math_rotate_forward3333(C_target_unrotated,Q(1:3,1:3,i))) lattice_C66_trans(1:6,1:6,i) = math_Mandel3333to66(math_rotate_forward3333(C_target_unrotated,Q(1:3,1:3,i)))
enddo enddo
end function end function lattice_C66_trans
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -1584,14 +1593,17 @@ function lattice_interaction_SlipSlip(Nslip,interactionValues,structure) result(
integer(pInt), dimension(:), allocatable :: NslipMax integer(pInt), dimension(:), allocatable :: NslipMax
integer(pInt), dimension(:,:), allocatable :: interactionTypes integer(pInt), dimension(:,:), allocatable :: interactionTypes
select case(structure) if (len_trim(structure) /= 3_pInt) &
call IO_error(137_pInt,ext_msg='lattice_interaction_SlipSlip: '//trim(structure))
select case(structure(1:3))
case('fcc') case('fcc')
interactionTypes = LATTICE_FCC_INTERACTIONSLIPSLIP interactionTypes = LATTICE_FCC_INTERACTIONSLIPSLIP
NslipMax = LATTICE_FCC_NSLIPSYSTEM NslipMax = LATTICE_FCC_NSLIPSYSTEM
case('bcc') case('bcc')
interactionTypes = LATTICE_BCC_INTERACTIONSLIPSLIP interactionTypes = LATTICE_BCC_INTERACTIONSLIPSLIP
NslipMax = LATTICE_BCC_NSLIPSYSTEM NslipMax = LATTICE_BCC_NSLIPSYSTEM
case('hex','hexagonal') ! ToDo: "No alias policy": long or short? case('hex')
interactionTypes = LATTICE_HEX_INTERACTIONSLIPSLIP interactionTypes = LATTICE_HEX_INTERACTIONSLIPSLIP
NslipMax = LATTICE_HEX_NSLIPSYSTEM NslipMax = LATTICE_HEX_NSLIPSYSTEM
case('bct') case('bct')
@ -1688,14 +1700,17 @@ function lattice_interaction_TwinTwin(Ntwin,interactionValues,structure) result(
20,20,20,20,20,20, 19,19,19,19,19,19, 18,18,18,18,18,18, 17,17,17,17,17,16 & 20,20,20,20,20,20, 19,19,19,19,19,19, 18,18,18,18,18,18, 17,17,17,17,17,16 &
],pInt),shape(HEX_INTERACTIONTWINTWIN),order=[2,1]) !< Twin-twin interaction types for hex ],pInt),shape(HEX_INTERACTIONTWINTWIN),order=[2,1]) !< Twin-twin interaction types for hex
select case(structure) if (len_trim(structure) /= 3_pInt) &
call IO_error(137_pInt,ext_msg='lattice_interaction_TwinTwin: '//trim(structure))
select case(structure(1:3))
case('fcc') case('fcc')
interactionTypes = FCC_INTERACTIONTWINTWIN interactionTypes = FCC_INTERACTIONTWINTWIN
NtwinMax = LATTICE_FCC_NTWINSYSTEM NtwinMax = LATTICE_FCC_NTWINSYSTEM
case('bcc') case('bcc')
interactionTypes = BCC_INTERACTIONTWINTWIN interactionTypes = BCC_INTERACTIONTWINTWIN
NtwinMax = LATTICE_BCC_NTWINSYSTEM NtwinMax = LATTICE_BCC_NTWINSYSTEM
case('hex','hexagonal') ! ToDo: "No alias policy": long or short? case('hex')
interactionTypes = HEX_INTERACTIONTWINTWIN interactionTypes = HEX_INTERACTIONTWINTWIN
NtwinMax = LATTICE_HEX_NTWINSYSTEM NtwinMax = LATTICE_HEX_NTWINSYSTEM
case default case default
@ -1740,7 +1755,10 @@ function lattice_interaction_TransTrans(Ntrans,interactionValues,structure) resu
2,2,2,2,2,2,2,2,2,1,1,1 & 2,2,2,2,2,2,2,2,2,1,1,1 &
],pInt),shape(FCC_INTERACTIONTRANSTRANS),order=[2,1]) !< Trans-trans interaction types for fcc ],pInt),shape(FCC_INTERACTIONTRANSTRANS),order=[2,1]) !< Trans-trans interaction types for fcc
if (trim(structure) == 'fcc') then if (len_trim(structure) /= 3_pInt) &
call IO_error(137_pInt,ext_msg='lattice_interaction_TransTrans: '//trim(structure))
if(structure(1:3) == 'fcc') then
interactionTypes = FCC_INTERACTIONTRANSTRANS interactionTypes = FCC_INTERACTIONTRANSTRANS
NtransMax = LATTICE_FCC_NTRANSSYSTEM NtransMax = LATTICE_FCC_NTRANSSYSTEM
else else
@ -1870,8 +1888,10 @@ function lattice_interaction_SlipTwin(Nslip,Ntwin,interactionValues,structure) r
! !
],pInt),shape(HEX_INTERACTIONSLIPTWIN),order=[2,1]) !< Slip-twin interaction types for hex ],pInt),shape(HEX_INTERACTIONSLIPTWIN),order=[2,1]) !< Slip-twin interaction types for hex
if (len_trim(structure) /= 3_pInt) &
select case(structure) call IO_error(137_pInt,ext_msg='lattice_interaction_SlipTwin: '//trim(structure))
select case(structure(1:3))
case('fcc') case('fcc')
interactionTypes = FCC_INTERACTIONSLIPTWIN interactionTypes = FCC_INTERACTIONSLIPTWIN
NslipMax = LATTICE_FCC_NSLIPSYSTEM NslipMax = LATTICE_FCC_NSLIPSYSTEM
@ -1880,7 +1900,7 @@ function lattice_interaction_SlipTwin(Nslip,Ntwin,interactionValues,structure) r
interactionTypes = BCC_INTERACTIONSLIPTWIN interactionTypes = BCC_INTERACTIONSLIPTWIN
NslipMax = LATTICE_BCC_NSLIPSYSTEM NslipMax = LATTICE_BCC_NSLIPSYSTEM
NtwinMax = LATTICE_BCC_NTWINSYSTEM NtwinMax = LATTICE_BCC_NTWINSYSTEM
case('hex','hexagonal') ! ToDo: "No alias policy": long or short? case('hex')
interactionTypes = HEX_INTERACTIONSLIPTWIN interactionTypes = HEX_INTERACTIONSLIPTWIN
NslipMax = LATTICE_HEX_NSLIPSYSTEM NslipMax = LATTICE_HEX_NSLIPSYSTEM
NtwinMax = LATTICE_HEX_NTWINSYSTEM NtwinMax = LATTICE_HEX_NTWINSYSTEM
@ -1936,7 +1956,10 @@ function lattice_interaction_SlipTrans(Nslip,Ntrans,interactionValues,structure)
4,4,4,4,4,4,4,4,4,4,4,4 & 4,4,4,4,4,4,4,4,4,4,4,4 &
],pInt),shape(FCC_INTERACTIONSLIPTRANS),order=[2,1]) !< Slip-trans interaction types for fcc ],pInt),shape(FCC_INTERACTIONSLIPTRANS),order=[2,1]) !< Slip-trans interaction types for fcc
select case(structure) if (len_trim(structure) /= 3_pInt) &
call IO_error(137_pInt,ext_msg='lattice_interaction_SlipTrans: '//trim(structure))
select case(structure(1:3))
case('fcc') case('fcc')
interactionTypes = FCC_INTERACTIONSLIPTRANS interactionTypes = FCC_INTERACTIONSLIPTRANS
NslipMax = LATTICE_FCC_NSLIPSYSTEM NslipMax = LATTICE_FCC_NSLIPSYSTEM
@ -2005,8 +2028,11 @@ function lattice_interaction_TwinSlip(Ntwin,Nslip,interactionValues,structure) r
4, 4, 4, 8, 8, 8, 12,12,12, 16,16,16,16,16,16, 20,20,20,20,20,20,20,20,20,20,20,20, 24,24,24,24,24,24, & 4, 4, 4, 8, 8, 8, 12,12,12, 16,16,16,16,16,16, 20,20,20,20,20,20,20,20,20,20,20,20, 24,24,24,24,24,24, &
4, 4, 4, 8, 8, 8, 12,12,12, 16,16,16,16,16,16, 20,20,20,20,20,20,20,20,20,20,20,20, 24,24,24,24,24,24 & 4, 4, 4, 8, 8, 8, 12,12,12, 16,16,16,16,16,16, 20,20,20,20,20,20,20,20,20,20,20,20, 24,24,24,24,24,24 &
],pInt),shape(HEX_INTERACTIONTWINSLIP),order=[2,1]) !< Twin-twin interaction types for hex ],pInt),shape(HEX_INTERACTIONTWINSLIP),order=[2,1]) !< Twin-twin interaction types for hex
if (len_trim(structure) /= 3_pInt) &
call IO_error(137_pInt,ext_msg='lattice_interaction_TwinSlip: '//trim(structure))
select case(structure) select case(structure(1:3))
case('fcc') case('fcc')
interactionTypes = FCC_INTERACTIONTWINSLIP interactionTypes = FCC_INTERACTIONTWINSLIP
NtwinMax = LATTICE_FCC_NTWINSYSTEM NtwinMax = LATTICE_FCC_NTWINSYSTEM
@ -2015,7 +2041,7 @@ function lattice_interaction_TwinSlip(Ntwin,Nslip,interactionValues,structure) r
interactionTypes = BCC_INTERACTIONTWINSLIP interactionTypes = BCC_INTERACTIONTWINSLIP
NtwinMax = LATTICE_BCC_NTWINSYSTEM NtwinMax = LATTICE_BCC_NTWINSYSTEM
NslipMax = LATTICE_BCC_NSLIPSYSTEM NslipMax = LATTICE_BCC_NSLIPSYSTEM
case('hex','hexagonal') ! ToDo: "No alias policy": long or short? case('hex')
interactionTypes = HEX_INTERACTIONTWINSLIP interactionTypes = HEX_INTERACTIONTWINSLIP
NtwinMax = LATTICE_HEX_NTWINSYSTEM NtwinMax = LATTICE_HEX_NTWINSYSTEM
NslipMax = LATTICE_HEX_NSLIPSYSTEM NslipMax = LATTICE_HEX_NSLIPSYSTEM
@ -2051,15 +2077,18 @@ function lattice_SchmidMatrix_slip(Nslip,structure,cOverA) result(SchmidMatrix)
real(pReal), dimension(:,:), allocatable :: slipSystems real(pReal), dimension(:,:), allocatable :: slipSystems
integer(pInt), dimension(:), allocatable :: NslipMax integer(pInt), dimension(:), allocatable :: NslipMax
integer(pInt) :: i integer(pInt) :: i
if (len_trim(structure) /= 3_pInt) &
call IO_error(137_pInt,ext_msg='lattice_SchmidMatrix_slip: '//trim(structure))
select case(structure) select case(structure(1:3))
case('fcc') case('fcc')
NslipMax = LATTICE_FCC_NSLIPSYSTEM NslipMax = LATTICE_FCC_NSLIPSYSTEM
slipSystems = LATTICE_FCC_SYSTEMSLIP slipSystems = LATTICE_FCC_SYSTEMSLIP
case('bcc') case('bcc')
NslipMax = LATTICE_BCC_NSLIPSYSTEM NslipMax = LATTICE_BCC_NSLIPSYSTEM
slipSystems = LATTICE_BCC_SYSTEMSLIP slipSystems = LATTICE_BCC_SYSTEMSLIP
case('hex','hexagonal') ! ToDo: "No alias policy": long or short? case('hex')
NslipMax = LATTICE_HEX_NSLIPSYSTEM NslipMax = LATTICE_HEX_NSLIPSYSTEM
slipSystems = LATTICE_HEX_SYSTEMSLIP slipSystems = LATTICE_HEX_SYSTEMSLIP
case('bct') case('bct')
@ -2109,14 +2138,17 @@ function lattice_SchmidMatrix_twin(Ntwin,structure,cOverA) result(SchmidMatrix)
integer(pInt), dimension(:), allocatable :: NtwinMax integer(pInt), dimension(:), allocatable :: NtwinMax
integer(pInt) :: i integer(pInt) :: i
select case(structure) if (len_trim(structure) /= 3_pInt) &
call IO_error(137_pInt,ext_msg='lattice_SchmidMatrix_twin: '//trim(structure))
select case(structure(1:3))
case('fcc') case('fcc')
NtwinMax = LATTICE_FCC_NTWINSYSTEM NtwinMax = LATTICE_FCC_NTWINSYSTEM
twinSystems = LATTICE_FCC_SYSTEMTWIN twinSystems = LATTICE_FCC_SYSTEMTWIN
case('bcc') case('bcc')
NtwinMax = LATTICE_BCC_NTWINSYSTEM NtwinMax = LATTICE_BCC_NTWINSYSTEM
twinSystems = LATTICE_BCC_SYSTEMTWIN twinSystems = LATTICE_BCC_SYSTEMTWIN
case('hex','hexagonal') !ToDo: "No alias policy": long or short? case('hex')
NtwinMax = LATTICE_HEX_NTWINSYSTEM NtwinMax = LATTICE_HEX_NTWINSYSTEM
twinSystems = LATTICE_HEX_SYSTEMTWIN twinSystems = LATTICE_HEX_SYSTEMTWIN
case default case default
@ -2162,11 +2194,17 @@ function lattice_SchmidMatrix_trans(Ntrans,structure_target,cOverA,a_bcc,a_fcc)
real(pReal), dimension(3,3,sum(Ntrans)) :: devNull real(pReal), dimension(3,3,sum(Ntrans)) :: devNull
real(pReal) :: a_bcc, a_fcc real(pReal) :: a_bcc, a_fcc
! ToDo: Error checking!!!!!!!!!!!!!!!!!!!
if (len_trim(structure_target) /= 3_pInt) &
call IO_error(137_pInt,ext_msg='lattice_SchmidMatrix_trans (target): '//trim(structure_target))
if (structure_target(1:3) /= 'bcc' .and. structure_target(1:3) /= 'hex') &
call IO_error(137_pInt,ext_msg='lattice_SchmidMatrix_trans (target): '//trim(structure_target))
!ToDo: add checks for CoverA_trans,a_fcc,a_bcc
call buildTransformationSystem(devNull,SchmidMatrix,Ntrans,cOverA,a_fcc,a_bcc) call buildTransformationSystem(devNull,SchmidMatrix,Ntrans,cOverA,a_fcc,a_bcc)
end function lattice_SchmidMatrix_trans
end function lattice_SchmidMatrix_trans
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -2189,8 +2227,11 @@ function lattice_SchmidMatrix_cleavage(Ncleavage,structure,cOverA) result(Schmid
real(pReal), dimension(:,:), allocatable :: cleavageSystems real(pReal), dimension(:,:), allocatable :: cleavageSystems
integer(pInt), dimension(:), allocatable :: NcleavageMax integer(pInt), dimension(:), allocatable :: NcleavageMax
integer(pInt) :: i integer(pInt) :: i
if (len_trim(structure) /= 3_pInt) &
call IO_error(137_pInt,ext_msg='lattice_SchmidMatrix_cleavage: '//trim(structure))
select case(structure) select case(structure(1:3))
case('iso') case('iso')
NcleavageMax = LATTICE_ISO_NCLEAVAGESYSTEM NcleavageMax = LATTICE_ISO_NCLEAVAGESYSTEM
cleavageSystems = LATTICE_ISO_SYSTEMCLEAVAGE cleavageSystems = LATTICE_ISO_SYSTEMCLEAVAGE
@ -2203,7 +2244,7 @@ function lattice_SchmidMatrix_cleavage(Ncleavage,structure,cOverA) result(Schmid
case('bcc') case('bcc')
NcleavageMax = LATTICE_BCC_NCLEAVAGESYSTEM NcleavageMax = LATTICE_BCC_NCLEAVAGESYSTEM
cleavageSystems = LATTICE_BCC_SYSTEMCLEAVAGE cleavageSystems = LATTICE_BCC_SYSTEMCLEAVAGE
case('hex','hexagonal') !ToDo: "No alias policy": long or short? case('hex')
NcleavageMax = LATTICE_HEX_NCLEAVAGESYSTEM NcleavageMax = LATTICE_HEX_NCLEAVAGESYSTEM
cleavageSystems = LATTICE_HEX_SYSTEMCLEAVAGE cleavageSystems = LATTICE_HEX_SYSTEMCLEAVAGE
case default case default
@ -2246,14 +2287,17 @@ function lattice_forestProjection(Nslip,structure,cOverA) result(projection)
integer(pInt), dimension(:), allocatable :: NslipMax integer(pInt), dimension(:), allocatable :: NslipMax
integer(pInt) :: i, j integer(pInt) :: i, j
select case(structure) if (len_trim(structure) /= 3_pInt) &
call IO_error(137_pInt,ext_msg='lattice_forestProjection: '//trim(structure))
select case(structure(1:3))
case('fcc') case('fcc')
NslipMax = LATTICE_FCC_NSLIPSYSTEM NslipMax = LATTICE_FCC_NSLIPSYSTEM
slipSystems = LATTICE_FCC_SYSTEMSLIP slipSystems = LATTICE_FCC_SYSTEMSLIP
case('bcc') case('bcc')
NslipMax = LATTICE_BCC_NSLIPSYSTEM NslipMax = LATTICE_BCC_NSLIPSYSTEM
slipSystems = LATTICE_BCC_SYSTEMSLIP slipSystems = LATTICE_BCC_SYSTEMSLIP
case('hex','hexagonal') ! ToDo: "No alias policy": long or short? case('hex')
NslipMax = LATTICE_HEX_NSLIPSYSTEM NslipMax = LATTICE_HEX_NSLIPSYSTEM
slipSystems = LATTICE_HEX_SYSTEMSLIP slipSystems = LATTICE_HEX_SYSTEMSLIP
case('bct') case('bct')
@ -2346,9 +2390,11 @@ function buildCoordinateSystem(active,complete,system,structure,cOverA)
f, & !< index of my family f, & !< index of my family
s !< index of my system in current family s !< index of my system in current family
if (trim(structure) == 'bct' .and. cOverA > 2.0_pReal) & if (len_trim(structure) /= 3_pInt) &
call IO_error(137_pInt,ext_msg='buildCoordinateSystem: '//trim(structure))
if (trim(structure(1:3)) == 'bct' .and. cOverA > 2.0_pReal) &
call IO_error(131_pInt,ext_msg='buildCoordinateSystem:'//trim(structure)) call IO_error(131_pInt,ext_msg='buildCoordinateSystem:'//trim(structure))
if (trim(structure) == 'hex' .and. (cOverA < 1.0_pReal .or. cOverA > 2.0_pReal)) & if (trim(structure(1:3)) == 'hex' .and. (cOverA < 1.0_pReal .or. cOverA > 2.0_pReal)) &
call IO_error(131_pInt,ext_msg='buildCoordinateSystem:'//trim(structure)) call IO_error(131_pInt,ext_msg='buildCoordinateSystem:'//trim(structure))
a = 0_pInt a = 0_pInt
@ -2357,7 +2403,7 @@ function buildCoordinateSystem(active,complete,system,structure,cOverA)
a = a + 1_pInt a = a + 1_pInt
c = sum(complete(1:f-1))+s c = sum(complete(1:f-1))+s
select case(trim(structure)) select case(trim(structure(1:3)))
case ('fcc','bcc','iso','ort','bct') case ('fcc','bcc','iso','ort','bct')
direction = system(1:3,c) direction = system(1:3,c)
@ -2391,7 +2437,7 @@ end function buildCoordinateSystem
!> @brief Helper function to define transformation systems !> @brief Helper function to define transformation systems
! Needed to calculate Schmid matrix and rotated stiffness matrices. ! Needed to calculate Schmid matrix and rotated stiffness matrices.
! @details: set c/a = 0.0 for fcc -> bcc transformation ! @details: set c/a = 0.0 for fcc -> bcc transformation
! set a_bcc = 0.0 for fcc -> bcc transformation ! set a_bcc = 0.0 for fcc -> hex transformation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc) subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc)
use prec, only: & use prec, only: &
@ -2493,7 +2539,6 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc)
if (size(Ntrans) < 1_pInt .or. size(Ntrans) > 1_pInt) print*, 'mist' ! ToDo if (size(Ntrans) < 1_pInt .or. size(Ntrans) > 1_pInt) print*, 'mist' ! ToDo
if (a_bcc > 0.0_pReal .and. dEq0(cOverA)) then ! fcc -> bcc transformation if (a_bcc > 0.0_pReal .and. dEq0(cOverA)) then ! fcc -> bcc transformation
if (a_bcc <= 0.0_pReal) print*, 'mist' ! ToDo
do i = 1_pInt,sum(Ntrans) do i = 1_pInt,sum(Ntrans)
R = math_axisAngleToR(lattice_fccTobcc_systemTrans(1:3,i), & R = math_axisAngleToR(lattice_fccTobcc_systemTrans(1:3,i), &
lattice_fccTobcc_systemTrans(4,i)*INRAD) lattice_fccTobcc_systemTrans(4,i)*INRAD)
@ -2525,6 +2570,8 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc)
Q(1:3,3,i) = z Q(1:3,3,i) = z
S(1:3,1:3,i) = math_mul33x33(Q(1:3,1:3,i), math_mul33x33(math_mul33x33(sd,ss), transpose(Q(1:3,1:3,i)))) - MATH_I3 ! ToDo: This is of interest for the Schmid matrix only S(1:3,1:3,i) = math_mul33x33(Q(1:3,1:3,i), math_mul33x33(math_mul33x33(sd,ss), transpose(Q(1:3,1:3,i)))) - MATH_I3 ! ToDo: This is of interest for the Schmid matrix only
enddo enddo
else
call IO_error(0_pInt) !ToDo: define error
endif endif
end subroutine buildTransformationSystem end subroutine buildTransformationSystem

52
src/plastic_disloUCLA.f90
View File

@ -93,7 +93,7 @@ module plastic_disloUCLA
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! containers for parameters and state ! containers for parameters and state
type(tParameters), allocatable, dimension(:), private :: param type(tParameters), allocatable, dimension(:), private :: param
type(tDisloUCLAState ), allocatable, dimension(:), private :: & type(tDisloUCLAState), allocatable, dimension(:), private :: &
dotState, & dotState, &
state state
type(tDisloUCLAdependentState), allocatable, dimension(:), private :: dependentState type(tDisloUCLAdependentState), allocatable, dimension(:), private :: dependentState
@ -164,7 +164,6 @@ subroutine plastic_disloUCLA_init()
outputID outputID
character(len=pStringLen) :: & character(len=pStringLen) :: &
structure = '',&
extmsg = '' extmsg = ''
character(len=65536), dimension(:), allocatable :: & character(len=65536), dimension(:), allocatable :: &
outputs outputs
@ -197,8 +196,6 @@ subroutine plastic_disloUCLA_init()
dst => dependentState(phase_plasticityInstance(p)), & dst => dependentState(phase_plasticityInstance(p)), &
config => config_phase(p)) config => config_phase(p))
structure = config%getString('lattice_structure')
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! optional parameters that need to be defined ! optional parameters that need to be defined
prm%mu = lattice_mu(p) prm%mu = lattice_mu(p)
@ -213,36 +210,37 @@ subroutine plastic_disloUCLA_init()
prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray) prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray)
prm%totalNslip = sum(prm%Nslip) prm%totalNslip = sum(prm%Nslip)
slipActive: if (prm%totalNslip > 0_pInt) then slipActive: if (prm%totalNslip > 0_pInt) then
prm%Schmid = lattice_SchmidMatrix_slip(prm%Nslip,structure(1:3),& prm%Schmid = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal)) config%getFloat('c/a',defaultVal=0.0_pReal))
if(structure=='bcc') then
prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',& if(trim(config%getString('lattice_structure')) == 'bcc') then
prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',&
defaultVal = emptyRealArray) defaultVal = emptyRealArray)
prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt) prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1_pInt) prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1_pInt)
else else
prm%nonSchmid_pos = prm%Schmid prm%nonSchmid_pos = prm%Schmid
prm%nonSchmid_neg = prm%Schmid prm%nonSchmid_neg = prm%Schmid
endif endif
prm%interaction_SlipSlip = lattice_interaction_SlipSlip(prm%Nslip, & prm%interaction_SlipSlip = lattice_interaction_SlipSlip(prm%Nslip, &
config%getFloats('interaction_slipslip'), & config%getFloats('interaction_slipslip'), &
structure(1:3)) config%getString('lattice_structure'))
prm%rho0 = config%getFloats('rhoedge0', requiredShape=shape(prm%Nslip)) prm%rho0 = config%getFloats('rhoedge0', requiredSize=size(prm%Nslip))
prm%rhoDip0 = config%getFloats('rhoedgedip0', requiredShape=shape(prm%Nslip)) prm%rhoDip0 = config%getFloats('rhoedgedip0', requiredSize=size(prm%Nslip))
prm%v0 = config%getFloats('v0', requiredShape=shape(prm%Nslip)) prm%v0 = config%getFloats('v0', requiredSize=size(prm%Nslip))
prm%burgers = config%getFloats('slipburgers', requiredShape=shape(prm%Nslip)) prm%burgers = config%getFloats('slipburgers', requiredSize=size(prm%Nslip))
prm%H0kp = config%getFloats('qedge', requiredShape=shape(prm%Nslip)) prm%H0kp = config%getFloats('qedge', requiredSize=size(prm%Nslip))
prm%clambda = config%getFloats('clambdaslip', requiredShape=shape(prm%Nslip)) prm%clambda = config%getFloats('clambdaslip', requiredSize=size(prm%Nslip))
prm%tau_Peierls = config%getFloats('tau_peierls', requiredShape=shape(prm%Nslip)) ! ToDo: Deprecated prm%tau_Peierls = config%getFloats('tau_peierls', requiredSize=size(prm%Nslip)) ! ToDo: Deprecated
prm%p = config%getFloats('p_slip', requiredShape=shape(prm%Nslip), & prm%p = config%getFloats('p_slip', requiredSize=size(prm%Nslip), &
defaultVal=[(1.0_pReal,i=1_pInt,size(prm%Nslip))]) defaultVal=[(1.0_pReal,i=1_pInt,size(prm%Nslip))])
prm%q = config%getFloats('q_slip', requiredShape=shape(prm%Nslip), & prm%q = config%getFloats('q_slip', requiredSize=size(prm%Nslip), &
defaultVal=[(1.0_pReal,i=1_pInt,size(prm%Nslip))]) defaultVal=[(1.0_pReal,i=1_pInt,size(prm%Nslip))])
prm%kink_height = config%getFloats('kink_height', requiredShape=shape(prm%Nslip)) prm%kink_height = config%getFloats('kink_height', requiredSize=size(prm%Nslip))
prm%w = config%getFloats('kink_width', requiredShape=shape(prm%Nslip)) prm%w = config%getFloats('kink_width', requiredSize=size(prm%Nslip))
prm%omega = config%getFloats('omega', requiredShape=shape(prm%Nslip)) prm%omega = config%getFloats('omega', requiredSize=size(prm%Nslip))
prm%B = config%getFloats('friction_coeff', requiredShape=shape(prm%Nslip)) prm%B = config%getFloats('friction_coeff', requiredSize=size(prm%Nslip))
prm%SolidSolutionStrength = config%getFloat('solidsolutionstrength') ! ToDo: Deprecated prm%SolidSolutionStrength = config%getFloat('solidsolutionstrength') ! ToDo: Deprecated
prm%grainSize = config%getFloat('grainsize') prm%grainSize = config%getFloat('grainsize')
@ -250,7 +248,7 @@ subroutine plastic_disloUCLA_init()
prm%Qsd = config%getFloat('qsd') prm%Qsd = config%getFloat('qsd')
prm%atomicVolume = config%getFloat('catomicvolume') * prm%burgers**3.0_pReal prm%atomicVolume = config%getFloat('catomicvolume') * prm%burgers**3.0_pReal
prm%minDipDistance = config%getFloat('cedgedipmindistance') * prm%burgers prm%minDipDistance = config%getFloat('cedgedipmindistance') * prm%burgers
prm%dipoleformation = config%getFloat('dipoleformationfactor') > 0.0_pReal !should be on by default, ToDo: change to /key/-key prm%dipoleformation = config%getFloat('dipoleformationfactor') > 0.0_pReal !should be on by default, ToDo: change to /key/-type key
! expand: family => system ! expand: family => system
prm%rho0 = math_expand(prm%rho0, prm%Nslip) prm%rho0 = math_expand(prm%rho0, prm%Nslip)

38
src/plastic_kinematichardening.f90
View File

@ -151,7 +151,6 @@ subroutine plastic_kinehardening_init
outputID outputID
character(len=pStringLen) :: & character(len=pStringLen) :: &
structure = '',&
extmsg = '' extmsg = ''
character(len=65536), dimension(:), allocatable :: & character(len=65536), dimension(:), allocatable :: &
outputs outputs
@ -187,8 +186,6 @@ subroutine plastic_kinehardening_init
endif endif
#endif #endif
structure = config%getString('lattice_structure')
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! optional parameters that need to be defined ! optional parameters that need to be defined
prm%aTolResistance = config%getFloat('atol_resistance',defaultVal=1.0_pReal) prm%aTolResistance = config%getFloat('atol_resistance',defaultVal=1.0_pReal)
@ -203,28 +200,29 @@ subroutine plastic_kinehardening_init
prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray) prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray)
prm%totalNslip = sum(prm%Nslip) prm%totalNslip = sum(prm%Nslip)
slipActive: if (prm%totalNslip > 0_pInt) then slipActive: if (prm%totalNslip > 0_pInt) then
prm%Schmid = lattice_SchmidMatrix_slip(prm%Nslip,structure(1:3),& prm%Schmid = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal)) config%getFloat('c/a',defaultVal=0.0_pReal))
if(structure=='bcc') then
prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',& if(trim(config%getString('lattice_structure')) == 'bcc') then
defaultVal = emptyRealArray) prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',&
prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt) defaultVal = emptyRealArray)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1_pInt) prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1_pInt)
else else
prm%nonSchmid_pos = prm%Schmid prm%nonSchmid_pos = prm%Schmid
prm%nonSchmid_neg = prm%Schmid prm%nonSchmid_neg = prm%Schmid
endif endif
prm%interaction_SlipSlip = lattice_interaction_SlipSlip(prm%Nslip, & prm%interaction_SlipSlip = lattice_interaction_SlipSlip(prm%Nslip, &
config%getFloats('interaction_slipslip'), & config%getFloats('interaction_slipslip'), &
structure(1:3)) config%getString('lattice_structure'))
prm%crss0 = config%getFloats('crss0', requiredShape=shape(prm%Nslip)) prm%crss0 = config%getFloats('crss0', requiredSize=size(prm%Nslip))
prm%tau1 = config%getFloats('tau1', requiredShape=shape(prm%Nslip)) prm%tau1 = config%getFloats('tau1', requiredSize=size(prm%Nslip))
prm%tau1_b = config%getFloats('tau1_b', requiredShape=shape(prm%Nslip)) prm%tau1_b = config%getFloats('tau1_b', requiredSize=size(prm%Nslip))
prm%theta0 = config%getFloats('theta0', requiredShape=shape(prm%Nslip)) prm%theta0 = config%getFloats('theta0', requiredSize=size(prm%Nslip))
prm%theta1 = config%getFloats('theta1', requiredShape=shape(prm%Nslip)) prm%theta1 = config%getFloats('theta1', requiredSize=size(prm%Nslip))
prm%theta0_b = config%getFloats('theta0_b', requiredShape=shape(prm%Nslip)) prm%theta0_b = config%getFloats('theta0_b', requiredSize=size(prm%Nslip))
prm%theta1_b = config%getFloats('theta1_b', requiredShape=shape(prm%Nslip)) prm%theta1_b = config%getFloats('theta1_b', requiredSize=size(prm%Nslip))
prm%gdot0 = config%getFloat('gdot0') prm%gdot0 = config%getFloat('gdot0')
prm%n = config%getFloat('n_slip') prm%n = config%getFloat('n_slip')

48
src/plastic_phenopowerlaw.f90
View File

@ -153,7 +153,6 @@ subroutine plastic_phenopowerlaw_init
outputID outputID
character(len=pStringLen) :: & character(len=pStringLen) :: &
structure = '',&
extmsg = '' extmsg = ''
character(len=65536), dimension(:), allocatable :: & character(len=65536), dimension(:), allocatable :: &
outputs outputs
@ -181,8 +180,6 @@ subroutine plastic_phenopowerlaw_init
stt => state(phase_plasticityInstance(p)), & stt => state(phase_plasticityInstance(p)), &
config => config_phase(p)) config => config_phase(p))
structure = config%getString('lattice_structure')
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! optional parameters that need to be defined ! optional parameters that need to be defined
prm%twinB = config%getFloat('twin_b',defaultVal=1.0_pReal) prm%twinB = config%getFloat('twin_b',defaultVal=1.0_pReal)
@ -204,30 +201,31 @@ subroutine plastic_phenopowerlaw_init
prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray) prm%Nslip = config%getInts('nslip',defaultVal=emptyIntArray)
prm%totalNslip = sum(prm%Nslip) prm%totalNslip = sum(prm%Nslip)
slipActive: if (prm%totalNslip > 0_pInt) then slipActive: if (prm%totalNslip > 0_pInt) then
prm%Schmid_slip = lattice_SchmidMatrix_slip(prm%Nslip,structure(1:3),& prm%Schmid_slip = lattice_SchmidMatrix_slip(prm%Nslip,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal)) config%getFloat('c/a',defaultVal=0.0_pReal))
if(structure=='bcc') then
prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',& if(trim(config%getString('lattice_structure')) == 'bcc') then
prm%nonSchmidCoeff = config%getFloats('nonschmid_coefficients',&
defaultVal = emptyRealArray) defaultVal = emptyRealArray)
prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt) prm%nonSchmid_pos = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,+1_pInt)
prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1_pInt) prm%nonSchmid_neg = lattice_nonSchmidMatrix(prm%Nslip,prm%nonSchmidCoeff,-1_pInt)
else else
prm%nonSchmid_pos = prm%Schmid_slip prm%nonSchmid_pos = prm%Schmid_slip
prm%nonSchmid_neg = prm%Schmid_slip prm%nonSchmid_neg = prm%Schmid_slip
endif endif
prm%interaction_SlipSlip = lattice_interaction_SlipSlip(prm%Nslip, & prm%interaction_SlipSlip = lattice_interaction_SlipSlip(prm%Nslip, &
config%getFloats('interaction_slipslip'), & config%getFloats('interaction_slipslip'), &
structure(1:3)) config%getString('lattice_structure'))
prm%xi_slip_0 = config%getFloats('tau0_slip', requiredSize=size(prm%Nslip)) prm%xi_slip_0 = config%getFloats('tau0_slip', requiredSize=size(prm%Nslip))
prm%xi_slip_sat = config%getFloats('tausat_slip', requiredSize=size(prm%Nslip)) prm%xi_slip_sat = config%getFloats('tausat_slip', requiredSize=size(prm%Nslip))
prm%H_int = config%getFloats('h_int', requiredSize=size(prm%Nslip), & prm%H_int = config%getFloats('h_int', requiredSize=size(prm%Nslip), &
defaultVal=[(0.0_pReal,i=1_pInt,size(prm%Nslip))]) defaultVal=[(0.0_pReal,i=1_pInt,size(prm%Nslip))])
prm%gdot0_slip = config%getFloat('gdot0_slip') prm%gdot0_slip = config%getFloat('gdot0_slip')
prm%n_slip = config%getFloat('n_slip') prm%n_slip = config%getFloat('n_slip')
prm%a_slip = config%getFloat('a_slip') prm%a_slip = config%getFloat('a_slip')
prm%h0_SlipSlip = config%getFloat('h0_slipslip') prm%h0_SlipSlip = config%getFloat('h0_slipslip')
! expand: family => system ! expand: family => system
prm%xi_slip_0 = math_expand(prm%xi_slip_0, prm%Nslip) prm%xi_slip_0 = math_expand(prm%xi_slip_0, prm%Nslip)
@ -250,12 +248,12 @@ subroutine plastic_phenopowerlaw_init
prm%Ntwin = config%getInts('ntwin', defaultVal=emptyIntArray) prm%Ntwin = config%getInts('ntwin', defaultVal=emptyIntArray)
prm%totalNtwin = sum(prm%Ntwin) prm%totalNtwin = sum(prm%Ntwin)
twinActive: if (prm%totalNtwin > 0_pInt) then twinActive: if (prm%totalNtwin > 0_pInt) then
prm%Schmid_twin = lattice_SchmidMatrix_twin(prm%Ntwin,structure(1:3),& prm%Schmid_twin = lattice_SchmidMatrix_twin(prm%Ntwin,config%getString('lattice_structure'),&
config%getFloat('c/a',defaultVal=0.0_pReal)) config%getFloat('c/a',defaultVal=0.0_pReal))
prm%interaction_TwinTwin = lattice_interaction_TwinTwin(prm%Ntwin,& prm%interaction_TwinTwin = lattice_interaction_TwinTwin(prm%Ntwin,&
config%getFloats('interaction_twintwin'), & config%getFloats('interaction_twintwin'), &
structure(1:3)) config%getString('lattice_structure'))
prm%gamma_twin_char = lattice_characteristicShear_twin(prm%Ntwin,structure(1:3),& prm%gamma_twin_char = lattice_characteristicShear_twin(prm%Ntwin,config%getString('lattice_structure'),&
config%getFloat('c/a')) config%getFloat('c/a'))
prm%xi_twin_0 = config%getFloats('tau0_twin',requiredSize=size(prm%Ntwin)) prm%xi_twin_0 = config%getFloats('tau0_twin',requiredSize=size(prm%Ntwin))
@ -282,10 +280,10 @@ subroutine plastic_phenopowerlaw_init
slipAndTwinActive: if (prm%totalNslip > 0_pInt .and. prm%totalNtwin > 0_pInt) then slipAndTwinActive: if (prm%totalNslip > 0_pInt .and. prm%totalNtwin > 0_pInt) then
prm%interaction_SlipTwin = lattice_interaction_SlipTwin(prm%Nslip,prm%Ntwin,& prm%interaction_SlipTwin = lattice_interaction_SlipTwin(prm%Nslip,prm%Ntwin,&
config%getFloats('interaction_sliptwin'), & config%getFloats('interaction_sliptwin'), &
structure(1:3)) config%getString('lattice_structure'))
prm%interaction_TwinSlip = lattice_interaction_TwinSlip(prm%Ntwin,prm%Nslip,& prm%interaction_TwinSlip = lattice_interaction_TwinSlip(prm%Ntwin,prm%Nslip,&
config%getFloats('interaction_twinslip'), & config%getFloats('interaction_twinslip'), &
structure(1:3)) config%getString('lattice_structure'))
else slipAndTwinActive else slipAndTwinActive
allocate(prm%interaction_SlipTwin(prm%totalNslip,prm%TotalNtwin)) ! at least one dimension is 0 allocate(prm%interaction_SlipTwin(prm%totalNslip,prm%TotalNtwin)) ! at least one dimension is 0
allocate(prm%interaction_TwinSlip(prm%totalNtwin,prm%TotalNslip)) ! at least one dimension is 0 allocate(prm%interaction_TwinSlip(prm%totalNtwin,prm%TotalNslip)) ! at least one dimension is 0