moved reading in of lattice type and elastic constants to lattice module

removed structure type for hex, fcc, bcc, now defining slip/twin systems for each phase found in material.config
constitutive modules will only be initialized if needed
homogenizedC function is only needed for models incorporating twinning in a physical way (titanmod and dislotwin)
This commit is contained in:
Martin Diehl 2014-03-08 20:50:31 +00:00
parent 3b20f1154e
commit 2b589c3d71
9 changed files with 2510 additions and 3113 deletions

View File

@ -154,12 +154,12 @@ subroutine constitutive_init
! parse plasticities from config file
if (.not. IO_open_jobFile_stat(FILEUNIT,material_localFileExt)) & ! no local material configuration present...
call IO_open_file(FILEUNIT,material_configFile) ! ... open material.config file
call constitutive_none_init(FILEUNIT)
call constitutive_j2_init(FILEUNIT)
call constitutive_phenopowerlaw_init(FILEUNIT)
call constitutive_titanmod_init(FILEUNIT)
call constitutive_dislotwin_init(FILEUNIT)
call constitutive_nonlocal_init(FILEUNIT)
if (any(phase_plasticity == PLASTICITY_NONE_ID)) call constitutive_none_init(FILEUNIT)
if (any(phase_plasticity == PLASTICITY_J2_ID)) call constitutive_j2_init(FILEUNIT)
if (any(phase_plasticity == PLASTICITY_PHENOPOWERLAW_ID)) call constitutive_phenopowerlaw_init(FILEUNIT)
if (any(phase_plasticity == PLASTICITY_TITANMOD_ID)) call constitutive_titanmod_init(FILEUNIT)
if (any(phase_plasticity == PLASTICITY_DISLOTWIN_ID)) call constitutive_dislotwin_init(FILEUNIT)
if (any(phase_plasticity == PLASTICITY_NONLOCAL_ID)) call constitutive_nonlocal_init(FILEUNIT)
close(FILEUNIT)
write(6,'(/,a)') ' <<<+- constitutive init -+>>>'
@ -352,7 +352,7 @@ subroutine constitutive_init
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_dislotwin_sizeDotState(instance)))
enddo
endif
constitutive_state0(g,i,e)%p = constitutive_dislotwin_stateInit(instance)
constitutive_state0(g,i,e)%p = constitutive_dislotwin_stateInit(instance,material_phase(g,i,e))
constitutive_aTolState(g,i,e)%p = constitutive_dislotwin_aTolState(instance)
constitutive_sizeState(g,i,e) = constitutive_dislotwin_sizeState(instance)
constitutive_sizeDotState(g,i,e) = constitutive_dislotwin_sizeDotState(instance)
@ -379,11 +379,11 @@ subroutine constitutive_init
allocate(constitutive_RKCK45dotState(s,g,i,e)%p(constitutive_titanmod_sizeDotState(instance)))
enddo
endif
constitutive_state0(g,i,e)%p = constitutive_titanmod_stateInit(instance)
constitutive_aTolState(g,i,e)%p = constitutive_titanmod_aTolState(instance)
constitutive_sizeState(g,i,e) = constitutive_titanmod_sizeState(instance)
constitutive_sizeDotState(g,i,e) = constitutive_titanmod_sizeDotState(instance)
constitutive_sizePostResults(g,i,e) = constitutive_titanmod_sizePostResults(instance)
constitutive_state0(g,i,e)%p = constitutive_titanmod_stateInit(instance,material_phase(g,i,e))
constitutive_aTolState(g,i,e)%p = constitutive_titanmod_aTolState(instance)
constitutive_sizeState(g,i,e) = constitutive_titanmod_sizeState(instance)
constitutive_sizeDotState(g,i,e) = constitutive_titanmod_sizeDotState(instance)
constitutive_sizePostResults(g,i,e) = constitutive_titanmod_sizePostResults(instance)
case (PLASTICITY_NONLOCAL_ID)
nonlocalConstitutionPresent = .true.
if(myNgrains/=1_pInt) call IO_error(252_pInt, e,i,g)
@ -465,25 +465,15 @@ pure function constitutive_homogenizedC(ipc,ip,el)
use material, only: &
phase_plasticity, &
material_phase, &
PLASTICITY_NONE_ID, &
PLASTICITY_J2_ID, &
PLASTICITY_PHENOPOWERLAW_ID, &
PLASTICITY_DISLOTWIN_ID, &
PLASTICITY_TITANMOD_ID, &
PLASTICITY_NONLOCAL_ID
use constitutive_none, only: &
constitutive_none_homogenizedC
use constitutive_j2, only: &
constitutive_j2_homogenizedC
use constitutive_phenopowerlaw, only: &
constitutive_phenopowerlaw_homogenizedC
PLASTICITY_DISLOTWIN_ID
use constitutive_dislotwin, only: &
constitutive_dislotwin_homogenizedC
use constitutive_titanmod, only: &
constitutive_titanmod_homogenizedC
use constitutive_nonlocal, only: &
constitutive_nonlocal_homogenizedC
use lattice, only: &
lattice_C66
implicit none
real(pReal), dimension(6,6) :: constitutive_homogenizedC
integer(pInt), intent(in) :: &
@ -493,23 +483,14 @@ pure function constitutive_homogenizedC(ipc,ip,el)
select case (phase_plasticity(material_phase(ipc,ip,el)))
case (PLASTICITY_NONE_ID)
constitutive_homogenizedC = constitutive_none_homogenizedC(ipc,ip,el)
case (PLASTICITY_J2_ID)
constitutive_homogenizedC = constitutive_j2_homogenizedC(ipc,ip,el)
case (PLASTICITY_PHENOPOWERLAW_ID)
constitutive_homogenizedC = constitutive_phenopowerlaw_homogenizedC(ipc,ip,el)
case (PLASTICITY_DISLOTWIN_ID)
constitutive_homogenizedC = constitutive_dislotwin_homogenizedC(constitutive_state,ipc,ip,el)
case (PLASTICITY_TITANMOD_ID)
constitutive_homogenizedC = constitutive_titanmod_homogenizedC(constitutive_state,ipc,ip,el)
case (PLASTICITY_NONLOCAL_ID)
constitutive_homogenizedC = constitutive_nonlocal_homogenizedC(ipc,ip,el)
case default
constitutive_homogenizedC = lattice_C66(1:6,1:6,material_phase(ipc,ip,el))
end select

File diff suppressed because it is too large Load Diff

View File

@ -66,8 +66,6 @@ module constitutive_j2
constitutive_j2_tausat_SinhFitC, & !< fitting parameter for normalized strain rate vs. stress function
constitutive_j2_tausat_SinhFitD !< fitting parameter for normalized strain rate vs. stress function
real(pReal), dimension(:,:,:), allocatable, private :: &
constitutive_j2_Cslip_66
enum, bind(c)
enumerator :: undefined_ID, &
flowstress_ID, &
@ -80,7 +78,6 @@ module constitutive_j2
constitutive_j2_init, &
constitutive_j2_stateInit, &
constitutive_j2_aTolState, &
constitutive_j2_homogenizedC, &
constitutive_j2_LpAndItsTangent, &
constitutive_j2_dotState, &
constitutive_j2_postResults
@ -128,7 +125,7 @@ subroutine constitutive_j2_init(fileUnit)
integer(pInt), parameter :: MAXNCHUNKS = 7_pInt
integer(pInt), dimension(1_pInt+2_pInt*MAXNCHUNKS) :: positions
integer(pInt) :: section = 0_pInt, maxNinstance, instance,o, mySize
integer(pInt) :: phase, maxNinstance, instance,o, mySize
character(len=65536) :: &
tag = '', &
line = ''
@ -152,7 +149,6 @@ subroutine constitutive_j2_init(fileUnit)
constitutive_j2_output = ''
allocate(constitutive_j2_outputID(maxval(phase_Noutput),maxNinstance), source=undefined_ID)
allocate(constitutive_j2_Noutput(maxNinstance), source=0_pInt)
allocate(constitutive_j2_Cslip_66(6,6,maxNinstance), source=0.0_pReal)
allocate(constitutive_j2_fTaylor(maxNinstance), source=0.0_pReal)
allocate(constitutive_j2_tau0(maxNinstance), source=0.0_pReal)
allocate(constitutive_j2_gdot0(maxNinstance), source=0.0_pReal)
@ -168,11 +164,12 @@ subroutine constitutive_j2_init(fileUnit)
allocate(constitutive_j2_tausat_SinhFitD(maxNinstance), source=0.0_pReal)
rewind(fileUnit)
phase = 0_pInt
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partPhase) ! wind forward to <phase>
line = IO_read(fileUnit)
enddo
do while (trim(line) /= IO_EOF) ! read through sections of phase part
parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') then ! stop at next part
@ -180,19 +177,19 @@ subroutine constitutive_j2_init(fileUnit)
exit
endif
if (IO_getTag(line,'[',']') /= '') then ! next section
section = section + 1_pInt ! advance section counter
if (phase_plasticity(section) == PLASTICITY_J2_ID) then
instance = phase_plasticityInstance(section)
constitutive_j2_Cslip_66(1:6,1:6,instance) = lattice_Cslip_66(1:6,1:6,section)
phase = phase + 1_pInt ! advance section counter
if (phase_plasticity(phase) == PLASTICITY_J2_ID) then
instance = phase_plasticityInstance(phase)
endif
cycle ! skip to next line
endif
if (section > 0_pInt ) then; if (phase_plasticity(section) == PLASTICITY_J2_ID) then ! one of my sections. Do not short-circuit here (.and. between if-statements), it's not safe in Fortran
instance = phase_plasticityInstance(section) ! which instance of my plasticity is present phase
if (phase > 0_pInt ) then; if (phase_plasticity(phase) == PLASTICITY_J2_ID) then ! one of my sections. Do not short-circuit here (.and. between if-statements), it's not safe in Fortran
instance = phase_plasticityInstance(phase) ! which instance of my plasticity is present phase
positions = IO_stringPos(line,MAXNCHUNKS)
tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key
select case(tag)
case ('plasticity','elasticity','lattice_structure','covera_ratio',&
case ('plasticity','elasticity','lattice_structure', &
'covera_ratio','c/a_ratio','c/a', &
'c11','c12','c13','c22','c23','c33','c44','c55','c66')
case ('(output)')
constitutive_j2_Noutput(instance) = constitutive_j2_Noutput(instance) + 1_pInt
@ -250,7 +247,7 @@ subroutine constitutive_j2_init(fileUnit)
call IO_error(210_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_J2_label//')')
end select
endif; endif
enddo
enddo parsingFile
instancesLoop: do instance = 1_pInt,maxNinstance
outputsLoop: do o = 1_pInt,constitutive_j2_Noutput(instance)
@ -292,42 +289,15 @@ end function constitutive_j2_stateInit
pure function constitutive_j2_aTolState(instance)
implicit none
integer(pInt), intent(in) :: instance !< number specifying the instance of the plasticity
real(pReal), dimension(1) :: constitutive_j2_aTolState
integer(pInt), intent(in) :: instance !< number specifying the instance of the plasticity
real(pReal), dimension(constitutive_j2_sizeState(instance)) :: &
constitutive_j2_aTolState
constitutive_j2_aTolState = constitutive_j2_aTolResistance(instance)
end function constitutive_j2_aTolState
!--------------------------------------------------------------------------------------------------
!> @brief returns the homogenized elasticity matrix
!--------------------------------------------------------------------------------------------------
pure function constitutive_j2_homogenizedC(ipc,ip,el)
use mesh, only: &
mesh_NcpElems, &
mesh_maxNips
use material, only: &
homogenization_maxNgrains,&
material_phase, &
phase_plasticityInstance
implicit none
real(pReal), dimension(6,6) :: &
constitutive_j2_homogenizedC
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
constitutive_j2_homogenizedC = constitutive_j2_Cslip_66(1:6,1:6,&
phase_plasticityInstance(material_phase(ipc,ip,el)))
end function constitutive_j2_homogenizedC
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!--------------------------------------------------------------------------------------------------

View File

@ -25,11 +25,8 @@
!--------------------------------------------------------------------------------------------------
module constitutive_none
use prec, only: &
pReal, &
pInt
use lattice, only: &
LATTICE_undefined_ID
implicit none
private
integer(pInt), dimension(:), allocatable, public, protected :: &
@ -40,15 +37,8 @@ module constitutive_none
integer(pInt), dimension(:,:), allocatable, target, public :: &
constitutive_none_sizePostResult !< size of each post result output
integer(kind(LATTICE_undefined_ID)), dimension(:), allocatable, public :: &
constitutive_none_structureID !< ID of the lattice structure
real(pReal), dimension(:,:,:), allocatable, private :: &
constitutive_none_Cslip_66
public :: &
constitutive_none_init, &
constitutive_none_homogenizedC
constitutive_none_init
contains
@ -63,43 +53,21 @@ subroutine constitutive_none_init(fileUnit)
debug_level, &
debug_constitutive, &
debug_levelBasic
use math, only: &
math_Mandel3333to66, &
math_Voigt66to3333
use IO, only: &
IO_read, &
IO_lc, &
IO_getTag, &
IO_isBlank, &
IO_stringPos, &
IO_stringValue, &
IO_floatValue, &
IO_error, &
IO_timeStamp, &
IO_EOF
IO_timeStamp
use material, only: &
homogenization_maxNgrains, &
phase_plasticity, &
phase_plasticityInstance, &
phase_Noutput, &
PLASTICITY_NONE_label, &
PLASTICITY_NONE_ID, &
MATERIAL_partPhase
use lattice
implicit none
integer(pInt), intent(in) :: fileUnit
integer(pInt), parameter :: MAXNCHUNKS = 7_pInt
integer(pInt), dimension(1_pInt+2_pInt*MAXNCHUNKS) :: positions
integer(pInt) :: section = 0_pInt, maxNinstance, instance
character(len=32) :: &
structure = ''
character(len=65536) :: &
tag = '', &
line = ''
integer(pInt), intent(in) :: fileUnit
integer(pInt) :: maxNinstance
write(6,'(/,a)') ' <<<+- constitutive_'//PLASTICITY_NONE_label//' init -+>>>'
write(6,'(a)') ' $Id$'
@ -115,106 +83,7 @@ subroutine constitutive_none_init(fileUnit)
allocate(constitutive_none_sizeDotState(maxNinstance), source=1_pInt)
allocate(constitutive_none_sizeState(maxNinstance), source=1_pInt)
allocate(constitutive_none_sizePostResults(maxNinstance), source=0_pInt)
allocate(constitutive_none_structureID(maxNinstance), source=LATTICE_undefined_ID)
allocate(constitutive_none_Cslip_66(6,6,maxNinstance), source=0.0_pReal)
rewind(fileUnit)
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partPhase) ! wind forward to <phase>
line = IO_read(fileUnit)
enddo
do while (trim(line) /= IO_EOF) ! read through sections of phase part
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') then ! stop at next part
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif
if (IO_getTag(line,'[',']') /= '') then ! next section
section = section + 1_pInt ! advance section counter
cycle ! skip to next line
endif
if (section > 0_pInt ) then ! do not short-circuit here (.and. with next if-statement). It's not safe in Fortran
if (phase_plasticity(section) == PLASTICITY_NONE_ID) then ! one of my sections
instance = phase_plasticityInstance(section) ! which instance of my plasticity is present phase
positions = IO_stringPos(line,MAXNCHUNKS)
tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key
select case(tag)
case ('plasticity','elasticity','covera_ratio')
case ('lattice_structure')
structure = IO_lc(IO_stringValue(line,positions,2_pInt))
select case(structure(1:3))
case(LATTICE_iso_label)
constitutive_none_structureID(instance) = LATTICE_iso_ID
case(LATTICE_fcc_label)
constitutive_none_structureID(instance) = LATTICE_fcc_ID
case(LATTICE_bcc_label)
constitutive_none_structureID(instance) = LATTICE_bcc_ID
case(LATTICE_hex_label)
constitutive_none_structureID(instance) = LATTICE_hex_ID
case(LATTICE_ort_label)
constitutive_none_structureID(instance) = LATTICE_ort_ID
end select
case ('c11')
constitutive_none_Cslip_66(1,1,instance) = IO_floatValue(line,positions,2_pInt)
case ('c12')
constitutive_none_Cslip_66(1,2,instance) = IO_floatValue(line,positions,2_pInt)
case ('c13')
constitutive_none_Cslip_66(1,3,instance) = IO_floatValue(line,positions,2_pInt)
case ('c22')
constitutive_none_Cslip_66(2,2,instance) = IO_floatValue(line,positions,2_pInt)
case ('c23')
constitutive_none_Cslip_66(2,3,instance) = IO_floatValue(line,positions,2_pInt)
case ('c33')
constitutive_none_Cslip_66(3,3,instance) = IO_floatValue(line,positions,2_pInt)
case ('c44')
constitutive_none_Cslip_66(4,4,instance) = IO_floatValue(line,positions,2_pInt)
case ('c55')
constitutive_none_Cslip_66(5,5,instance) = IO_floatValue(line,positions,2_pInt)
case ('c66')
constitutive_none_Cslip_66(6,6,instance) = IO_floatValue(line,positions,2_pInt)
case default
call IO_error(210_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_NONE_label//')')
end select
endif
endif
enddo
instancesLoop: do instance = 1_pInt,maxNinstance
constitutive_none_Cslip_66(1:6,1:6,instance) = &
lattice_symmetrizeC66(constitutive_none_structureID(instance),constitutive_none_Cslip_66(1:6,1:6,instance))
constitutive_none_Cslip_66(1:6,1:6,instance) = & ! Literature data is Voigt, DAMASK uses Mandel
math_Mandel3333to66(math_Voigt66to3333(constitutive_none_Cslip_66(1:6,1:6,instance)))
enddo instancesLoop
end subroutine constitutive_none_init
!--------------------------------------------------------------------------------------------------
!> @brief returns the homogenized elasticity matrix
!--------------------------------------------------------------------------------------------------
pure function constitutive_none_homogenizedC(ipc,ip,el)
use prec, only: &
p_vec
use mesh, only: &
mesh_NcpElems, &
mesh_maxNips
use material, only: &
homogenization_maxNgrains, &
material_phase, &
phase_plasticityInstance
implicit none
real(pReal), dimension(6,6) :: &
constitutive_none_homogenizedC
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
constitutive_none_homogenizedC = constitutive_none_Cslip_66(1:6,1:6,&
phase_plasticityInstance(material_phase(ipc,ip,el)))
end function constitutive_none_homogenizedC
end module constitutive_none

File diff suppressed because it is too large Load Diff

View File

@ -34,8 +34,7 @@ module constitutive_phenopowerlaw
integer(pInt), dimension(:), allocatable, public, protected :: &
constitutive_phenopowerlaw_sizeDotState, &
constitutive_phenopowerlaw_sizeState, &
constitutive_phenopowerlaw_sizePostResults, & !< cumulative size of post results
constitutive_phenopowerlaw_structure
constitutive_phenopowerlaw_sizePostResults !< cumulative size of post results
integer(pInt), dimension(:,:), allocatable, target, public :: &
constitutive_phenopowerlaw_sizePostResult !< size of each post result output
@ -87,8 +86,8 @@ module constitutive_phenopowerlaw
constitutive_phenopowerlaw_hardeningMatrix_SlipSlip, &
constitutive_phenopowerlaw_hardeningMatrix_SlipTwin, &
constitutive_phenopowerlaw_hardeningMatrix_TwinSlip, &
constitutive_phenopowerlaw_hardeningMatrix_TwinTwin, &
constitutive_phenopowerlaw_Cslip_66
constitutive_phenopowerlaw_hardeningMatrix_TwinTwin
enum, bind(c)
enumerator :: undefined_ID, &
resistance_slip_ID, &
@ -109,7 +108,6 @@ module constitutive_phenopowerlaw
constitutive_phenopowerlaw_init, &
constitutive_phenopowerlaw_stateInit, &
constitutive_phenopowerlaw_aTolState, &
constitutive_phenopowerlaw_homogenizedC, &
constitutive_phenopowerlaw_LpAndItsTangent, &
constitutive_phenopowerlaw_dotState, &
constitutive_phenopowerlaw_postResults
@ -160,16 +158,13 @@ subroutine constitutive_phenopowerlaw_init(fileUnit)
integer(pInt), parameter :: MAXNCHUNKS = LATTICE_maxNinteraction + 1_pInt
integer(pInt), dimension(1_pInt+2_pInt*MAXNCHUNKS) :: positions
integer(pInt), dimension(7) :: configNchunks
integer(pInt) :: &
maxNinstance, &
instance,j,k, f,o, &
instance,phase,j,k, f,o, &
Nchunks_SlipSlip, Nchunks_SlipTwin, Nchunks_TwinSlip, Nchunks_TwinTwin, &
Nchunks_SlipFamilies, Nchunks_TwinFamilies, Nchunks_nonSchmid, &
structID, index_myFamily, index_otherFamily, &
mySize=0_pInt, section = 0_pInt
character(len=32) :: &
structure = ''
index_myFamily, index_otherFamily, &
mySize=0_pInt
character(len=65536) :: &
tag = '', &
line = ''
@ -184,15 +179,7 @@ subroutine constitutive_phenopowerlaw_init(fileUnit)
if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
Nchunks_SlipFamilies = lattice_maxNslipFamily
Nchunks_TwinFamilies = lattice_maxNtwinFamily
Nchunks_SlipSlip = lattice_maxNinteraction
Nchunks_SlipTwin = lattice_maxNinteraction
Nchunks_TwinSlip = lattice_maxNinteraction
Nchunks_TwinTwin = lattice_maxNinteraction
Nchunks_nonSchmid = lattice_maxNnonSchmid
allocate(constitutive_phenopowerlaw_sizeDotState(maxNinstance), source=0_pInt)
allocate(constitutive_phenopowerlaw_sizeState(maxNinstance), source=0_pInt)
allocate(constitutive_phenopowerlaw_sizePostResults(maxNinstance), source=0_pInt)
@ -202,12 +189,10 @@ subroutine constitutive_phenopowerlaw_init(fileUnit)
constitutive_phenopowerlaw_output = ''
allocate(constitutive_phenopowerlaw_outputID(maxval(phase_Noutput),maxNinstance),source=undefined_ID)
allocate(constitutive_phenopowerlaw_Noutput(maxNinstance), source=0_pInt)
allocate(constitutive_phenopowerlaw_structure(maxNinstance), source=0_pInt)
allocate(constitutive_phenopowerlaw_Nslip(lattice_maxNslipFamily,maxNinstance), source=0_pInt)
allocate(constitutive_phenopowerlaw_Ntwin(lattice_maxNtwinFamily,maxNinstance), source=0_pInt)
allocate(constitutive_phenopowerlaw_totalNslip(maxNinstance), source=0_pInt)
allocate(constitutive_phenopowerlaw_totalNtwin(maxNinstance), source=0_pInt)
allocate(constitutive_phenopowerlaw_Cslip_66(6,6,maxNinstance), source=0.0_pReal)
allocate(constitutive_phenopowerlaw_gdot0_slip(maxNinstance), source=0.0_pReal)
allocate(constitutive_phenopowerlaw_n_slip(maxNinstance), source=0.0_pReal)
allocate(constitutive_phenopowerlaw_tau0_slip(lattice_maxNslipFamily,maxNinstance), &
@ -243,40 +228,38 @@ subroutine constitutive_phenopowerlaw_init(fileUnit)
source=0.0_pReal)
rewind(fileUnit)
phase = 0_pInt
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partPhase) ! wind forward to <phase>
line = IO_read(fileUnit)
enddo
do while (trim(line) /= IO_EOF) ! read through sections of phase part
parsingFile: do while (trim(line) /= IO_EOF) ! read through sections of phase part
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') then ! stop at next part
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif
if (IO_getTag(line,'[',']') /= '') then ! next section
section = section + 1_pInt ! advance section counter
if (phase_plasticity(section) == PLASTICITY_PHENOPOWERLAW_ID) then
instance = phase_plasticityInstance(section)
constitutive_phenopowerlaw_Cslip_66(1:6,1:6,instance) = lattice_Cslip_66(1:6,1:6,section)
constitutive_phenopowerlaw_structure(instance) = lattice_structure(section)
configNchunks = lattice_configNchunks(lattice_structureID(section))
Nchunks_SlipFamilies = configNchunks(1)
Nchunks_TwinFamilies = configNchunks(2)
Nchunks_SlipSlip = configNchunks(3)
Nchunks_SlipTwin = configNchunks(4)
Nchunks_TwinSlip = configNchunks(5)
Nchunks_TwinTwin = configNchunks(6)
Nchunks_nonSchmid = configNchunks(7)
if (IO_getTag(line,'[',']') /= '') then ! next phase
phase = phase + 1_pInt ! advance phase section counter
if (phase_plasticity(phase) == PLASTICITY_PHENOPOWERLAW_ID) then
Nchunks_SlipFamilies = count(lattice_NslipSystem(:,phase) > 0_pInt)
Nchunks_TwinFamilies = count(lattice_NtwinSystem(:,phase) > 0_pInt)
Nchunks_SlipSlip = maxval(lattice_interactionSlipSlip(:,:,phase))
Nchunks_SlipTwin = maxval(lattice_interactionSlipTwin(:,:,phase))
Nchunks_TwinSlip = maxval(lattice_interactionTwinSlip(:,:,phase))
Nchunks_TwinTwin = maxval(lattice_interactionTwinTwin(:,:,phase))
Nchunks_nonSchmid = lattice_NnonSchmid(phase)
endif
cycle ! skip to next line
endif
if (section > 0_pInt ) then; if (phase_plasticity(section) == PLASTICITY_PHENOPOWERLAW_ID) then ! one of my sections. Do not short-circuit here (.and. between if-statements), it's not safe in Fortran
instance = phase_plasticityInstance(section) ! which instance of my plasticity is present phase
if (phase > 0_pInt ) then; if (phase_plasticity(phase) == PLASTICITY_PHENOPOWERLAW_ID) then ! one of my phases. Do not short-circuit here (.and. between if-statements), it's not safe in Fortran
instance = phase_plasticityInstance(phase) ! which instance of my plasticity is present phase
positions = IO_stringPos(line,MAXNCHUNKS)
tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key
select case(tag)
case ('plasticity','elasticity','lattice_structure','covera_ratio',&
case ('plasticity','elasticity','lattice_structure', &
'covera_ratio','c/a_ratio','c/a', &
'c11','c12','c13','c22','c23','c33','c44','c55','c66')
case ('(output)')
constitutive_phenopowerlaw_Noutput(instance) = constitutive_phenopowerlaw_Noutput(instance) + 1_pInt
@ -401,47 +384,49 @@ subroutine constitutive_phenopowerlaw_init(fileUnit)
call IO_error(210_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_PHENOPOWERLAW_label//')')
end select
endif; endif
enddo
enddo parsingFile
sanityChecks: do instance = 1_pInt,maxNinstance
constitutive_phenopowerlaw_Nslip(1:lattice_maxNslipFamily,instance) = &
min(lattice_NslipSystem(1:lattice_maxNslipFamily,constitutive_phenopowerlaw_structure(instance)),& ! limit active slip systems per family to min of available and requested
sanityChecks: do phase = 1_pInt, size(phase_plasticity)
myPhase: if (phase_plasticity(phase) == PLASTICITY_phenopowerlaw_ID) then
instance = phase_plasticityInstance(phase)
constitutive_phenopowerlaw_Nslip(1:lattice_maxNslipFamily,instance) = &
min(lattice_NslipSystem(1:lattice_maxNslipFamily,phase),& ! limit active slip systems per family to min of available and requested
constitutive_phenopowerlaw_Nslip(1:lattice_maxNslipFamily,instance))
constitutive_phenopowerlaw_Ntwin(1:lattice_maxNtwinFamily,instance) = &
min(lattice_NtwinSystem(1:lattice_maxNtwinFamily,constitutive_phenopowerlaw_structure(instance)),& ! limit active twin systems per family to min of available and requested
constitutive_phenopowerlaw_Ntwin(1:lattice_maxNtwinFamily,instance) = &
min(lattice_NtwinSystem(1:lattice_maxNtwinFamily,phase),& ! limit active twin systems per family to min of available and requested
constitutive_phenopowerlaw_Ntwin(:,instance))
constitutive_phenopowerlaw_totalNslip(instance) = sum(constitutive_phenopowerlaw_Nslip(:,instance)) ! how many slip systems altogether
constitutive_phenopowerlaw_totalNtwin(instance) = sum(constitutive_phenopowerlaw_Ntwin(:,instance)) ! how many twin systems altogether
if (any(constitutive_phenopowerlaw_tau0_slip(:,instance) < 0.0_pReal .and. &
constitutive_phenopowerlaw_Nslip(:,instance) > 0)) &
call IO_error(211_pInt,el=instance,ext_msg='tau0_slip ('//PLASTICITY_PHENOPOWERLAW_label//')')
if (constitutive_phenopowerlaw_gdot0_slip(instance) <= 0.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='gdot0_slip ('//PLASTICITY_PHENOPOWERLAW_label//')')
if (constitutive_phenopowerlaw_n_slip(instance) <= 0.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='n_slip ('//PLASTICITY_PHENOPOWERLAW_label//')')
if (any(constitutive_phenopowerlaw_tausat_slip(:,instance) <= 0.0_pReal .and. &
constitutive_phenopowerlaw_Nslip(:,instance) > 0)) &
call IO_error(211_pInt,el=instance,ext_msg='tausat_slip ('//PLASTICITY_PHENOPOWERLAW_label//')')
if (any(constitutive_phenopowerlaw_a_slip(instance) == 0.0_pReal .and. &
constitutive_phenopowerlaw_Nslip(:,instance) > 0)) &
call IO_error(211_pInt,el=instance,ext_msg='a_slip ('//PLASTICITY_PHENOPOWERLAW_label//')')
if (any(constitutive_phenopowerlaw_tau0_twin(:,instance) < 0.0_pReal .and. &
constitutive_phenopowerlaw_Ntwin(:,instance) > 0)) &
call IO_error(211_pInt,el=instance,ext_msg='tau0_twin ('//PLASTICITY_PHENOPOWERLAW_label//')')
if ( constitutive_phenopowerlaw_gdot0_twin(instance) <= 0.0_pReal .and. &
any(constitutive_phenopowerlaw_Ntwin(:,instance) > 0)) &
call IO_error(211_pInt,el=instance,ext_msg='gdot0_twin ('//PLASTICITY_PHENOPOWERLAW_label//')')
if ( constitutive_phenopowerlaw_n_twin(instance) <= 0.0_pReal .and. &
any(constitutive_phenopowerlaw_Ntwin(:,instance) > 0)) &
call IO_error(211_pInt,el=instance,ext_msg='n_twin ('//PLASTICITY_PHENOPOWERLAW_label//')')
if (constitutive_phenopowerlaw_aTolResistance(instance) <= 0.0_pReal) &
constitutive_phenopowerlaw_aTolResistance(instance) = 1.0_pReal ! default absolute tolerance 1 Pa
if (constitutive_phenopowerlaw_aTolShear(instance) <= 0.0_pReal) &
constitutive_phenopowerlaw_aTolShear(instance) = 1.0e-6_pReal ! default absolute tolerance 1e-6
if (constitutive_phenopowerlaw_aTolTwinfrac(instance) <= 0.0_pReal) &
constitutive_phenopowerlaw_aTolTwinfrac(instance) = 1.0e-6_pReal ! default absolute tolerance 1e-6
constitutive_phenopowerlaw_totalNslip(instance) = sum(constitutive_phenopowerlaw_Nslip(:,instance)) ! how many slip systems altogether
constitutive_phenopowerlaw_totalNtwin(instance) = sum(constitutive_phenopowerlaw_Ntwin(:,instance)) ! how many twin systems altogether
if (any(constitutive_phenopowerlaw_tau0_slip(:,instance) < 0.0_pReal .and. &
constitutive_phenopowerlaw_Nslip(:,instance) > 0)) &
call IO_error(211_pInt,el=instance,ext_msg='tau0_slip ('//PLASTICITY_PHENOPOWERLAW_label//')')
if (constitutive_phenopowerlaw_gdot0_slip(instance) <= 0.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='gdot0_slip ('//PLASTICITY_PHENOPOWERLAW_label//')')
if (constitutive_phenopowerlaw_n_slip(instance) <= 0.0_pReal) &
call IO_error(211_pInt,el=instance,ext_msg='n_slip ('//PLASTICITY_PHENOPOWERLAW_label//')')
if (any(constitutive_phenopowerlaw_tausat_slip(:,instance) <= 0.0_pReal .and. &
constitutive_phenopowerlaw_Nslip(:,instance) > 0)) &
call IO_error(211_pInt,el=instance,ext_msg='tausat_slip ('//PLASTICITY_PHENOPOWERLAW_label//')')
if (any(constitutive_phenopowerlaw_a_slip(instance) == 0.0_pReal .and. &
constitutive_phenopowerlaw_Nslip(:,instance) > 0)) &
call IO_error(211_pInt,el=instance,ext_msg='a_slip ('//PLASTICITY_PHENOPOWERLAW_label//')')
if (any(constitutive_phenopowerlaw_tau0_twin(:,instance) < 0.0_pReal .and. &
constitutive_phenopowerlaw_Ntwin(:,instance) > 0)) &
call IO_error(211_pInt,el=instance,ext_msg='tau0_twin ('//PLASTICITY_PHENOPOWERLAW_label//')')
if ( constitutive_phenopowerlaw_gdot0_twin(instance) <= 0.0_pReal .and. &
any(constitutive_phenopowerlaw_Ntwin(:,instance) > 0)) &
call IO_error(211_pInt,el=instance,ext_msg='gdot0_twin ('//PLASTICITY_PHENOPOWERLAW_label//')')
if ( constitutive_phenopowerlaw_n_twin(instance) <= 0.0_pReal .and. &
any(constitutive_phenopowerlaw_Ntwin(:,instance) > 0)) &
call IO_error(211_pInt,el=instance,ext_msg='n_twin ('//PLASTICITY_PHENOPOWERLAW_label//')')
if (constitutive_phenopowerlaw_aTolResistance(instance) <= 0.0_pReal) &
constitutive_phenopowerlaw_aTolResistance(instance) = 1.0_pReal ! default absolute tolerance 1 Pa
if (constitutive_phenopowerlaw_aTolShear(instance) <= 0.0_pReal) &
constitutive_phenopowerlaw_aTolShear(instance) = 1.0e-6_pReal ! default absolute tolerance 1e-6
if (constitutive_phenopowerlaw_aTolTwinfrac(instance) <= 0.0_pReal) &
constitutive_phenopowerlaw_aTolTwinfrac(instance) = 1.0e-6_pReal ! default absolute tolerance 1e-6
endif myPhase
enddo sanityChecks
!--------------------------------------------------------------------------------------------------
@ -459,95 +444,95 @@ subroutine constitutive_phenopowerlaw_init(fileUnit)
maxval(constitutive_phenopowerlaw_totalNtwin),&
maxNinstance), source=0.0_pReal)
instancesLoop: do instance = 1_pInt,maxNinstance
outputsLoop: do o = 1_pInt,constitutive_phenopowerlaw_Noutput(instance)
select case(constitutive_phenopowerlaw_outputID(o,instance))
case(resistance_slip_ID, &
shearrate_slip_ID, &
accumulatedshear_slip_ID, &
resolvedstress_slip_ID &
)
mySize = constitutive_phenopowerlaw_totalNslip(instance)
case(resistance_twin_ID, &
shearrate_twin_ID, &
accumulatedshear_twin_ID, &
resolvedstress_twin_ID &
)
mySize = constitutive_phenopowerlaw_totalNtwin(instance)
case(totalshear_ID, &
totalvolfrac_ID &
)
mySize = 1_pInt
case default
end select
outputFound: if (mySize > 0_pInt) then
constitutive_phenopowerlaw_sizePostResult(o,instance) = mySize
constitutive_phenopowerlaw_sizePostResults(instance) = constitutive_phenopowerlaw_sizePostResults(instance) + mySize
endif outputFound
enddo outputsLoop
constitutive_phenopowerlaw_sizeDotState(instance) = constitutive_phenopowerlaw_totalNslip(instance)+ &
constitutive_phenopowerlaw_totalNtwin(instance)+ &
2_pInt + &
constitutive_phenopowerlaw_totalNslip(instance)+ &
constitutive_phenopowerlaw_totalNtwin(instance) ! s_slip, s_twin, sum(gamma), sum(f), accshear_slip, accshear_twin
constitutive_phenopowerlaw_sizeState(instance) = constitutive_phenopowerlaw_sizeDotState(instance)
structID = constitutive_phenopowerlaw_structure(instance)
do f = 1_pInt,lattice_maxNslipFamily ! >>> interaction slip -- X
index_myFamily = sum(constitutive_phenopowerlaw_Nslip(1:f-1_pInt,instance))
do j = 1_pInt,constitutive_phenopowerlaw_Nslip(f,instance) ! loop over (active) systems in my family (slip)
do o = 1_pInt,lattice_maxNslipFamily
index_otherFamily = sum(constitutive_phenopowerlaw_Nslip(1:o-1_pInt,instance))
do k = 1_pInt,constitutive_phenopowerlaw_Nslip(o,instance) ! loop over (active) systems in other family (slip)
constitutive_phenopowerlaw_hardeningMatrix_SlipSlip(index_myFamily+j,index_otherFamily+k,instance) = &
constitutive_phenopowerlaw_interaction_SlipSlip(lattice_interactionSlipSlip( &
sum(lattice_NslipSystem(1:f-1,structID))+j, &
sum(lattice_NslipSystem(1:o-1,structID))+k, &
structID), instance )
enddo; enddo
do o = 1_pInt,lattice_maxNtwinFamily
index_otherFamily = sum(constitutive_phenopowerlaw_Ntwin(1:o-1_pInt,instance))
do k = 1_pInt,constitutive_phenopowerlaw_Ntwin(o,instance) ! loop over (active) systems in other family (twin)
constitutive_phenopowerlaw_hardeningMatrix_SlipTwin(index_myFamily+j,index_otherFamily+k,instance) = &
constitutive_phenopowerlaw_interaction_SlipTwin(lattice_interactionSlipTwin( &
sum(lattice_NslipSystem(1:f-1_pInt,structID))+j, &
sum(lattice_NtwinSystem(1:o-1_pInt,structID))+k, &
structID), instance )
enddo; enddo
enddo; enddo
do f = 1_pInt,lattice_maxNtwinFamily ! >>> interaction twin -- X
index_myFamily = sum(constitutive_phenopowerlaw_Ntwin(1:f-1_pInt,instance))
do j = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,instance) ! loop over (active) systems in my family (twin)
do o = 1_pInt,lattice_maxNslipFamily
index_otherFamily = sum(constitutive_phenopowerlaw_Nslip(1:o-1_pInt,instance))
do k = 1_pInt,constitutive_phenopowerlaw_Nslip(o,instance) ! loop over (active) systems in other family (slip)
constitutive_phenopowerlaw_hardeningMatrix_TwinSlip(index_myFamily+j,index_otherFamily+k,instance) = &
constitutive_phenopowerlaw_interaction_TwinSlip(lattice_interactionTwinSlip( &
sum(lattice_NtwinSystem(1:f-1_pInt,structID))+j, &
sum(lattice_NslipSystem(1:o-1_pInt,structID))+k, &
structID), instance )
enddo; enddo
do o = 1_pInt,lattice_maxNtwinFamily
index_otherFamily = sum(constitutive_phenopowerlaw_Ntwin(1:o-1_pInt,instance))
do k = 1_pInt,constitutive_phenopowerlaw_Ntwin(o,instance) ! loop over (active) systems in other family (twin)
constitutive_phenopowerlaw_hardeningMatrix_TwinTwin(index_myFamily+j,index_otherFamily+k,instance) = &
constitutive_phenopowerlaw_interaction_TwinTwin(lattice_interactionTwinTwin( &
sum(lattice_NtwinSystem(1:f-1_pInt,structID))+j, &
sum(lattice_NtwinSystem(1:o-1_pInt,structID))+k, &
structID), instance )
enddo; enddo
enddo; enddo
enddo instancesLoop
initializeInstances: do phase = 1_pInt, size(phase_plasticity)
if (phase_plasticity(phase) == PLASTICITY_phenopowerlaw_ID) then
instance = phase_plasticityInstance(phase)
outputsLoop: do o = 1_pInt,constitutive_phenopowerlaw_Noutput(instance)
select case(constitutive_phenopowerlaw_outputID(o,instance))
case(resistance_slip_ID, &
shearrate_slip_ID, &
accumulatedshear_slip_ID, &
resolvedstress_slip_ID &
)
mySize = constitutive_phenopowerlaw_totalNslip(instance)
case(resistance_twin_ID, &
shearrate_twin_ID, &
accumulatedshear_twin_ID, &
resolvedstress_twin_ID &
)
mySize = constitutive_phenopowerlaw_totalNtwin(instance)
case(totalshear_ID, &
totalvolfrac_ID &
)
mySize = 1_pInt
case default
end select
outputFound: if (mySize > 0_pInt) then
constitutive_phenopowerlaw_sizePostResult(o,instance) = mySize
constitutive_phenopowerlaw_sizePostResults(instance) = constitutive_phenopowerlaw_sizePostResults(instance) + mySize
endif outputFound
enddo outputsLoop
constitutive_phenopowerlaw_sizeDotState(instance) = constitutive_phenopowerlaw_totalNslip(instance)+ &
constitutive_phenopowerlaw_totalNtwin(instance)+ &
2_pInt + &
constitutive_phenopowerlaw_totalNslip(instance)+ &
constitutive_phenopowerlaw_totalNtwin(instance) ! s_slip, s_twin, sum(gamma), sum(f), accshear_slip, accshear_twin
constitutive_phenopowerlaw_sizeState(instance) = constitutive_phenopowerlaw_sizeDotState(instance)
do f = 1_pInt,lattice_maxNslipFamily ! >>> interaction slip -- X
index_myFamily = sum(constitutive_phenopowerlaw_Nslip(1:f-1_pInt,instance))
do j = 1_pInt,constitutive_phenopowerlaw_Nslip(f,instance) ! loop over (active) systems in my family (slip)
do o = 1_pInt,lattice_maxNslipFamily
index_otherFamily = sum(constitutive_phenopowerlaw_Nslip(1:o-1_pInt,instance))
do k = 1_pInt,constitutive_phenopowerlaw_Nslip(o,instance) ! loop over (active) systems in other family (slip)
constitutive_phenopowerlaw_hardeningMatrix_SlipSlip(index_myFamily+j,index_otherFamily+k,instance) = &
constitutive_phenopowerlaw_interaction_SlipSlip(lattice_interactionSlipSlip( &
sum(lattice_NslipSystem(1:f-1,phase))+j, &
sum(lattice_NslipSystem(1:o-1,phase))+k, &
phase), instance )
enddo; enddo
do o = 1_pInt,lattice_maxNtwinFamily
index_otherFamily = sum(constitutive_phenopowerlaw_Ntwin(1:o-1_pInt,instance))
do k = 1_pInt,constitutive_phenopowerlaw_Ntwin(o,instance) ! loop over (active) systems in other family (twin)
constitutive_phenopowerlaw_hardeningMatrix_SlipTwin(index_myFamily+j,index_otherFamily+k,instance) = &
constitutive_phenopowerlaw_interaction_SlipTwin(lattice_interactionSlipTwin( &
sum(lattice_NslipSystem(1:f-1_pInt,phase))+j, &
sum(lattice_NtwinSystem(1:o-1_pInt,phase))+k, &
phase), instance )
enddo; enddo
enddo; enddo
do f = 1_pInt,lattice_maxNtwinFamily ! >>> interaction twin -- X
index_myFamily = sum(constitutive_phenopowerlaw_Ntwin(1:f-1_pInt,instance))
do j = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,instance) ! loop over (active) systems in my family (twin)
do o = 1_pInt,lattice_maxNslipFamily
index_otherFamily = sum(constitutive_phenopowerlaw_Nslip(1:o-1_pInt,instance))
do k = 1_pInt,constitutive_phenopowerlaw_Nslip(o,instance) ! loop over (active) systems in other family (slip)
constitutive_phenopowerlaw_hardeningMatrix_TwinSlip(index_myFamily+j,index_otherFamily+k,instance) = &
constitutive_phenopowerlaw_interaction_TwinSlip(lattice_interactionTwinSlip( &
sum(lattice_NtwinSystem(1:f-1_pInt,phase))+j, &
sum(lattice_NslipSystem(1:o-1_pInt,phase))+k, &
phase), instance )
enddo; enddo
do o = 1_pInt,lattice_maxNtwinFamily
index_otherFamily = sum(constitutive_phenopowerlaw_Ntwin(1:o-1_pInt,instance))
do k = 1_pInt,constitutive_phenopowerlaw_Ntwin(o,instance) ! loop over (active) systems in other family (twin)
constitutive_phenopowerlaw_hardeningMatrix_TwinTwin(index_myFamily+j,index_otherFamily+k,instance) = &
constitutive_phenopowerlaw_interaction_TwinTwin(lattice_interactionTwinTwin( &
sum(lattice_NtwinSystem(1:f-1_pInt,phase))+j, &
sum(lattice_NtwinSystem(1:o-1_pInt,phase))+k, &
phase), instance )
enddo; enddo
enddo; enddo
endif
enddo initializeInstances
end subroutine constitutive_phenopowerlaw_init
@ -562,7 +547,7 @@ pure function constitutive_phenopowerlaw_stateInit(instance)
implicit none
integer(pInt), intent(in) :: &
instance !< number specifying the instance of the plasticity
instance !< number specifying the instance of the plasticity
real(pReal), dimension(constitutive_phenopowerlaw_sizeDotState(instance)) :: &
constitutive_phenopowerlaw_stateInit
integer(pInt) :: &
@ -594,9 +579,9 @@ end function constitutive_phenopowerlaw_stateInit
pure function constitutive_phenopowerlaw_aTolState(instance)
implicit none
integer(pInt), intent(in) :: instance !< number specifying the instance of the plasticity
integer(pInt), intent(in) :: instance !< number specifying the instance of the plasticity
real(pReal), dimension(constitutive_phenopowerlaw_sizeState(instance)) :: &
real(pReal), dimension(constitutive_phenopowerlaw_sizeState(instance)) :: &
constitutive_phenopowerlaw_aTolState
constitutive_phenopowerlaw_aTolState(1:constitutive_phenopowerlaw_totalNslip(instance)+ &
@ -617,34 +602,6 @@ real(pReal), dimension(constitutive_phenopowerlaw_sizeState(instance)) :: &
end function constitutive_phenopowerlaw_aTolState
!--------------------------------------------------------------------------------------------------
!> @brief returns the homogenized elasticity matrix
!--------------------------------------------------------------------------------------------------
pure function constitutive_phenopowerlaw_homogenizedC(ipc,ip,el)
use prec, only: &
p_vec
use mesh, only: &
mesh_NcpElems, &
mesh_maxNips
use material, only: &
homogenization_maxNgrains, &
material_phase, &
phase_plasticityInstance
implicit none
real(pReal), dimension(6,6) :: &
constitutive_phenopowerlaw_homogenizedC
integer(pInt), intent(in) :: &
ipc, & !< component-ID of integration point
ip, & !< integration point
el !< element
constitutive_phenopowerlaw_homogenizedC = constitutive_phenopowerlaw_Cslip_66(1:6,1:6,&
phase_plasticityInstance(material_phase(ipc,ip,el)))
end function constitutive_phenopowerlaw_homogenizedC
!--------------------------------------------------------------------------------------------------
!> @brief calculates plastic velocity gradient and its tangent
!--------------------------------------------------------------------------------------------------
@ -690,7 +647,7 @@ pure subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar
integer(pInt) :: &
instance, &
nSlip, &
nTwin,structID,index_Gamma,index_F,index_myFamily, &
nTwin,phase,index_Gamma,index_F,index_myFamily, &
f,i,j,k,l,m,n
real(pReal), dimension(3,3,3,3) :: &
dLp_dTstar3333 !< derivative of Lp with respect to Tstar as 4th order tensor
@ -701,12 +658,11 @@ pure subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar
real(pReal), dimension(constitutive_phenopowerlaw_totalNtwin(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
gdot_twin,dgdot_dtautwin,tau_twin
instance = phase_plasticityInstance(material_phase(ipc,ip,el))
structID = constitutive_phenopowerlaw_structure(instance)
phase = material_phase(ipc,ip,el)
instance = phase_plasticityInstance(phase)
nSlip = constitutive_phenopowerlaw_totalNslip(instance)
nTwin = constitutive_phenopowerlaw_totalNtwin(instance)
index_Gamma = nSlip + nTwin + 1_pInt
index_F = nSlip + nTwin + 2_pInt
@ -716,25 +672,25 @@ pure subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar
j = 0_pInt
slipFamiliesLoop: do f = 1_pInt,lattice_maxNslipFamily
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,structID)) ! at which index starts my family
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,phase)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,instance) ! process each (active) slip system in family
j = j+1_pInt
!--------------------------------------------------------------------------------------------------
! Calculation of Lp
tau_slip_pos(j) = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,structID))
tau_slip_pos(j) = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,phase))
tau_slip_neg(j) = tau_slip_pos(j)
nonSchmid_tensor(1:3,1:3,1) = lattice_Sslip(1:3,1:3,1,index_myFamily+i,structID)
nonSchmid_tensor(1:3,1:3,1) = lattice_Sslip(1:3,1:3,1,index_myFamily+i,phase)
nonSchmid_tensor(1:3,1:3,2) = nonSchmid_tensor(1:3,1:3,1)
do k = 1,lattice_NnonSchmid(structID)
do k = 1,lattice_NnonSchmid(phase)
tau_slip_pos(j) = tau_slip_pos(j) + constitutive_phenopowerlaw_nonSchmidCoeff(k,instance)* &
dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k,index_myFamily+i,structID))
dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k,index_myFamily+i,phase))
tau_slip_neg(j) = tau_slip_neg(j) + constitutive_phenopowerlaw_nonSchmidCoeff(k,instance)* &
dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,structID))
dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,phase))
nonSchmid_tensor(1:3,1:3,1) = nonSchmid_tensor(1:3,1:3,1) + constitutive_phenopowerlaw_nonSchmidCoeff(k,instance)*&
lattice_Sslip(1:3,1:3,2*k,index_myFamily+i,structID)
lattice_Sslip(1:3,1:3,2*k,index_myFamily+i,phase)
nonSchmid_tensor(1:3,1:3,2) = nonSchmid_tensor(1:3,1:3,2) + constitutive_phenopowerlaw_nonSchmidCoeff(k,instance)*&
lattice_Sslip(1:3,1:3,2*k+1,index_myFamily+i,structID)
lattice_Sslip(1:3,1:3,2*k+1,index_myFamily+i,phase)
enddo
gdot_slip_pos(j) = 0.5_pReal*constitutive_phenopowerlaw_gdot0_slip(instance)* &
((abs(tau_slip_pos(j))/state(ipc,ip,el)%p(j))**constitutive_phenopowerlaw_n_slip(instance))*&
@ -743,7 +699,7 @@ pure subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar
((abs(tau_slip_neg(j))/state(ipc,ip,el)%p(j))**constitutive_phenopowerlaw_n_slip(instance))*&
sign(1.0_pReal,tau_slip_neg(j))
Lp = Lp + (1.0_pReal-state(ipc,ip,el)%p(index_F))*& ! 1-F
(gdot_slip_pos(j)+gdot_slip_neg(j))*lattice_Sslip(1:3,1:3,1,index_myFamily+i,structID)
(gdot_slip_pos(j)+gdot_slip_neg(j))*lattice_Sslip(1:3,1:3,1,index_myFamily+i,phase)
!--------------------------------------------------------------------------------------------------
! Calculation of the tangent of Lp
@ -751,7 +707,7 @@ pure subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar
dgdot_dtauslip_pos(j) = gdot_slip_pos(j)*constitutive_phenopowerlaw_n_slip(instance)/tau_slip_pos(j)
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
dgdot_dtauslip_pos(j)*lattice_Sslip(k,l,1,index_myFamily+i,structID)* &
dgdot_dtauslip_pos(j)*lattice_Sslip(k,l,1,index_myFamily+i,phase)* &
nonSchmid_tensor(m,n,1)
endif
@ -759,7 +715,7 @@ pure subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar
dgdot_dtauslip_neg(j) = gdot_slip_neg(j)*constitutive_phenopowerlaw_n_slip(instance)/tau_slip_neg(j)
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
dgdot_dtauslip_neg(j)*lattice_Sslip(k,l,1,index_myFamily+i,structID)* &
dgdot_dtauslip_neg(j)*lattice_Sslip(k,l,1,index_myFamily+i,phase)* &
nonSchmid_tensor(m,n,2)
endif
enddo
@ -767,18 +723,18 @@ pure subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar
j = 0_pInt
twinFamiliesLoop: do f = 1_pInt,lattice_maxNtwinFamily
index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,structID)) ! at which index starts my family
index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,phase)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,instance) ! process each (active) twin system in family
j = j+1_pInt
!--------------------------------------------------------------------------------------------------
! Calculation of Lp
tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,structID))
tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,phase))
gdot_twin(j) = (1.0_pReal-state(ipc,ip,el)%p(index_F))*& ! 1-F
constitutive_phenopowerlaw_gdot0_twin(instance)*&
(abs(tau_twin(j))/state(ipc,ip,el)%p(nSlip+j))**&
constitutive_phenopowerlaw_n_twin(instance)*max(0.0_pReal,sign(1.0_pReal,tau_twin(j)))
Lp = Lp + gdot_twin(j)*lattice_Stwin(1:3,1:3,index_myFamily+i,structID)
Lp = Lp + gdot_twin(j)*lattice_Stwin(1:3,1:3,index_myFamily+i,phase)
!--------------------------------------------------------------------------------------------------
! Calculation of the tangent of Lp
@ -786,8 +742,8 @@ pure subroutine constitutive_phenopowerlaw_LpAndItsTangent(Lp,dLp_dTstar99,Tstar
dgdot_dtautwin(j) = gdot_twin(j)*constitutive_phenopowerlaw_n_twin(instance)/tau_twin(j)
forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
dgdot_dtautwin(j)*lattice_Stwin(k,l,index_myFamily+i,structID)* &
lattice_Stwin(m,n,index_myFamily+i,structID)
dgdot_dtautwin(j)*lattice_Stwin(k,l,index_myFamily+i,phase)* &
lattice_Stwin(m,n,index_myFamily+i,phase)
endif
enddo
enddo twinFamiliesLoop
@ -834,7 +790,7 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,state,ipc,ip,el)
constitutive_phenopowerlaw_dotState
integer(pInt) :: &
instance,structID, &
instance,phase, &
nSlip,nTwin, &
f,i,j,k, &
index_Gamma,index_F,index_myFamily, &
@ -848,9 +804,8 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,state,ipc,ip,el)
real(pReal), dimension(constitutive_phenopowerlaw_totalNtwin(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
gdot_twin,tau_twin,left_TwinSlip,left_TwinTwin,right_SlipTwin,right_TwinTwin
instance = phase_plasticityInstance(material_phase(ipc,ip,el))
structID = constitutive_phenopowerlaw_structure(instance)
phase = material_phase(ipc,ip,el)
instance = phase_plasticityInstance(phase)
nSlip = constitutive_phenopowerlaw_totalNslip(instance)
nTwin = constitutive_phenopowerlaw_totalNtwin(instance)
@ -878,8 +833,8 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,state,ipc,ip,el)
ssat_offset = constitutive_phenopowerlaw_spr(instance)*sqrt(state(ipc,ip,el)%p(index_F))
j = 0_pInt
slipFamiliesLoop1: do f = 1_pInt,lattice_maxNslipFamily
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,structID)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,instance) ! process each (active) slip system in family
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,phase)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,instance) ! process each (active) slip system in family
j = j+1_pInt
left_SlipSlip(j) = 1.0_pReal ! no system-dependent left part
left_SlipTwin(j) = 1.0_pReal ! no system-dependent left part
@ -892,13 +847,13 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,state,ipc,ip,el)
!--------------------------------------------------------------------------------------------------
! Calculation of dot gamma
tau_slip_pos(j) = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,structID))
tau_slip_pos(j) = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,phase))
tau_slip_neg(j) = tau_slip_pos(j)
do k = 1,lattice_NnonSchmid(structID)
do k = 1,lattice_NnonSchmid(phase)
tau_slip_pos(j) = tau_slip_pos(j) + constitutive_phenopowerlaw_nonSchmidCoeff(k,instance)* &
dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k,index_myFamily+i,structID))
dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k,index_myFamily+i,phase))
tau_slip_neg(j) = tau_slip_neg(j) + constitutive_phenopowerlaw_nonSchmidCoeff(k,instance)* &
dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,structID))
dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,phase))
enddo
gdot_slip(j) = constitutive_phenopowerlaw_gdot0_slip(instance)*0.5_pReal* &
((abs(tau_slip_pos(j))/state(ipc,ip,el)%p(j))**constitutive_phenopowerlaw_n_slip(instance) &
@ -909,8 +864,8 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,state,ipc,ip,el)
j = 0_pInt
twinFamiliesLoop1: do f = 1_pInt,lattice_maxNtwinFamily
index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,structID)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,instance) ! process each (active) twin system in family
index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,phase)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,instance) ! process each (active) twin system in family
j = j+1_pInt
left_TwinSlip(j) = 1.0_pReal ! no system-dependent right part
left_TwinTwin(j) = 1.0_pReal ! no system-dependent right part
@ -919,7 +874,7 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,state,ipc,ip,el)
!--------------------------------------------------------------------------------------------------
! Calculation of dot vol frac
tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,structID))
tau_twin(j) = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,phase))
gdot_twin(j) = (1.0_pReal-state(ipc,ip,el)%p(index_F))*& ! 1-F
constitutive_phenopowerlaw_gdot0_twin(instance)*&
(abs(tau_twin(j))/state(ipc,ip,el)%p(nSlip+j))**&
@ -931,7 +886,7 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,state,ipc,ip,el)
! calculate the overall hardening based on above
j = 0_pInt
slipFamiliesLoop2: do f = 1_pInt,lattice_maxNslipFamily
do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,instance) ! process each (active) slip system in family
do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,instance) ! process each (active) slip system in family
j = j+1_pInt
constitutive_phenopowerlaw_dotState(j) = & ! evolution of slip resistance j
c_SlipSlip * left_SlipSlip(j) * &
@ -948,8 +903,8 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,state,ipc,ip,el)
j = 0_pInt
twinFamiliesLoop2: do f = 1_pInt,lattice_maxNtwinFamily
index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,structID)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,instance) ! process each (active) twin system in family
index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,phase)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,instance) ! process each (active) twin system in family
j = j+1_pInt
constitutive_phenopowerlaw_dotState(j+nSlip) = & ! evolution of twin resistance j
c_TwinSlip * left_TwinSlip(j) * &
@ -960,7 +915,7 @@ function constitutive_phenopowerlaw_dotState(Tstar_v,state,ipc,ip,el)
right_TwinTwin*gdot_twin) ! dot gamma_twin modulated by right-side twin factor
if (state(ipc,ip,el)%p(index_F) < 0.98_pReal) & ! ensure twin volume fractions stays below 1.0
constitutive_phenopowerlaw_dotState(index_F) = constitutive_phenopowerlaw_dotState(index_F) + &
gdot_twin(j)/lattice_shearTwin(index_myFamily+i,structID)
gdot_twin(j)/lattice_shearTwin(index_myFamily+i,phase)
constitutive_phenopowerlaw_dotState(offset_accshear_twin+j) = abs(gdot_twin(j))
enddo
enddo twinFamiliesLoop2
@ -1008,16 +963,15 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,state,ipc,ip,el)
constitutive_phenopowerlaw_postResults
integer(pInt) :: &
instance,structID, &
instance,phase, &
nSlip,nTwin, &
o,f,i,c,j,k, &
index_Gamma,index_F,index_accshear_slip,index_accshear_twin,index_myFamily
real(pReal) :: &
tau_slip_pos,tau_slip_neg,tau
instance = phase_plasticityInstance(material_phase(ipc,ip,el))
structID = constitutive_phenopowerlaw_structure(instance)
phase = material_phase(ipc,ip,el)
instance = phase_plasticityInstance(phase)
nSlip = constitutive_phenopowerlaw_totalNslip(instance)
nTwin = constitutive_phenopowerlaw_totalNtwin(instance)
@ -1044,16 +998,16 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,state,ipc,ip,el)
case (shearrate_slip_ID)
j = 0_pInt
slipFamiliesLoop1: do f = 1_pInt,lattice_maxNslipFamily
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,structID)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,instance) ! process each (active) slip system in family
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,phase)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,instance) ! process each (active) slip system in family
j = j + 1_pInt
tau_slip_pos = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,structID))
tau_slip_pos = dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,phase))
tau_slip_neg = tau_slip_pos
do k = 1,lattice_NnonSchmid(structID)
do k = 1,lattice_NnonSchmid(phase)
tau_slip_pos = tau_slip_pos + constitutive_phenopowerlaw_nonSchmidCoeff(k,instance)* &
dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k,index_myFamily+i,structID))
dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k,index_myFamily+i,phase))
tau_slip_neg = tau_slip_neg + constitutive_phenopowerlaw_nonSchmidCoeff(k,instance)* &
dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,structID))
dot_product(Tstar_v,lattice_Sslip_v(1:6,2*k+1,index_myFamily+i,phase))
enddo
constitutive_phenopowerlaw_postResults(c+j) = constitutive_phenopowerlaw_gdot0_slip(instance)*0.5_pReal* &
((abs(tau_slip_pos)/state(ipc,ip,el)%p(j))**constitutive_phenopowerlaw_n_slip(instance) &
@ -1066,11 +1020,11 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,state,ipc,ip,el)
case (resolvedstress_slip_ID)
j = 0_pInt
slipFamiliesLoop2: do f = 1_pInt,lattice_maxNslipFamily
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,structID)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,instance) ! process each (active) slip system in family
index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,phase)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Nslip(f,instance) ! process each (active) slip system in family
j = j + 1_pInt
constitutive_phenopowerlaw_postResults(c+j) = &
dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,structID))
dot_product(Tstar_v,lattice_Sslip_v(1:6,1,index_myFamily+i,phase))
enddo
enddo slipFamiliesLoop2
c = c + nSlip
@ -1093,11 +1047,11 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,state,ipc,ip,el)
case (shearrate_twin_ID)
j = 0_pInt
twinFamiliesLoop1: do f = 1_pInt,lattice_maxNtwinFamily
index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,structID)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,instance) ! process each (active) twin system in family
index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,phase)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,instance) ! process each (active) twin system in family
j = j + 1_pInt
tau = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,structID))
constitutive_phenopowerlaw_postResults(c+j) = (1.0_pReal-state(ipc,ip,el)%p(index_F))*& ! 1-F
tau = dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,phase))
constitutive_phenopowerlaw_postResults(c+j) = (1.0_pReal-state(ipc,ip,el)%p(index_F))*& ! 1-F
constitutive_phenopowerlaw_gdot0_twin(instance)*&
(abs(tau)/state(ipc,ip,el)%p(j+nSlip))**&
constitutive_phenopowerlaw_n_twin(instance)*max(0.0_pReal,sign(1.0_pReal,tau))
@ -1108,11 +1062,11 @@ pure function constitutive_phenopowerlaw_postResults(Tstar_v,state,ipc,ip,el)
case (resolvedstress_twin_ID)
j = 0_pInt
twinFamiliesLoop2: do f = 1_pInt,lattice_maxNtwinFamily
index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,structID)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,instance) ! process each (active) twin system in family
index_myFamily = sum(lattice_NtwinSystem(1:f-1_pInt,phase)) ! at which index starts my family
do i = 1_pInt,constitutive_phenopowerlaw_Ntwin(f,instance) ! process each (active) twin system in family
j = j + 1_pInt
constitutive_phenopowerlaw_postResults(c+j) = &
dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,structID))
dot_product(Tstar_v,lattice_Stwin_v(1:6,index_myFamily+i,phase))
enddo
enddo twinFamiliesLoop2
c = c + nTwin

File diff suppressed because it is too large Load Diff

View File

@ -40,8 +40,6 @@ module crystallite
crystallite_sizePostResults !< description not available
integer(pInt), dimension(:,:), allocatable, private :: &
crystallite_sizePostResult !< description not available
integer(pInt), dimension(:,:,:), allocatable, private :: &
crystallite_symmetryID !< crystallographic symmetry 1=cubic 2=hexagonal, needed in all orientation calcs
real(pReal), dimension(:,:), allocatable, public :: &
crystallite_temperature !< temperature (same on all components on one IP)
@ -189,16 +187,9 @@ subroutine crystallite_init(temperature)
IO_EOF
use material
use lattice, only: &
lattice_symmetryType, &
lattice_structureID
lattice_structure
use constitutive, only: &
constitutive_microstructure
use constitutive_dislotwin, only: &
constitutive_dislotwin_structureID
use constitutive_titanmod, only: &
constitutive_titanmod_structureID
use constitutive_nonlocal, only: &
constitutive_nonlocal_structureID
implicit none
real(pReal), intent(in) :: temperature
@ -220,9 +211,8 @@ subroutine crystallite_init(temperature)
j, &
p, &
output = 0_pInt, &
mySize, &
myPhase, &
myMat
mySize
character(len=65536) :: &
tag = '', &
line= ''
@ -272,7 +262,6 @@ subroutine crystallite_init(temperature)
allocate(crystallite_orientation0(4,gMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_rotation(4,gMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_disorientation(4,nMax,gMax,iMax,eMax), source=0.0_pReal)
allocate(crystallite_symmetryID(gMax,iMax,eMax), source=0_pInt)
allocate(crystallite_localPlasticity(gMax,iMax,eMax), source=.true.)
allocate(crystallite_requested(gMax,iMax,eMax), source=.false.)
allocate(crystallite_todo(gMax,iMax,eMax), source=.false.)
@ -434,33 +423,6 @@ subroutine crystallite_init(temperature)
crystallite_partionedF0 = crystallite_F0
crystallite_partionedF = crystallite_F0
!--------------------------------------------------------------------------------------------------
! Initialize crystallite_symmetryID
do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(mesh_element(3,e))
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e)
do g = 1_pInt,myNgrains
myPhase = material_phase(g,i,e)
myMat = phase_plasticityInstance(myPhase)
select case (phase_plasticity(myPhase))
case (PLASTICITY_PHENOPOWERLAW_ID)
crystallite_symmetryID(g,i,e) = lattice_symmetryType(lattice_structureID(myPhase))
case (PLASTICITY_TITANMOD_ID)
crystallite_symmetryID(g,i,e) = &
lattice_symmetryType(constitutive_titanmod_structureID(myMat))
case (PLASTICITY_DISLOTWIN_ID)
crystallite_symmetryID(g,i,e) = &
lattice_symmetryType(constitutive_dislotwin_structureID(myMat))
case (PLASTICITY_NONLOCAL_ID)
crystallite_symmetryID(g,i,e) = &
lattice_symmetryType(constitutive_nonlocal_structureID(myMat))
case default
crystallite_symmetryID(g,i,e) = 0_pInt !< @ToDo: does this happen for j2 material?
end select
enddo
enddo
enddo
call crystallite_orientations()
crystallite_orientation0 = crystallite_orientation ! store initial orientations for calculation of grain rotations
@ -496,7 +458,6 @@ subroutine crystallite_init(temperature)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedFp0: ', shape(crystallite_partionedFp0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_partionedLp0: ', shape(crystallite_partionedLp0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_subF: ', shape(crystallite_subF)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_symmetryID: ', shape(crystallite_symmetryID)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_subF0: ', shape(crystallite_subF0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_subFe0: ', shape(crystallite_subFe0)
write(6,'(a35,1x,7(i8,1x))') 'crystallite_subFp0: ', shape(crystallite_subFp0)
@ -3060,7 +3021,6 @@ logical function crystallite_integrateStress(&
logical error ! flag indicating an error
integer(pInt) NiterationStress, & ! number of stress integrations
ierr, & ! error indicator for LAPACK
n, &
o, &
p, &
jacoCounter ! counter to check for Jacobian update
@ -3342,7 +3302,8 @@ subroutine crystallite_orientations
FE_geomtype, &
FE_celltype
use lattice, only: &
lattice_qDisorientation
lattice_qDisorientation, &
lattice_structure
use constitutive_nonlocal, only: &
constitutive_nonlocal_structure, &
constitutive_nonlocal_updateCompatibility
@ -3387,8 +3348,7 @@ subroutine crystallite_orientations
orientation = math_RtoQ(transpose(R))
endif
crystallite_rotation(1:4,g,i,e) = lattice_qDisorientation(crystallite_orientation0(1:4,g,i,e), & ! active rotation from ori0
orientation, & ! to current orientation
0_pInt ) ! we don't want symmetry here
orientation) ! to current orientation (with no symmetry)
crystallite_orientation(1:4,g,i,e) = orientation
enddo
enddo
@ -3406,9 +3366,7 @@ subroutine crystallite_orientations
myPhase = material_phase(1,i,e) ! get my phase
if (.not. phase_localPlasticity(myPhase)) then ! if nonlocal model
myInstance = phase_plasticityInstance(myPhase)
myStructure = constitutive_nonlocal_structure(myInstance) ! get my crystal structure
! --- calculate disorientation between me and my neighbor ---
do n = 1_pInt,FE_NipNeighbors(FE_celltype(FE_geomtype(mesh_element(2,e)))) ! loop through my neighbors
@ -3423,7 +3381,7 @@ subroutine crystallite_orientations
crystallite_disorientation(:,n,1,i,e) = &
lattice_qDisorientation( crystallite_orientation(1:4,1,i,e), &
crystallite_orientation(1:4,1,neighboring_i,neighboring_e), &
crystallite_symmetryID(1,i,e)) ! calculate disorientation
lattice_structure(myPhase)) ! calculate disorientation for given symmetry
else ! for neighbor with different phase
crystallite_disorientation(:,n,1,i,e) = [0.0_pReal, 1.0_pReal, 0.0_pReal, 0.0_pReal] ! 180 degree rotation about 100 axis
endif

View File

@ -22,6 +22,7 @@
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief defines lattice structure definitions, slip and twin system definitions, Schimd matrix
!> calculation and non-Schmid behavior
!--------------------------------------------------------------------------------------------------
@ -75,24 +76,14 @@ module lattice
real(pReal), allocatable, dimension(:,:), protected, public :: &
lattice_shearTwin !< characteristic twin shear
integer(pInt), private :: &
lattice_Nhexagonal, & !< total # of hexagonal lattice structure (from tag CoverA_ratio)
lattice_Nstructure !< total # of lattice structures (1: fcc,2: bcc,3+: hexagonal)
integer(pInt), dimension(:,:), pointer, private :: &
interactionSlipSlip, &
interactionSlipTwin, &
interactionTwinSlip, &
interactionTwinTwin
integer(pInt), allocatable, dimension(:), protected, public :: &
lattice_NnonSchmid !< total # of non-Schmid contributions for each structure
!--------------------------------------------------------------------------------------------------
! fcc (1)
! fcc
integer(pInt), dimension(lattice_maxNslipFamily), parameter, public :: &
lattice_fcc_NslipSystem = int([12, 0, 0, 0, 0, 0],pInt) !< total # of slip systems per family for fcc
LATTICE_fcc_NslipSystem = int([12, 0, 0, 0, 0, 0],pInt) !< total # of slip systems per family for fcc
integer(pInt), dimension(lattice_maxNtwinFamily), parameter, public :: &
lattice_fcc_NtwinSystem = int([12, 0, 0, 0],pInt) !< total # of twin systems per family for fcc
@ -102,9 +93,6 @@ module lattice
lattice_fcc_Ntwin = 12_pInt, & ! sum(lattice_fcc_NtwinSystem) !< total # of twin systems for fcc
lattice_fcc_NnonSchmid = 0_pInt !< total # of non-Schmid contributions for fcc
integer(pInt), private :: &
lattice_fcc_Nstructure = 0_pInt
real(pReal), dimension(3+3,lattice_fcc_Nslip), parameter, private :: &
lattice_fcc_systemSlip = reshape(real([&
! Slip direction Plane normal
@ -142,7 +130,7 @@ module lattice
lattice_fcc_shearTwin = 0.5_pReal*sqrt(2.0_pReal) !< Twin system <112>{111} ??? Sorted according to Eisenlohr & Hantcherli
integer(pInt), dimension(2_pInt,lattice_fcc_Ntwin), parameter, public :: &
lattice_fcc_corellationTwinSlip = reshape(int( [&
lattice_fcc_twinNucleationSlipPair = reshape(int( [&
2,3, &
1,3, &
1,2, &
@ -157,7 +145,7 @@ module lattice
10,11 &
],pInt),[2_pInt,lattice_fcc_Ntwin])
integer(pInt), dimension(lattice_fcc_Nslip,lattice_fcc_Nslip), target, public :: &
integer(pInt), dimension(lattice_fcc_Nslip,lattice_fcc_Nslip), parameter, public :: &
lattice_fcc_interactionSlipSlip = reshape(int( [&
1,2,2,4,6,5,3,5,5,4,5,6, & ! ---> slip
2,1,2,6,4,5,5,4,6,5,3,5, & ! |
@ -178,7 +166,7 @@ module lattice
!< 4: Hirth locks
!< 5: glissile junctions
!< 6: Lomer locks
integer(pInt), dimension(lattice_fcc_Nslip,lattice_fcc_Ntwin), target, public :: &
integer(pInt), dimension(lattice_fcc_Nslip,lattice_fcc_Ntwin), parameter, public :: &
lattice_fcc_interactionSlipTwin = reshape(int( [&
1,1,1,3,3,3,2,2,2,3,3,3, & ! ---> twin
1,1,1,3,3,3,3,3,3,2,2,2, & ! |
@ -196,10 +184,10 @@ module lattice
!< 1: coplanar interaction
!< 2: screw trace between slip system and twin habit plane (easy cross slip)
!< 3: other interaction
integer(pInt), dimension(lattice_fcc_Ntwin,lattice_fcc_Nslip), target, public :: &
integer(pInt), dimension(lattice_fcc_Ntwin,lattice_fcc_Nslip), parameter, public :: &
lattice_fcc_interactionTwinSlip = 0_pInt !< Twin--Slip interaction types for fcc
integer(pInt), dimension(lattice_fcc_Ntwin,lattice_fcc_Ntwin), target, public :: &
integer(pInt), dimension(lattice_fcc_Ntwin,lattice_fcc_Ntwin), parameter,public :: &
lattice_fcc_interactionTwinTwin = reshape(int( [&
1,1,1,2,2,2,2,2,2,2,2,2, & ! ---> twin
1,1,1,2,2,2,2,2,2,2,2,2, & ! |
@ -218,7 +206,7 @@ module lattice
!--------------------------------------------------------------------------------------------------
! bcc (2)
! bcc
integer(pInt), dimension(lattice_maxNslipFamily), parameter, public :: &
lattice_bcc_NslipSystem = int([ 12, 12, 0, 0, 0, 0], pInt) !< total # of slip systems per family for bcc
@ -230,9 +218,6 @@ module lattice
lattice_bcc_Ntwin = 12_pInt, & ! sum(lattice_bcc_NtwinSystem) !< total # of twin systems for bcc
lattice_bcc_NnonSchmid = 6_pInt !< # of non-Schmid contributions for bcc. 6 known non schmid contributions for BCC (A. Koester, A. Ma, A. Hartmaier 2012)
integer(pInt), private :: &
lattice_bcc_Nstructure = 0_pInt
real(pReal), dimension(3+3,lattice_bcc_Nslip), parameter, private :: &
lattice_bcc_systemSlip = reshape(real([&
! Slip direction Plane normal
@ -309,7 +294,7 @@ module lattice
real(pReal), dimension(lattice_bcc_Ntwin), parameter, private :: &
lattice_bcc_shearTwin = 0.5_pReal*sqrt(2.0_pReal)
integer(pInt), dimension(lattice_bcc_Nslip,lattice_bcc_Nslip), target, public :: &
integer(pInt), dimension(lattice_bcc_Nslip,lattice_bcc_Nslip), parameter, public :: &
lattice_bcc_interactionSlipSlip = reshape(int( [&
1,2,6,6,5,4,4,3,4,3,5,4, 6,6,4,3,3,4,6,6,4,3,6,6, & ! ---> slip
2,1,6,6,4,3,5,4,5,4,4,3, 6,6,3,4,4,3,6,6,3,4,6,6, & ! |
@ -343,7 +328,7 @@ module lattice
!< 4: mixed-asymmetrical junction
!< 5: mixed-symmetrical junction
!< 6: edge junction
integer(pInt), dimension(lattice_bcc_Nslip,lattice_bcc_Ntwin), target, public :: &
integer(pInt), dimension(lattice_bcc_Nslip,lattice_bcc_Ntwin), parameter, public :: &
lattice_bcc_interactionSlipTwin = reshape(int( [&
3,3,3,2,2,3,3,3,3,2,3,3, & ! ---> twin
3,3,2,3,3,2,3,3,2,3,3,3, & ! |
@ -374,10 +359,10 @@ module lattice
!< 1: coplanar interaction
!< 2: screw trace between slip system and twin habit plane (easy cross slip)
!< 3: other interaction
integer(pInt), dimension(lattice_bcc_Ntwin,lattice_bcc_Nslip), target, public :: &
integer(pInt), dimension(lattice_bcc_Ntwin,lattice_bcc_Nslip), parameter, public :: &
lattice_bcc_interactionTwinSlip = 0_pInt !< Twin--slip interaction types for bcc @todo not implemented yet
integer(pInt), dimension(lattice_bcc_Ntwin,lattice_bcc_Ntwin), target, public :: &
integer(pInt), dimension(lattice_bcc_Ntwin,lattice_bcc_Ntwin), parameter, public :: &
lattice_bcc_interactionTwinTwin = reshape(int( [&
1,3,3,3,3,3,3,2,3,3,2,3, & ! ---> twin
3,1,3,3,3,3,2,3,3,3,3,2, & ! |
@ -398,7 +383,7 @@ module lattice
!--------------------------------------------------------------------------------------------------
! hex (3+)
! hex
integer(pInt), dimension(lattice_maxNslipFamily), parameter, public :: &
lattice_hex_NslipSystem = int([ 3, 3, 3, 6, 12, 6],pInt) !< # of slip systems per family for hex
@ -409,9 +394,6 @@ module lattice
lattice_hex_Nslip = 33_pInt, & ! sum(lattice_hex_NslipSystem), !< total # of slip systems for hex
lattice_hex_Ntwin = 24_pInt, & ! sum(lattice_hex_NtwinSystem) !< total # of twin systems for hex
lattice_hex_NnonSchmid = 0_pInt !< # of non-Schmid contributions for hex
integer(pInt), private :: &
lattice_hex_Nstructure = 0_pInt
real(pReal), dimension(4+4,lattice_hex_Nslip), parameter, private :: &
lattice_hex_systemSlip = reshape(real([&
@ -517,7 +499,7 @@ module lattice
4 &
],pInt),[lattice_hex_Ntwin])
integer(pInt), dimension(lattice_hex_Nslip,lattice_hex_Nslip), target, public :: &
integer(pInt), dimension(lattice_hex_Nslip,lattice_hex_Nslip), parameter, public :: &
lattice_hex_interactionSlipSlip = reshape(int( [&
1, 2, 2, 3, 3, 3, 7, 7, 7, 13,13,13,13,13,13, 21,21,21,21,21,21,21,21,21,21,21,21, 31,31,31,31,31,31, & ! ---> slip
2, 1, 2, 3, 3, 3, 7, 7, 7, 13,13,13,13,13,13, 21,21,21,21,21,21,21,21,21,21,21,21, 31,31,31,31,31,31, & ! |
@ -560,7 +542,7 @@ module lattice
!
],pInt),[lattice_hex_Nslip,lattice_hex_Nslip],order=[2,1]) !< Slip--slip interaction types for hex (32? in total)
integer(pInt), dimension(lattice_hex_Nslip,lattice_hex_Ntwin), target, public :: &
integer(pInt), dimension(lattice_hex_Nslip,lattice_hex_Ntwin), parameter, public :: &
lattice_hex_interactionSlipTwin = reshape(int( [&
1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, & ! --> twin
1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 4, 4, 4, 4, 4, 4, & ! |
@ -603,7 +585,7 @@ module lattice
!
],pInt),[lattice_hex_Nslip,lattice_hex_Ntwin],order=[2,1]) !< Slip--twin interaction types for hex (isotropic, 24 in total)
integer(pInt), dimension(lattice_hex_Ntwin,lattice_hex_Nslip), target, public :: &
integer(pInt), dimension(lattice_hex_Ntwin,lattice_hex_Nslip), parameter, public :: &
lattice_hex_interactionTwinSlip = reshape(int( [&
1, 1, 1, 5, 5, 5, 9, 9, 9, 13,13,13,13,13,13, 17,17,17,17,17,17,17,17,17,17,17,17, 21,21,21,21,21,21, & ! --> slip
1, 1, 1, 5, 5, 5, 9, 9, 9, 13,13,13,13,13,13, 17,17,17,17,17,17,17,17,17,17,17,17, 21,21,21,21,21,21, & ! |
@ -634,7 +616,7 @@ module lattice
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),[lattice_hex_Ntwin,lattice_hex_Nslip],order=[2,1]) !< Twin--twin interaction types for hex (isotropic, 20 in total)
integer(pInt), dimension(lattice_hex_Ntwin,lattice_hex_Ntwin), target, public :: &
integer(pInt), dimension(lattice_hex_Ntwin,lattice_hex_Ntwin), parameter, public :: &
lattice_hex_interactionTwinTwin = reshape(int( [&
1, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 7, 7, 7, 7, 7, 7, 13,13,13,13,13,13, & ! ---> twin
2, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 7, 7, 7, 7, 7, 7, 13,13,13,13,13,13, & ! |
@ -664,6 +646,13 @@ module lattice
20,20,20,20,20,20, 19,19,19,19,19,19, 18,18,18,18,18,18, 17,17,17,17,16,17, &
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),[lattice_hex_Ntwin,lattice_hex_Ntwin],order=[2,1]) !< Twin--slip interaction types for hex (isotropic, 16 in total)
real(pReal), dimension(:,:,:), allocatable, public, protected :: &
lattice_C66
real(pReal), dimension(:,:,:,:,:), allocatable, public, protected :: &
lattice_C3333
real(pReal), dimension(:), allocatable, public, protected :: &
lattice_mu, &
lattice_nu
enum, bind(c)
enumerator :: LATTICE_undefined_ID, &
LATTICE_iso_ID, &
@ -672,12 +661,8 @@ module lattice
LATTICE_hex_ID, &
LATTICE_ort_ID
end enum
integer(pInt), dimension(:), allocatable, public, protected :: &
lattice_structure
integer(kind(LATTICE_undefined_ID)), dimension(:), allocatable, public, protected :: &
lattice_structureID
real(pReal), dimension(:,:,:), allocatable, public, protected :: &
lattice_Cslip_66
lattice_structure
integer(pInt), dimension(2), parameter, private :: &
@ -789,26 +774,12 @@ real(pReal), dimension(4,36), parameter, private :: &
! [ 1.0,0.0,0.0,0.0 ],
! ]
character(len=*), parameter, public :: &
LATTICE_iso_label = 'iso', &
LATTICE_fcc_label = 'fcc', &
LATTICE_bcc_label = 'bcc', &
LATTICE_hex_label = 'hex', &
LATTICE_ort_label = 'ort'
public :: &
lattice_init, &
lattice_initializeStructure, &
lattice_symmetryType, &
lattice_symmetrizeC66, &
lattice_configNchunks, &
lattice_qDisorientation, &
LATTICE_undefined_ID, &
LATTICE_iso_ID, &
LATTICE_fcc_ID, &
LATTICE_bcc_ID, &
LATTICE_hex_ID, &
LATTICE_ort_ID
LATTICE_hex_ID
contains
@ -844,23 +815,17 @@ subroutine lattice_init
debug_level, &
debug_lattice, &
debug_levelBasic
use math, only: &
math_Mandel3333to66, &
math_Voigt66to3333
implicit none
integer(pInt), parameter :: FILEUNIT = 200_pInt
integer(pInt) :: Nsections
character(len=32) :: &
structure = ''
integer(pInt) :: Nphases
character(len=65536) :: &
tag = '', &
line = ''
integer(pInt), parameter :: MAXNCHUNKS = 2_pInt
integer(pInt), dimension(1+2*MAXNCHUNKS) :: positions
integer(pInt) :: section = 0_pInt,i
real(pReal), dimension(:), allocatable :: CoverA
real(pReal), dimension(:), allocatable :: CoverA !< c/a ratio for hex type lattice
write(6,'(/,a)') ' <<<+- lattice init -+>>>'
write(6,'(a)') ' $Id$'
@ -897,15 +862,40 @@ subroutine lattice_init
! read from material configuration file
if (.not. IO_open_jobFile_stat(FILEUNIT,material_localFileExt)) & ! no local material configuration present...
call IO_open_file(FILEUNIT,material_configFile) ! ... open material.config file
Nsections = IO_countSections(FILEUNIT,material_partPhase)
lattice_Nstructure = 2_pInt + sum(IO_countTagInPart(FILEUNIT,material_partPhase,'covera_ratio',Nsections)) ! fcc + bcc + all hex
Nphases = IO_countSections(FILEUNIT,material_partPhase)
allocate(lattice_structure(Nsections), source=0_pInt)
allocate(lattice_structureID(Nsections), source=LATTICE_undefined_ID)
allocate(lattice_Cslip_66(6,6,Nsections),source=0.0_pReal)
allocate(CoverA(Nsections), source=0.0_pReal)
allocate(lattice_structure(Nphases),source = LATTICE_undefined_ID)
allocate(lattice_C66(6,6,Nphases), source=0.0_pReal)
allocate(lattice_C3333(3,3,3,3,Nphases), source=0.0_pReal)
allocate(lattice_mu(Nphases), source=0.0_pReal)
allocate(lattice_nu(Nphases), source=0.0_pReal)
allocate(lattice_NnonSchmid(Nphases), source=0_pInt)
allocate(lattice_Sslip(3,3,1+2*lattice_maxNnonSchmid,lattice_maxNslip,Nphases),source=0.0_pReal)
allocate(lattice_Sslip_v(6,1+2*lattice_maxNnonSchmid,lattice_maxNslip,Nphases),source=0.0_pReal)
allocate(lattice_sd(3,lattice_maxNslip,Nphases),source=0.0_pReal)
allocate(lattice_st(3,lattice_maxNslip,Nphases),source=0.0_pReal)
allocate(lattice_sn(3,lattice_maxNslip,Nphases),source=0.0_pReal)
allocate(lattice_Qtwin(3,3,lattice_maxNtwin,Nphases),source=0.0_pReal)
allocate(lattice_Stwin(3,3,lattice_maxNtwin,Nphases),source=0.0_pReal)
allocate(lattice_Stwin_v(6,lattice_maxNtwin,Nphases),source=0.0_pReal)
allocate(lattice_td(3,lattice_maxNtwin,Nphases),source=0.0_pReal)
allocate(lattice_tt(3,lattice_maxNtwin,Nphases),source=0.0_pReal)
allocate(lattice_tn(3,lattice_maxNtwin,Nphases),source=0.0_pReal)
allocate(lattice_shearTwin(lattice_maxNtwin,Nphases),source=0.0_pReal)
allocate(lattice_NslipSystem(lattice_maxNslipFamily,Nphases),source=0_pInt)
allocate(lattice_NtwinSystem(lattice_maxNtwinFamily,Nphases),source=0_pInt)
allocate(lattice_interactionSlipSlip(lattice_maxNslip,lattice_maxNslip,Nphases),source=0_pInt)! other:me
allocate(lattice_interactionSlipTwin(lattice_maxNslip,lattice_maxNtwin,Nphases),source=0_pInt)! other:me
allocate(lattice_interactionTwinSlip(lattice_maxNtwin,lattice_maxNslip,Nphases),source=0_pInt)! other:me
allocate(lattice_interactionTwinTwin(lattice_maxNtwin,lattice_maxNtwin,Nphases),source=0_pInt)! other:me
allocate(CoverA(Nphases),source=0.0_pReal)
rewind(fileUnit)
line = '' ! to have it initialized
section = 0_pInt ! - " -
@ -928,80 +918,52 @@ subroutine lattice_init
tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key
select case(tag)
case ('lattice_structure')
structure = IO_lc(IO_stringValue(line,positions,2_pInt))
select case(structure(1:3))
case(LATTICE_iso_label)
lattice_structureID(section) = LATTICE_iso_ID
case(LATTICE_fcc_label)
lattice_structureID(section) = LATTICE_fcc_ID
case(LATTICE_bcc_label)
lattice_structureID(section) = LATTICE_bcc_ID
case(LATTICE_hex_label)
lattice_structureID(section) = LATTICE_hex_ID
case(LATTICE_ort_label)
lattice_structureID(section) = LATTICE_ort_ID
select case(trim(IO_lc(IO_stringValue(line,positions,2_pInt))))
case('iso','isotropic')
lattice_structure(section) = LATTICE_iso_ID
case('fcc')
lattice_structure(section) = LATTICE_fcc_ID
case('bcc')
lattice_structure(section) = LATTICE_bcc_ID
case('hex','hexagonal')
lattice_structure(section) = LATTICE_hex_ID
case('ort','orthorombic')
lattice_structure(section) = LATTICE_ort_ID
case default
!there should be an error here
!there will be an error here
end select
case ('c11')
lattice_Cslip_66(1,1,section) = IO_floatValue(line,positions,2_pInt)
lattice_C66(1,1,section) = IO_floatValue(line,positions,2_pInt)
case ('c12')
lattice_Cslip_66(1,2,section) = IO_floatValue(line,positions,2_pInt)
lattice_C66(1,2,section) = IO_floatValue(line,positions,2_pInt)
case ('c13')
lattice_Cslip_66(1,3,section) = IO_floatValue(line,positions,2_pInt)
lattice_C66(1,3,section) = IO_floatValue(line,positions,2_pInt)
case ('c22')
lattice_Cslip_66(2,2,section) = IO_floatValue(line,positions,2_pInt)
lattice_C66(2,2,section) = IO_floatValue(line,positions,2_pInt)
case ('c23')
lattice_Cslip_66(2,3,section) = IO_floatValue(line,positions,2_pInt)
lattice_C66(2,3,section) = IO_floatValue(line,positions,2_pInt)
case ('c33')
lattice_Cslip_66(3,3,section) = IO_floatValue(line,positions,2_pInt)
lattice_C66(3,3,section) = IO_floatValue(line,positions,2_pInt)
case ('c44')
lattice_Cslip_66(4,4,section) = IO_floatValue(line,positions,2_pInt)
lattice_C66(4,4,section) = IO_floatValue(line,positions,2_pInt)
case ('c55')
lattice_Cslip_66(5,5,section) = IO_floatValue(line,positions,2_pInt)
lattice_C66(5,5,section) = IO_floatValue(line,positions,2_pInt)
case ('c66')
lattice_Cslip_66(6,6,section) = IO_floatValue(line,positions,2_pInt)
case ('covera_ratio')
lattice_C66(6,6,section) = IO_floatValue(line,positions,2_pInt)
case ('covera_ratio','c/a_ratio','c/a')
CoverA(section) = IO_floatValue(line,positions,2_pInt)
if (CoverA(section) < 1.0_pReal .or. CoverA(section) > 2.0_pReal) call IO_error(206_pInt) ! checking physical significance of c/a
end select
endif
enddo
if (iand(debug_level(debug_lattice),debug_levelBasic) /= 0_pInt) then
write(6,'(a16,1x,i5)') ' # phases:',Nsections
write(6,'(a16,1x,i5,/)') ' # structures:',lattice_Nstructure
write(6,'(a16,1x,i5)') ' # phases:',Nphases
endif
allocate(lattice_NnonSchmid(lattice_Nstructure), source=0_pInt)
allocate(lattice_Sslip(3,3,1+2*lattice_maxNnonSchmid,lattice_maxNslip,lattice_Nstructure),source= 0.0_pReal)
allocate(lattice_Sslip_v(6,1+2*lattice_maxNnonSchmid,lattice_maxNslip,lattice_Nstructure),source = 0.0_pReal)
allocate(lattice_sd(3,lattice_maxNslip,lattice_Nstructure),source=0.0_pReal)
allocate(lattice_st(3,lattice_maxNslip,lattice_Nstructure),source=0.0_pReal)
allocate(lattice_sn(3,lattice_maxNslip,lattice_Nstructure),source=0.0_pReal)
allocate(lattice_Qtwin(3,3,lattice_maxNtwin,lattice_Nstructure),source= 0.0_pReal)
allocate(lattice_Stwin(3,3,lattice_maxNtwin,lattice_Nstructure),source= 0.0_pReal)
allocate(lattice_Stwin_v(6,lattice_maxNtwin,lattice_Nstructure),source= 0.0_pReal)
allocate(lattice_td(3,lattice_maxNtwin,lattice_Nstructure),source= 0.0_pReal)
allocate(lattice_tt(3,lattice_maxNtwin,lattice_Nstructure),source= 0.0_pReal)
allocate(lattice_tn(3,lattice_maxNtwin,lattice_Nstructure),source= 0.0_pReal)
allocate(lattice_shearTwin(lattice_maxNtwin,lattice_Nstructure),source= 0.0_pReal)
allocate(lattice_NslipSystem(lattice_maxNslipFamily,lattice_Nstructure), source=0_pInt)
allocate(lattice_NtwinSystem(lattice_maxNtwinFamily,lattice_Nstructure), source=0_pInt)
allocate(lattice_interactionSlipSlip(lattice_maxNslip,lattice_maxNslip,lattice_Nstructure), source=0_pInt)! other:me
allocate(lattice_interactionSlipTwin(lattice_maxNslip,lattice_maxNtwin,lattice_Nstructure), source=0_pInt)! other:me
allocate(lattice_interactionTwinSlip(lattice_maxNtwin,lattice_maxNslip,lattice_Nstructure), source=0_pInt)! other:me
allocate(lattice_interactionTwinTwin(lattice_maxNtwin,lattice_maxNtwin,lattice_Nstructure), source=0_pInt)! other:me
do i = 1_pInt,Nsections
lattice_structure(i) = lattice_initializeStructure(lattice_structureID(i), CoverA(i)) ! get structure
lattice_Cslip_66(1:6,1:6,i) = lattice_symmetrizeC66(lattice_structureID(i),lattice_Cslip_66(1:6,1:6,i))
lattice_Cslip_66(1:6,1:6,i) = math_Mandel3333to66(math_Voigt66to3333(lattice_Cslip_66(1:6,1:6,i))) ! Literature data is Voigt, DAMASK uses Mandel
enddo
do i = 1_pInt,Nphases
call lattice_initializeStructure(i, CoverA(i))
enddo
deallocate(CoverA)
@ -1011,7 +973,9 @@ end subroutine lattice_init
!--------------------------------------------------------------------------------------------------
!> @brief Calculation of Schmid matrices, etc.
!--------------------------------------------------------------------------------------------------
integer(pInt) function lattice_initializeStructure(struct_ID,CoverA)
subroutine lattice_initializeStructure(myPhase,CoverA)
use prec, only: &
tol_math_check
use math, only: &
math_vectorproduct, &
math_tensorproduct, &
@ -1020,72 +984,74 @@ integer(pInt) function lattice_initializeStructure(struct_ID,CoverA)
math_trace33, &
math_symmetric33, &
math_Mandel33to6, &
math_Mandel3333to66, &
math_Voigt66to3333, &
math_axisAngleToR, &
INRAD
use IO, only: &
IO_error
implicit none
integer(kind(LATTICE_fcc_ID)), intent(in) :: struct_ID
integer(pInt), intent(in) :: myPhase
real(pReal), intent(in) :: CoverA
real(pReal), dimension(3) :: sdU = 0.0_pReal, &
snU = 0.0_pReal, &
np = 0.0_pReal, &
nn = 0.0_pReal
real(pReal), dimension(3,lattice_maxNslip) :: sd = 0.0_pReal, &
sn = 0.0_pReal
real(pReal), dimension(3,3,2,lattice_maxNnonSchmid,lattice_maxNslip) :: sns = 0.0_pReal
real(pReal), dimension(3,lattice_maxNtwin) :: td = 0.0_pReal, &
tn = 0.0_pReal
real(pReal), dimension(lattice_maxNtwin) :: ts = 0.0_pReal
integer(pInt), dimension(lattice_maxNslipFamily) :: myNslipSystem = 0_pInt
integer(pInt), dimension(lattice_maxNtwinFamily) :: myNtwinSystem = 0_pInt
integer(pInt) :: i,j,myNslip,myNtwin,myStructure = 0_pInt
logical :: processMe
processMe = .false.
real(pReal), dimension(3) :: &
sdU, snU, &
np, nn
real(pReal), dimension(3,lattice_maxNslip) :: &
sd, sn
real(pReal), dimension(3,3,2,lattice_maxNnonSchmid,lattice_maxNslip) :: &
sns
real(pReal), dimension(3,lattice_maxNtwin) :: &
td, tn
real(pReal), dimension(lattice_maxNtwin) :: &
ts
integer(pInt) :: &
i,j, &
myNslip, myNtwin
select case(struct_ID)
lattice_C66(1:6,1:6,myPhase) = lattice_symmetrizeC66(lattice_structure(myPhase),lattice_C66(1:6,1:6,myPhase))
lattice_mu(myPhase) = 0.2_pReal * (lattice_C66(1,1,myPhase) - lattice_C66(1,2,myPhase) + 3.0_pReal*lattice_C66(4,4,myPhase)) ! (C11iso-C12iso)/2 with C11iso=(3*C11+2*C12+4*C44)/5 and C12iso=(C11+4*C12-2*C44)/5
lattice_nu(myPhase) = (lattice_C66(1,1,myPhase) + 4.0_pReal*lattice_C66(1,2,myPhase) - 2.0_pReal*lattice_C66(4,4,myPhase)) &
/ (4.0_pReal*lattice_C66(1,1,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
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
select case(lattice_structure(myPhase))
!--------------------------------------------------------------------------------------------------
! fcc
case (LATTICE_fcc_ID)
myStructure = 1_pInt
myNslipSystem = lattice_fcc_NslipSystem ! size of slip system families
myNtwinSystem = lattice_fcc_NtwinSystem ! size of twin system families
myNslip = lattice_fcc_Nslip ! overall number of slip systems
myNtwin = lattice_fcc_Ntwin ! overall number of twin systems
lattice_fcc_Nstructure = lattice_fcc_Nstructure + 1_pInt ! count fcc instances
if (lattice_fcc_Nstructure == 1_pInt) then ! me is first fcc structure
processMe = .true.
lattice_NnonSchmid(myStructure) = lattice_fcc_NnonSchmid ! Currently no known non Schmid contributions for FCC (to be changed later)
do i = 1_pInt,myNslip ! assign slip system vectors
myNslip = lattice_fcc_Nslip
myNtwin = lattice_fcc_Ntwin
do i = 1_pInt,lattice_fcc_Nslip ! assign slip system vectors
sd(1:3,i) = lattice_fcc_systemSlip(1:3,i)
sn(1:3,i) = lattice_fcc_systemSlip(4:6,i)
do j = 1_pInt,lattice_fcc_NnonSchmid
sns(1:3,1:3,1,j,i) = 0.0_pReal
sns(1:3,1:3,2,j,i) = 0.0_pReal
enddo
enddo
do i = 1_pInt,myNtwin ! assign twin system vectors and shears
do i = 1_pInt,lattice_fcc_Ntwin ! assign twin system vectors and shears
td(1:3,i) = lattice_fcc_systemTwin(1:3,i)
tn(1:3,i) = lattice_fcc_systemTwin(4:6,i)
ts(i) = lattice_fcc_shearTwin(i)
enddo
interactionSlipSlip => lattice_fcc_interactionSlipSlip
interactionSlipTwin => lattice_fcc_interactionSlipTwin
interactionTwinSlip => lattice_fcc_interactionTwinSlip
interactionTwinTwin => lattice_fcc_interactionTwinTwin
endif
print*, shape(lattice_NslipSystem),shape(lattice_fcc_NslipSystem)
lattice_NslipSystem(1:lattice_maxNslipFamily,myPhase) = lattice_fcc_NslipSystem
lattice_NtwinSystem(1:lattice_maxNtwinFamily,myPhase) = lattice_fcc_NtwinSystem
lattice_NnonSchmid(myPhase) = lattice_fcc_NnonSchmid
lattice_interactionSlipSlip(1:lattice_fcc_Nslip,1:lattice_fcc_Nslip,myPhase) = &
lattice_fcc_interactionSlipSlip
lattice_interactionSlipTwin(1:lattice_fcc_Nslip,1:lattice_fcc_Ntwin,myPhase) = &
lattice_fcc_interactionSlipTwin
lattice_interactionTwinSlip(1:lattice_fcc_Ntwin,1:lattice_fcc_Nslip,myPhase) = &
lattice_fcc_interactionTwinSlip
lattice_interactionTwinTwin(1:lattice_fcc_Ntwin,1:lattice_fcc_Ntwin,myPhase) = &
lattice_fcc_interactionTwinTwin
!--------------------------------------------------------------------------------------------------
! bcc
case (LATTICE_bcc_ID)
myStructure = 2_pInt
myNslipSystem = lattice_bcc_NslipSystem ! size of slip system families
myNtwinSystem = lattice_bcc_NtwinSystem ! size of twin system families
myNslip = lattice_bcc_Nslip ! overall number of slip systems
myNtwin = lattice_bcc_Ntwin ! overall number of twin systems
lattice_bcc_Nstructure = lattice_bcc_Nstructure + 1_pInt ! count bcc instances
if (lattice_bcc_Nstructure == 1_pInt) then ! me is first bcc structure
processMe = .true.
lattice_NnonSchmid(myStructure) = lattice_bcc_NnonSchmid
do i = 1_pInt,myNslip ! assign slip system vectors
myNslip = lattice_bcc_Nslip
myNtwin = lattice_bcc_Ntwin
do i = 1_pInt,lattice_bcc_Nslip ! assign slip system vectors
sd(1:3,i) = lattice_bcc_systemSlip(1:3,i)
sn(1:3,i) = lattice_bcc_systemSlip(4:6,i)
sdU = sd(1:3,i) / math_norm3(sd(1:3,i))
@ -1105,43 +1071,37 @@ integer(pInt) function lattice_initializeStructure(struct_ID,CoverA)
sns(1:3,1:3,1,6,i) = math_tensorproduct(sdU, sdU)
sns(1:3,1:3,2,6,i) = math_tensorproduct(-sdU, -sdU)
enddo
do i = 1_pInt,myNtwin ! assign twin system vectors and shears
do i = 1_pInt,lattice_bcc_Ntwin ! assign twin system vectors and shears
td(1:3,i) = lattice_bcc_systemTwin(1:3,i)
tn(1:3,i) = lattice_bcc_systemTwin(4:6,i)
ts(i) = lattice_bcc_shearTwin(i)
enddo
interactionSlipSlip => lattice_bcc_interactionSlipSlip
interactionSlipTwin => lattice_bcc_interactionSlipTwin
interactionTwinSlip => lattice_bcc_interactionTwinSlip
interactionTwinTwin => lattice_bcc_interactionTwinTwin
endif
lattice_NslipSystem(1:lattice_maxNslipFamily,myPhase) = lattice_bcc_NslipSystem
lattice_NtwinSystem(1:lattice_maxNtwinFamily,myPhase) = lattice_bcc_NtwinSystem
lattice_NnonSchmid(myPhase) = lattice_bcc_NnonSchmid
lattice_interactionSlipSlip(1:lattice_bcc_Nslip,1:lattice_bcc_Nslip,myPhase) = &
lattice_bcc_interactionSlipSlip
lattice_interactionSlipTwin(1:lattice_bcc_Nslip,1:lattice_bcc_Ntwin,myPhase) = &
lattice_bcc_interactionSlipTwin
lattice_interactionTwinSlip(1:lattice_bcc_Ntwin,1:lattice_bcc_Nslip,myPhase) = &
lattice_bcc_interactionTwinSlip
lattice_interactionTwinTwin(1:lattice_bcc_Ntwin,1:lattice_bcc_Ntwin,myPhase) = &
lattice_bcc_interactionTwinTwin
!--------------------------------------------------------------------------------------------------
! hex (including conversion from miller-bravais (a1=a2=a3=c) to miller (a, b, c) indices)
case (LATTICE_hex_ID)
if (CoverA < 1.0_pReal .or. CoverA > 2.0_pReal) call IO_error(206_pInt) ! checking physical significance of c/a
lattice_hex_Nstructure = lattice_hex_Nstructure + 1_pInt ! count instances of hex structures
myStructure = 2_pInt + lattice_hex_Nstructure ! 3,4,5,.. for hex
myNslipSystem = lattice_hex_NslipSystem ! size of slip system families
myNtwinSystem = lattice_hex_NtwinSystem ! size of twin system families
myNslip = lattice_hex_Nslip ! overall number of slip systems
myNtwin = lattice_hex_Ntwin ! overall number of twin systems
processMe = .true.
lattice_NnonSchmid(myStructure) = lattice_hex_NnonSchmid ! Currently no known non Schmid contributions for hex (to be changed later)
! converting from 4 axes coordinate system (a1=a2=a3=c) to ortho-hexagonal system (a, b, c)
do i = 1_pInt,myNslip
myNslip = lattice_hex_Nslip
myNtwin = lattice_hex_Ntwin
do i = 1_pInt,lattice_hex_Nslip ! assign slip system vectors
sd(1,i) = lattice_hex_systemSlip(1,i)*1.5_pReal ! direction [uvtw]->[3u/2 (u+2v)*sqrt(3)/2 w*(c/a)]
sd(2,i) = (lattice_hex_systemSlip(1,i)+2.0_pReal*lattice_hex_systemSlip(2,i))*(0.5_pReal*sqrt(3.0_pReal))
sd(3,i) = lattice_hex_systemSlip(4,i)*CoverA
sn(1,i) = lattice_hex_systemSlip(5,i) ! plane (hkil)->(h (h+2k)/sqrt(3) l/(c/a))
sn(2,i) = (lattice_hex_systemSlip(5,i)+2.0_pReal*lattice_hex_systemSlip(6,i))/sqrt(3.0_pReal)
sn(3,i) = lattice_hex_systemSlip(8,i)/CoverA
do j = 1_pInt,lattice_hex_NnonSchmid
sns(1:3,1:3,1,j,i) = 0.0_pReal
sns(1:3,1:3,2,j,i) = 0.0_pReal
enddo
enddo
do i = 1_pInt,myNtwin
enddo
do i = 1_pInt,lattice_hex_Ntwin ! assign twin system vectors and shears
td(1,i) = lattice_hex_systemTwin(1,i)*1.5_pReal
td(2,i) = (lattice_hex_systemTwin(1,i)+2.0_pReal*lattice_hex_systemTwin(2,i))*(0.5_pReal*sqrt(3.0_pReal))
td(3,i) = lattice_hex_systemTwin(4,i)*CoverA
@ -1159,100 +1119,99 @@ integer(pInt) function lattice_initializeStructure(struct_ID,CoverA)
ts(i) = 2.0_pReal*(CoverA*CoverA-2.0_pReal)/3.0_pReal/CoverA
end select
enddo
lattice_NslipSystem(1:lattice_maxNslipFamily,myPhase) = lattice_hex_NslipSystem
lattice_NtwinSystem(1:lattice_maxNtwinFamily,myPhase) = lattice_hex_NtwinSystem
lattice_NnonSchmid(myPhase) = lattice_hex_NnonSchmid
lattice_interactionSlipSlip(1:lattice_hex_Nslip,1:lattice_hex_Nslip,myPhase) = &
lattice_hex_interactionSlipSlip
lattice_interactionSlipTwin(1:lattice_hex_Nslip,1:lattice_hex_Ntwin,myPhase) = &
lattice_hex_interactionSlipTwin
lattice_interactionTwinSlip(1:lattice_hex_Ntwin,1:lattice_hex_Nslip,myPhase) = &
lattice_hex_interactionTwinSlip
lattice_interactionTwinTwin(1:lattice_hex_Ntwin,1:lattice_hex_Ntwin,myPhase) = &
lattice_hex_interactionTwinTwin
interactionSlipSlip => lattice_hex_interactionSlipSlip
interactionSlipTwin => lattice_hex_interactionSlipTwin
interactionTwinSlip => lattice_hex_interactionTwinSlip
interactionTwinTwin => lattice_hex_interactionTwinTwin
!--------------------------------------------------------------------------------------------------
! orthorombic and isotropic (no crystal plasticity)
case (LATTICE_ort_ID, LATTICE_iso_ID)
myNslip = 0_pInt
myNtwin = 0_pInt
!--------------------------------------------------------------------------------------------------
! something went wrong
case default
processMe = .false.
myStructure = 0_pInt
print*, 'error'
end select
if (processMe) then
if (myStructure > lattice_Nstructure) &
call IO_error(666_pInt,myStructure,ext_msg = 'structure index out of bounds') ! check for memory leakage
do i = 1_pInt,myNslip ! store slip system vectors and Schmid matrix for my structure
lattice_sd(1:3,i,myStructure) = sd(1:3,i)/math_norm3(sd(1:3,i)) ! make unit vector
lattice_sn(1:3,i,myStructure) = sn(1:3,i)/math_norm3(sn(1:3,i)) ! make unit vector
lattice_st(1:3,i,myStructure) = math_vectorproduct(lattice_sd(1:3,i,myStructure), &
lattice_sn(1:3,i,myStructure))
lattice_Sslip(1:3,1:3,1,i,myStructure) = math_tensorproduct(lattice_sd(1:3,i,myStructure), &
lattice_sn(1:3,i,myStructure))
do j = 1_pInt,lattice_NnonSchmid(myStructure)
lattice_Sslip(1:3,1:3,2*j ,i,myStructure) = sns(1:3,1:3,1,j,i)
lattice_Sslip(1:3,1:3,2*j+1,i,myStructure) = sns(1:3,1:3,2,j,i)
do i = 1_pInt,myNslip ! store slip system vectors and Schmid matrix for my structure
lattice_sd(1:3,i,myPhase) = sd(1:3,i)/math_norm3(sd(1:3,i)) ! make unit vector
lattice_sn(1:3,i,myPhase) = sn(1:3,i)/math_norm3(sn(1:3,i)) ! make unit vector
lattice_st(1:3,i,myPhase) = math_vectorproduct(lattice_sd(1:3,i,myPhase), &
lattice_sn(1:3,i,myPhase))
lattice_Sslip(1:3,1:3,1,i,myPhase) = math_tensorproduct(lattice_sd(1:3,i,myPhase), &
lattice_sn(1:3,i,myPhase))
do j = 1_pInt,lattice_NnonSchmid(myPhase)
lattice_Sslip(1:3,1:3,2*j ,i,myPhase) = sns(1:3,1:3,1,j,i)
lattice_Sslip(1:3,1:3,2*j+1,i,myPhase) = sns(1:3,1:3,2,j,i)
enddo
do j = 1_pInt,1_pInt+2_pInt*lattice_NnonSchmid(myStructure)
lattice_Sslip_v(1:6,j,i,myStructure) = &
math_Mandel33to6(math_symmetric33(lattice_Sslip(1:3,1:3,j,i,myStructure)))
do j = 1_pInt,1_pInt+2_pInt*lattice_NnonSchmid(myPhase)
lattice_Sslip_v(1:6,j,i,myPhase) = &
math_Mandel33to6(math_symmetric33(lattice_Sslip(1:3,1:3,j,i,myPhase)))
enddo
if (abs(math_trace33(lattice_Sslip(1:3,1:3,1,i,myStructure))) > 1.0e-8_pReal) &
call IO_error(0_pInt,myStructure,i,0_pInt,ext_msg = 'dilatational slip Schmid matrix')
if (abs(math_trace33(lattice_Sslip(1:3,1:3,1,i,myPhase))) > tol_math_check) &
call IO_error(0_pInt,myPhase,i,0_pInt,ext_msg = 'dilatational slip Schmid matrix')
enddo
do i = 1_pInt,myNtwin ! store twin system vectors and Schmid plus rotation matrix for my structure
lattice_td(1:3,i,myStructure) = td(1:3,i)/math_norm3(td(1:3,i)) ! make unit vector
lattice_tn(1:3,i,myStructure) = tn(1:3,i)/math_norm3(tn(1:3,i)) ! make unit vector
lattice_tt(1:3,i,myStructure) = math_vectorproduct(lattice_td(1:3,i,myStructure), &
lattice_tn(1:3,i,myStructure))
lattice_Stwin(1:3,1:3,i,myStructure) = math_tensorproduct(lattice_td(1:3,i,myStructure), &
lattice_tn(1:3,i,myStructure))
lattice_Stwin_v(1:6,i,myStructure) = math_Mandel33to6(math_symmetric33(lattice_Stwin(1:3,1:3,i,myStructure)))
lattice_Qtwin(1:3,1:3,i,myStructure) = math_axisAngleToR(tn(1:3,i),180.0_pReal*INRAD)
lattice_shearTwin(i,myStructure) = ts(i)
if (abs(math_trace33(lattice_Stwin(1:3,1:3,i,myStructure))) > 1.0e-8_pReal) &
call IO_error(0_pInt,myStructure,i,0_pInt,ext_msg = 'dilatational twin Schmid matrix')
enddo
lattice_NslipSystem(1:lattice_maxNslipFamily,myStructure) = myNslipSystem ! number of slip systems in each family
lattice_NtwinSystem(1:lattice_maxNtwinFamily,myStructure) = myNtwinSystem ! number of twin systems in each family
lattice_interactionSlipSlip(1:myNslip,1:myNslip,myStructure) = interactionSlipSlip(1:myNslip,1:myNslip)
lattice_interactionSlipTwin(1:myNslip,1:myNtwin,myStructure) = interactionSlipTwin(1:myNslip,1:myNtwin)
lattice_interactionTwinSlip(1:myNtwin,1:myNslip,myStructure) = interactionTwinSlip(1:myNtwin,1:myNslip)
lattice_interactionTwinTwin(1:myNtwin,1:myNtwin,myStructure) = interactionTwinTwin(1:myNtwin,1:myNtwin)
endif
lattice_initializeStructure = myStructure ! report my structure index back
end function lattice_initializeStructure
!--------------------------------------------------------------------------------------------------
!> @brief Maps structure to symmetry type
!> @details fcc(1) and bcc(2) are cubic(1) hex(3+) is hexagonal(2)
!--------------------------------------------------------------------------------------------------
integer(pInt) pure function lattice_symmetryType(struct_ID)
implicit none
integer(kind(LATTICE_fcc_ID)), intent(in) :: struct_ID
select case(struct_ID)
case (LATTICE_fcc_ID,LATTICE_bcc_ID)
lattice_symmetryType = 1_pInt
case (LATTICE_hex_ID)
lattice_symmetryType = 2_pInt
case default
lattice_symmetryType = 0_pInt
end select
return
lattice_td(1:3,i,myPhase) = td(1:3,i)/math_norm3(td(1:3,i)) ! make unit vector
lattice_tn(1:3,i,myPhase) = tn(1:3,i)/math_norm3(tn(1:3,i)) ! make unit vector
lattice_tt(1:3,i,myPhase) = math_vectorproduct(lattice_td(1:3,i,myPhase), &
lattice_tn(1:3,i,myPhase))
lattice_Stwin(1:3,1:3,i,myPhase) = math_tensorproduct(lattice_td(1:3,i,myPhase), &
lattice_tn(1:3,i,myPhase))
lattice_Stwin_v(1:6,i,myPhase) = math_Mandel33to6(math_symmetric33(lattice_Stwin(1:3,1:3,i,myPhase)))
lattice_Qtwin(1:3,1:3,i,myPhase) = math_axisAngleToR(tn(1:3,i),180.0_pReal*INRAD)
lattice_shearTwin(i,myPhase) = ts(i)
if (abs(math_trace33(lattice_Stwin(1:3,1:3,i,myPhase))) > tol_math_check) &
call IO_error(301_pInt,myPhase,ext_msg = 'dilatational twin Schmid matrix')
enddo
end function lattice_symmetryType
print*, lattice_Stwin
print*, lattice_C66
print*, lattice_C3333
print*, lattice_mu
print*, lattice_nu
print*, lattice_sd
print*, lattice_sn
print*, lattice_st
print*, lattice_sslip
print*, lattice_td
print*, lattice_tn
print*, lattice_tt
print*, lattice_sd
print*, lattice_sn
print*, lattice_stwin
print*, lattice_qtwin
print*, lattice_qtwin
print*, lattice_sheartwin
end subroutine lattice_initializeStructure
!--------------------------------------------------------------------------------------------------
!> @brief Symmetrizes stiffness matrix according to lattice type
!--------------------------------------------------------------------------------------------------
pure function lattice_symmetrizeC66(struct_ID,C66)
pure function lattice_symmetrizeC66(struct,C66)
implicit none
integer(kind(LATTICE_fcc_ID)), intent(in) :: struct_ID
integer(kind(LATTICE_undefined_ID)), intent(in) :: struct
real(pReal), dimension(6,6), intent(in) :: C66
real(pReal), dimension(6,6) :: lattice_symmetrizeC66
integer(pInt) :: j,k
lattice_symmetrizeC66 = 0.0_pReal
select case(struct_ID)
select case(struct)
case (LATTICE_iso_ID)
forall(k=1_pInt:3_pInt)
forall(j=1_pInt:3_pInt) lattice_symmetrizeC66(k,j) = C66(1,2)
@ -1291,6 +1250,8 @@ pure function lattice_symmetrizeC66(struct_ID,C66)
lattice_symmetrizeC66(4,4) = C66(4,4)
lattice_symmetrizeC66(5,5) = C66(5,5)
lattice_symmetrizeC66(6,6) = C66(6,6)
case default
lattice_symmetrizeC66 = C66
end select
end function lattice_symmetrizeC66
@ -1299,25 +1260,25 @@ pure function lattice_symmetrizeC66(struct_ID,C66)
!--------------------------------------------------------------------------------------------------
!> @brief figures whether unit quat falls into stereographic standard triangle
!--------------------------------------------------------------------------------------------------
logical pure function lattice_qInSST(Q, symmetryType)
logical pure function lattice_qInSST(Q, struct)
use math, only: &
math_qToRodrig
implicit none
real(pReal), dimension(4), intent(in) :: Q ! orientation
integer(pInt), intent(in) :: symmetryType ! Type of crystal symmetry; 1:cubic, 2:hexagonal
integer(kind(LATTICE_undefined_ID)), intent(in) :: struct ! lattice structure
real(pReal), dimension(3) :: Rodrig ! Rodrigues vector of Q
Rodrig = math_qToRodrig(Q)
if (any(Rodrig/=Rodrig)) then
lattice_qInSST = .false.
else
select case (symmetryType)
case (1_pInt)
select case (struct)
case (LATTICE_bcc_ID,LATTICE_fcc_ID)
lattice_qInSST = Rodrig(1) > Rodrig(2) .and. &
Rodrig(2) > Rodrig(3) .and. &
Rodrig(3) > 0.0_pReal
case (2_pInt)
case (LATTICE_hex_ID)
lattice_qInSST = Rodrig(1) > sqrt(3.0_pReal)*Rodrig(2) .and. &
Rodrig(2) > 0.0_pReal .and. &
Rodrig(3) > 0.0_pReal
@ -1332,103 +1293,67 @@ end function lattice_qInSST
!--------------------------------------------------------------------------------------------------
!> @brief calculates the disorientation for 2 unit quaternions
!--------------------------------------------------------------------------------------------------
function lattice_qDisorientation(Q1, Q2, symmetryType)
pure function lattice_qDisorientation(Q1, Q2, struct)
use prec, only: &
tol_math_check
use IO, only: &
IO_error
tol_math_check
use math, only: &
math_qMul, &
math_qConj
implicit none
real(pReal), dimension(4) :: lattice_qDisorientation
real(pReal), dimension(4), intent(in) :: Q1, & ! 1st orientation
real(pReal), dimension(4), intent(in) :: &
Q1, & ! 1st orientation
Q2 ! 2nd orientation
integer(pInt), intent(in) :: symmetryType ! Type of crystal symmetry; 1:cubic, 2:hexagonal
! integer(kind(LATTICE_undefined_ID)), optional, intent(in) :: & ! if given, symmetries between the two orientation will be considered
! struct
integer(kind(LATTICE_undefined_ID)), optional, intent(in) :: & ! if given, symmetries between the two orientation will be considered
struct
real(pReal), dimension(4) :: dQ,dQsymA,mis
integer(pInt) :: i,j,k,s
integer(pInt) :: i,j,k,s,symmetry
integer(kind(LATTICE_undefined_ID)) :: myStruct
!--------------------------------------------------------------------------------------------------
! check if a structure with known symmetries is given
if (present(struct)) then
myStruct = struct
select case (struct)
case(LATTICE_fcc_ID,LATTICE_bcc_ID)
symmetry = 1_pInt
case(LATTICE_hex_ID)
symmetry = 2_pInt
case default
symmetry = 0_pInt
end select
else
symmetry = 0_pInt
myStruct = LATTICE_undefined_ID
endif
!--------------------------------------------------------------------------------------------------
! calculate misorientation, for cubic(1) and hexagonal(2) structure find symmetries
dQ = math_qMul(math_qConj(Q1),Q2)
lattice_qDisorientation = dQ
select case (symmetryType)
case (0_pInt)
if (lattice_qDisorientation(1) < 0.0_pReal) &
lattice_qDisorientation = -lattice_qDisorientation ! keep omega within 0 to 180 deg
select case(symmetry)
case (1_pInt,2_pInt)
s = sum(lattice_NsymOperations(1:symmetryType-1_pInt))
s = sum(lattice_NsymOperations(1:symmetry-1_pInt))
do i = 1_pInt,2_pInt
dQ = math_qConj(dQ) ! switch order of "from -- to"
do j = 1_pInt,lattice_NsymOperations(symmetryType) ! run through first crystal's symmetries
do j = 1_pInt,lattice_NsymOperations(symmetry) ! run through first crystal's symmetries
dQsymA = math_qMul(lattice_symOperations(1:4,s+j),dQ) ! apply sym
do k = 1_pInt,lattice_NsymOperations(symmetryType) ! run through 2nd crystal's symmetries
do k = 1_pInt,lattice_NsymOperations(symmetry) ! run through 2nd crystal's symmetries
mis = math_qMul(dQsymA,lattice_symOperations(1:4,s+k)) ! apply sym
if (mis(1) < 0.0_pReal) & ! want positive angle
mis = -mis
if (mis(1)-lattice_qDisorientation(1) > -tol_math_check .and. &
lattice_qInSST(mis,symmetryType)) &
lattice_qDisorientation = mis ! found better one
if (mis(1)-lattice_qDisorientation(1) > -tol_math_check &
.and. lattice_qInSST(mis,LATTICE_undefined_ID)) lattice_qDisorientation = mis ! found better one
enddo; enddo; enddo
case default
call IO_error(450_pInt,symmetryType) ! complain about unknown symmetry
case (0_pInt)
if (lattice_qDisorientation(1) < 0.0_pReal) lattice_qDisorientation = -lattice_qDisorientation ! keep omega within 0 to 180 deg
end select
end function lattice_qDisorientation
!--------------------------------------------------------------------------------------------------
!> @brief Number of parameters to expect in material.config section
! NslipFamilies
! NtwinFamilies
! SlipSlipInteraction
! SlipTwinInteraction
! TwinSlipInteraction
! TwinTwinInteraction
! NnonSchmid
!--------------------------------------------------------------------------------------------------
function lattice_configNchunks(struct_ID)
use prec, only: &
pInt
implicit none
integer(pInt), dimension(7) :: lattice_configNchunks
integer(kind(LATTICE_fcc_ID)) :: struct_ID
select case(struct_ID)
case (LATTICE_fcc_ID)
lattice_configNchunks(1) = count(lattice_fcc_NslipSystem > 0_pInt)
lattice_configNchunks(2) = count(lattice_fcc_NtwinSystem > 0_pInt)
lattice_configNchunks(3) = maxval(lattice_fcc_interactionSlipSlip)
lattice_configNchunks(4) = maxval(lattice_fcc_interactionSlipTwin)
lattice_configNchunks(5) = maxval(lattice_fcc_interactionTwinSlip)
lattice_configNchunks(6) = maxval(lattice_fcc_interactionTwinTwin)
lattice_configNchunks(7) = lattice_fcc_NnonSchmid
case (LATTICE_bcc_ID)
lattice_configNchunks(1) = count(lattice_bcc_NslipSystem > 0_pInt)
lattice_configNchunks(2) = count(lattice_bcc_NtwinSystem > 0_pInt)
lattice_configNchunks(3) = maxval(lattice_bcc_interactionSlipSlip)
lattice_configNchunks(4) = maxval(lattice_bcc_interactionSlipTwin)
lattice_configNchunks(5) = maxval(lattice_bcc_interactionTwinSlip)
lattice_configNchunks(6) = maxval(lattice_bcc_interactionTwinTwin)
lattice_configNchunks(7) = lattice_bcc_NnonSchmid
case (LATTICE_hex_ID)
lattice_configNchunks(1) = count(lattice_hex_NslipSystem > 0_pInt)
lattice_configNchunks(2) = count(lattice_hex_NtwinSystem > 0_pInt)
lattice_configNchunks(3) = maxval(lattice_hex_interactionSlipSlip)
lattice_configNchunks(4) = maxval(lattice_hex_interactionSlipTwin)
lattice_configNchunks(5) = maxval(lattice_hex_interactionTwinSlip)
lattice_configNchunks(6) = maxval(lattice_hex_interactionTwinTwin)
lattice_configNchunks(7) = lattice_hex_NnonSchmid
end select
end function lattice_configNchunks
end module lattice