sanity checks + documentation
This commit is contained in:
Martin Diehl 2018-12-12 00:29:19 +01:00
parent bf2b074787
commit 1446e9f4ab
2 changed files with 199 additions and 201 deletions

View File

@ -1236,6 +1236,10 @@ subroutine IO_error(error_ID,el,ip,g,instance,ext_msg)
msg = 'zero entry on stiffness diagonal' msg = 'zero entry on stiffness diagonal'
case (136_pInt) case (136_pInt)
msg = 'zero entry on stiffness diagonal for transformed phase' msg = 'zero entry on stiffness diagonal for transformed phase'
case (137_pInt)
msg = 'not defined for lattice structure'
case (138_pInt)
msg = 'not enough interaction parameters given'
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! errors related to the parsing of material.config ! errors related to the parsing of material.config

View File

@ -3,8 +3,8 @@
!> @author Philip Eisenlohr, 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 Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
!> @author Martin Diehl, 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 !> @brief contains lattice structure definitions including Schmid matrices for slip, twin, trans,
!> calculation and non-Schmid behavior ! and cleavage as well as interaction among the various systems
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module lattice module lattice
use prec, only: & use prec, only: &
@ -1395,6 +1395,7 @@ end subroutine lattice_initializeStructure
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Symmetrizes stiffness matrix according to lattice type !> @brief Symmetrizes stiffness matrix according to lattice type
!> @details J. A. Rayne and B. S. Chandrasekhar Phys. Rev. 120, 1658 Erratum Phys. Rev. 122, 1962
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function lattice_symmetrizeC66(struct,C66) pure function lattice_symmetrizeC66(struct,C66)
@ -1457,7 +1458,7 @@ pure function lattice_symmetrizeC66(struct,C66)
lattice_symmetrizeC66(3,2) = C66(1,3) lattice_symmetrizeC66(3,2) = C66(1,3)
lattice_symmetrizeC66(4,4) = C66(4,4) lattice_symmetrizeC66(4,4) = C66(4,4)
lattice_symmetrizeC66(5,5) = C66(4,4) lattice_symmetrizeC66(5,5) = C66(4,4)
lattice_symmetrizeC66(6,6) = C66(6,6) !J. A. Rayne and B. S. Chandrasekhar Phys. Rev. 120, 1658 Erratum Phys. Rev. 122, 1962 lattice_symmetrizeC66(6,6) = C66(6,6)
case default case default
lattice_symmetrizeC66 = C66 lattice_symmetrizeC66 = C66
end select end select
@ -1643,7 +1644,7 @@ end function lattice_qDisorientation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Provides characteristtic shear for twinning !> @brief Characteristic shear for twinning
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_characteristicShear_Twin(Ntwin,structure,CoverA) result(characteristicShear) function lattice_characteristicShear_Twin(Ntwin,structure,CoverA) result(characteristicShear)
use IO, only: & use IO, only: &
@ -1651,22 +1652,15 @@ function lattice_characteristicShear_Twin(Ntwin,structure,CoverA) result(charact
implicit none implicit none
integer(pInt), dimension(:), intent(in) :: Ntwin !< number of active twin systems per family integer(pInt), dimension(:), intent(in) :: Ntwin !< number of active twin systems per family
character(len=3), intent(in) :: structure character(len=3), intent(in) :: structure !< lattice structure
real(pReal), intent(in), optional :: & real(pReal), intent(in) :: cOverA !< c/a ratio
cOverA
real(pReal), dimension(sum(Ntwin)) :: characteristicShear real(pReal), dimension(sum(Ntwin)) :: characteristicShear
integer(pInt) :: & integer(pInt) :: &
ir, & !< index in reduced list a, & !< index of active system
ig, & !< index in full list c, & !< index in complete system list
mf, & !< index of my family mf, & !< index of my family
ms !< index of my system in current family ms !< index of my system in current family
real(pReal), dimension(LATTICE_FCC_NTWIN), parameter :: &
FCC_SHEARTWIN = 0.5_pReal*sqrt(2.0_pReal)
real(pReal), dimension(LATTICE_BCC_NTWIN), parameter :: &
BCC_SHEARTWIN = 0.5_pReal*sqrt(2.0_pReal)
integer(pInt), dimension(LATTICE_HEX_NTWIN), parameter :: & integer(pInt), dimension(LATTICE_HEX_NTWIN), parameter :: &
HEX_SHEARTWIN = reshape(int( [& HEX_SHEARTWIN = reshape(int( [&
1, & ! <-10.1>{10.2} 1, & ! <-10.1>{10.2}
@ -1693,32 +1687,31 @@ function lattice_characteristicShear_Twin(Ntwin,structure,CoverA) result(charact
4, & 4, &
4, & 4, &
4 & 4 &
],pInt),[LATTICE_HEX_NTWIN]) ! indicator to formula further below ],pInt),[LATTICE_HEX_NTWIN]) ! indicator to formulas below
ir = 0_pInt a = 0_pInt
myFamilies: do mf = 1_pInt,size(Ntwin,1) myFamilies: do mf = 1_pInt,size(Ntwin,1)
mySystems: do ms = 1_pInt,Ntwin(mf) mySystems: do ms = 1_pInt,Ntwin(mf)
ir = ir + 1_pInt a = a + 1_pInt
select case(structure) select case(trim(structure))
case('fcc') case('fcc','bcc')
ig = sum(LATTICE_FCC_NTWINSYSTEM(1:mf-1))+ms characteristicShear(a) = 0.5_pReal*sqrt(2.0_pReal)
characteristicShear(ir) = FCC_SHEARTWIN(ig)
case('bcc')
ig = sum(LATTICE_BCC_NTWINSYSTEM(1:mf-1))+ms
characteristicShear(ir) = BCC_SHEARTWIN(ig)
case('hex') case('hex')
if (.not. present(CoverA)) call IO_error(0_pInt) if (cOverA < 1.0_pReal .or. cOverA > 2.0_pReal) &
ig = sum(LATTICE_HEX_NTWINSYSTEM(1:mf-1))+ms call IO_error(131_pInt,ext_msg='lattice_characteristicShear_Twin')
select case(HEX_SHEARTWIN(ig)) ! from Christian & Mahajan 1995 p.29 c = sum(LATTICE_HEX_NTWINSYSTEM(1:mf-1))+ms
select case(HEX_SHEARTWIN(c)) ! from Christian & Mahajan 1995 p.29
case (1_pInt) ! <-10.1>{10.2} case (1_pInt) ! <-10.1>{10.2}
characteristicShear(ir) = (3.0_pReal-cOverA*cOverA)/sqrt(3.0_pReal)/CoverA characteristicShear(a) = (3.0_pReal-cOverA**2.0_pReal)/sqrt(3.0_pReal)/CoverA
case (2_pInt) ! <11.6>{-1-1.1} case (2_pInt) ! <11.6>{-1-1.1}
characteristicShear(ir) = 1.0_pReal/cOverA characteristicShear(a) = 1.0_pReal/cOverA
case (3_pInt) ! <10.-2>{10.1} case (3_pInt) ! <10.-2>{10.1}
characteristicShear(ir) = (4.0_pReal*cOverA*cOverA-9.0_pReal)/sqrt(48.0_pReal)/cOverA characteristicShear(a) = (4.0_pReal*cOverA**2.0_pReal-9.0_pReal)/sqrt(48.0_pReal)/cOverA
case (4_pInt) ! <11.-3>{11.2} case (4_pInt) ! <11.-3>{11.2}
characteristicShear(ir) = 2.0_pReal*(cOverA*cOverA-2.0_pReal)/3.0_pReal/cOverA characteristicShear(a) = 2.0_pReal*(cOverA**2.0_pReal-2.0_pReal)/3.0_pReal/cOverA
end select end select
case default
call IO_error(137_pInt,ext_msg='lattice_characteristicShear_Twin: '//trim(structure))
end select end select
enddo mySystems enddo mySystems
enddo myFamilies enddo myFamilies
@ -1727,7 +1720,7 @@ end function lattice_characteristicShear_Twin
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Calculates rotated elasticity matrices for twinning !> @brief Rotated elasticity matrices for twinning in Mandel notation
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_C66_twin(Ntwin,C66,structure,CoverA) function lattice_C66_twin(Ntwin,C66,structure,CoverA)
use IO, only: & use IO, only: &
@ -1742,8 +1735,8 @@ function lattice_C66_twin(Ntwin,C66,structure,CoverA)
implicit none implicit none
integer(pInt), dimension(:), intent(in) :: Ntwin !< number of active twin systems per family integer(pInt), dimension(:), intent(in) :: Ntwin !< number of active twin systems per family
character(len=*), intent(in) :: structure !< lattice structure character(len=*), intent(in) :: structure !< lattice structure
real(pReal), dimension(6,6), intent(in) :: C66 real(pReal), dimension(6,6), intent(in) :: C66 !< unrotated parent stiffness matrix
real(pReal), intent(in) :: cOverA real(pReal), intent(in) :: cOverA !< c/a ratio
real(pReal), dimension(6,6,sum(Ntwin)) :: lattice_C66_twin real(pReal), dimension(6,6,sum(Ntwin)) :: lattice_C66_twin
real(pReal), dimension(3,3,sum(Ntwin)) :: coordinateSystem real(pReal), dimension(3,3,sum(Ntwin)) :: coordinateSystem
@ -1751,16 +1744,20 @@ function lattice_C66_twin(Ntwin,C66,structure,CoverA)
real(pReal), dimension(3,3) :: R real(pReal), dimension(3,3) :: R
integer(pInt) :: i integer(pInt) :: i
select case(structure) select case(trim(structure))
case('fcc') case('fcc')
coordinateSystem = buildCoordinateSystem(Ntwin,LATTICE_FCC_NSLIPSYSTEM,LATTICE_FCC_SYSTEMTWIN,structure,cOverA) coordinateSystem = buildCoordinateSystem(Ntwin,LATTICE_FCC_NSLIPSYSTEM,LATTICE_FCC_SYSTEMTWIN,&
trim(structure),0.0_pReal)
case('bcc') case('bcc')
coordinateSystem = buildCoordinateSystem(Ntwin,LATTICE_BCC_NSLIPSYSTEM,LATTICE_BCC_SYSTEMTWIN,structure,cOverA) coordinateSystem = buildCoordinateSystem(Ntwin,LATTICE_BCC_NSLIPSYSTEM,LATTICE_BCC_SYSTEMTWIN,&
trim(structure),0.0_pReal)
case('hex','hexagonal') !ToDo: "No alias policy": long or short? case('hex','hexagonal') !ToDo: "No alias policy": long or short?
coordinateSystem = buildCoordinateSystem(Ntwin,LATTICE_HEX_NSLIPSYSTEM,LATTICE_HEX_SYSTEMTWIN,'hex',cOverA) coordinateSystem = buildCoordinateSystem(Ntwin,LATTICE_HEX_NSLIPSYSTEM,LATTICE_HEX_SYSTEMTWIN,&
'hex',cOverA)
case default case default
call IO_error(130_pInt,ext_msg=trim(structure)//' (lattice_C66_twin)') call IO_error(137_pInt,ext_msg='lattice_C66_twin: '//trim(structure))
end select end select
do i = 1, sum(Ntwin) do i = 1, sum(Ntwin)
R = math_axisAngleToR(coordinateSystem(1:3,2,i), 180.0_pReal * INRAD) ! ToDo: Why always 180 deg? R = math_axisAngleToR(coordinateSystem(1:3,2,i), 180.0_pReal * INRAD) ! ToDo: Why always 180 deg?
lattice_C66_twin(1:6,1:6,i) = math_Mandel3333to66(math_rotate_forward3333(math_Mandel66to3333(C66),R)) lattice_C66_twin(1:6,1:6,i) = math_Mandel3333to66(math_rotate_forward3333(math_Mandel66to3333(C66),R))
@ -1769,8 +1766,8 @@ end function lattice_C66_twin
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Calculates rotated elasticity matrices for transformation !> @brief Rotated elasticity matrices for transformation in Mandel notation
!> ToDo: Completely untested and incomplete !> ToDo: Completely untested and incomplete and undocumented
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_C66_trans(Ntrans,C_parent66,structure_parent, & function lattice_C66_trans(Ntrans,C_parent66,structure_parent, &
C_target66,structure_target, & C_target66,structure_target, &
@ -1840,16 +1837,14 @@ lattice_C66_trans = 0.0_pReal
end function end function
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Non-schmid tensor !> @brief Non-schmid projections for bcc with up to 6 coefficients
!> ToDo: Clean description needed ! Koester et al. 2012, Acta Materialia 60 (2012) 38943901, eq. (17)
! Schmid matrices with non-Schmid contributions according to Koester_etal2012, Acta Materialia 60 (2012) ! Gröger et al. 2008, Acta Materialia 56 (2008) 54125425, table 1
! 38943901, eq. (17) ("n1" is replaced by either "np" or "nn" according to either positive or negative slip direction)
! "np" and "nn" according to Gröger_etal2008, Acta Materialia 56 (2008) 54125425, table 1
! (corresponds to their "n1" for positive and negative slip direction respectively)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSchmidMatrix) function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSchmidMatrix)
use IO, only: &
IO_error
use math, only: & use math, only: &
INRAD, & INRAD, &
math_tensorproduct33, & math_tensorproduct33, &
@ -1858,20 +1853,21 @@ function lattice_nonSchmidMatrix(Nslip,nonSchmidCoefficients,sense) result(nonSc
math_axisAngleToR math_axisAngleToR
implicit none implicit none
integer(pInt), dimension(:), intent(in) :: Nslip !< number of active slip systems per family integer(pInt), dimension(:), intent(in) :: Nslip !< number of active slip systems per family
real(pReal), dimension(:), intent(in) :: nonSchmidCoefficients real(pReal), dimension(:), intent(in) :: nonSchmidCoefficients !< non-Schmid coefficients for projections
integer(pInt), intent(in) :: sense !< sense (-1,+1) integer(pInt), intent(in) :: sense !< sense (-1,+1)
real(pReal), dimension(1:3,1:3,sum(Nslip)) :: nonSchmidMatrix real(pReal), dimension(1:3,1:3,sum(Nslip)) :: nonSchmidMatrix
real(pReal), dimension(1:3,1:3,sum(Nslip)) :: coordinateSystem real(pReal), dimension(1:3,1:3,sum(Nslip)) :: coordinateSystem !< coordinate system of slip system
real(pReal), dimension(:), allocatable :: direction real(pReal), dimension(:), allocatable :: &
real(pReal), dimension(:), allocatable :: normal,np direction, normal, np
integer(pInt) :: i integer(pInt) :: i
if (abs(sense) /= 1_pInt) write(6,*) 'mist' if (abs(sense) /= 1_pInt) call IO_error(0_pInt,ext_msg='lattice_nonSchmidMatrix')
coordinateSystem = buildCoordinateSystem(Nslip,LATTICE_BCC_NSLIPSYSTEM,LATTICE_BCC_SYSTEMSLIP,'bcc',0.0_pReal)
coordinateSystem(1:3,1,1:sum(Nslip)) = coordinateSystem(1:3,1,1:sum(Nslip)) *real(sense,pReal) coordinateSystem = buildCoordinateSystem(Nslip,LATTICE_BCC_NSLIPSYSTEM,LATTICE_BCC_SYSTEMSLIP,&
nonSchmidMatrix = lattice_SchmidMatrix_slip(Nslip,'bcc',0.0_pReal) 'bcc',0.0_pReal)
coordinateSystem(1:3,1,1:sum(Nslip)) = coordinateSystem(1:3,1,1:sum(Nslip)) *real(sense,pReal) ! convert unidirectional coordinate system
nonSchmidMatrix = lattice_SchmidMatrix_slip(Nslip,'bcc',0.0_pReal) ! Schmid contribution
do i = 1_pInt,sum(Nslip) do i = 1_pInt,sum(Nslip)
direction = coordinateSystem(1:3,1,i) direction = coordinateSystem(1:3,1,i)
@ -1895,8 +1891,8 @@ end function lattice_nonSchmidMatrix
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Populates slip-slip interaction matrix !> @brief Slip-slip interaction matrix
!> details: only active slip systems are considered !> details only active slip systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_interaction_SlipSlip(Nslip,interactionValues,structure) result(interactionMatrix) function lattice_interaction_SlipSlip(Nslip,interactionValues,structure) result(interactionMatrix)
use IO, only: & use IO, only: &
@ -1904,7 +1900,7 @@ function lattice_interaction_SlipSlip(Nslip,interactionValues,structure) result(
implicit none implicit none
integer(pInt), dimension(:), intent(in) :: Nslip !< number of active slip systems per family integer(pInt), dimension(:), intent(in) :: Nslip !< number of active slip systems per family
real(pReal), dimension(:), intent(in) :: interactionValues !< interaction values slip-slip real(pReal), dimension(:), intent(in) :: interactionValues !< values for slip-slip interaction
character(len=*), intent(in) :: structure !< lattice structure character(len=*), intent(in) :: structure !< lattice structure
real(pReal), dimension(sum(Nslip),sum(Nslip)) :: interactionMatrix real(pReal), dimension(sum(Nslip),sum(Nslip)) :: interactionMatrix
@ -1925,20 +1921,17 @@ function lattice_interaction_SlipSlip(Nslip,interactionValues,structure) result(
interactionTypes = LATTICE_BCT_INTERACTIONSLIPSLIP interactionTypes = LATTICE_BCT_INTERACTIONSLIPSLIP
NslipMax = LATTICE_BCT_NSLIPSYSTEM NslipMax = LATTICE_BCT_NSLIPSYSTEM
case default case default
call IO_error(132_pInt,ext_msg=trim(structure)//' (slip slip interaction)') call IO_error(137_pInt,ext_msg='lattice_interaction_SlipSlip: '//trim(structure))
end select end select
!if (size(interactionValues) > maxval(interactionTypes)) &
! call IO_error(0_pInt) ! ToDo
interactionMatrix = buildInteraction(Nslip,Nslip,NslipMax,NslipMax,interactionValues,interactionTypes) interactionMatrix = buildInteraction(Nslip,Nslip,NslipMax,NslipMax,interactionValues,interactionTypes)
end function lattice_interaction_SlipSlip end function lattice_interaction_SlipSlip
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Populates twin-twin interaction matrix !> @brief Twin-twin interaction matrix
!> details: only active twin systems are considered !> details only active twin systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_interaction_TwinTwin(Ntwin,interactionValues,structure) result(interactionMatrix) function lattice_interaction_TwinTwin(Ntwin,interactionValues,structure) result(interactionMatrix)
use IO, only: & use IO, only: &
@ -1946,7 +1939,7 @@ function lattice_interaction_TwinTwin(Ntwin,interactionValues,structure) result(
implicit none implicit none
integer(pInt), dimension(:), intent(in) :: Ntwin !< number of active twin systems per family integer(pInt), dimension(:), intent(in) :: Ntwin !< number of active twin systems per family
real(pReal), dimension(:), intent(in) :: interactionValues !< interaction values twin-twin real(pReal), dimension(:), intent(in) :: interactionValues !< values for twin-twin interaction
character(len=*), intent(in) :: structure !< lattice structure character(len=*), intent(in) :: structure !< lattice structure
real(pReal), dimension(sum(Ntwin),sum(Ntwin)) :: interactionMatrix real(pReal), dimension(sum(Ntwin),sum(Ntwin)) :: interactionMatrix
@ -1967,7 +1960,7 @@ function lattice_interaction_TwinTwin(Ntwin,interactionValues,structure) result(
2,2,2,2,2,2,2,2,2,1,1,1, & 2,2,2,2,2,2,2,2,2,1,1,1, &
2,2,2,2,2,2,2,2,2,1,1,1, & 2,2,2,2,2,2,2,2,2,1,1,1, &
2,2,2,2,2,2,2,2,2,1,1,1 & 2,2,2,2,2,2,2,2,2,1,1,1 &
],pInt),shape(FCC_INTERACTIONTWINTWIN),order=[2,1]) !< Twin--twin interaction types for fcc ],pInt),shape(FCC_INTERACTIONTWINTWIN),order=[2,1]) !< Twin-twin interaction types for fcc
integer(pInt), dimension(LATTICE_BCC_NTWIN,LATTICE_BCC_NTWIN), parameter :: & integer(pInt), dimension(LATTICE_BCC_NTWIN,LATTICE_BCC_NTWIN), parameter :: &
BCC_INTERACTIONTWINTWIN = reshape(int( [& BCC_INTERACTIONTWINTWIN = reshape(int( [&
@ -1983,7 +1976,7 @@ function lattice_interaction_TwinTwin(Ntwin,interactionValues,structure) result(
3,3,3,2,2,3,3,3,3,1,3,3, & 3,3,3,2,2,3,3,3,3,1,3,3, &
2,3,3,3,3,3,3,2,3,3,1,3, & 2,3,3,3,3,3,3,2,3,3,1,3, &
3,2,3,3,3,3,2,3,3,3,3,1 & 3,2,3,3,3,3,2,3,3,3,3,1 &
],pInt),shape(BCC_INTERACTIONTWINTWIN),order=[2,1]) !< Twin--twin interaction types for bcc ],pInt),shape(BCC_INTERACTIONTWINTWIN),order=[2,1]) !< Twin-twin interaction types for bcc
!< 1: self interaction !< 1: self interaction
!< 2: collinear interaction !< 2: collinear interaction
!< 3: other interaction !< 3: other interaction
@ -2016,7 +2009,7 @@ function lattice_interaction_TwinTwin(Ntwin,interactionValues,structure) result(
20,20,20,20,20,20, 19,19,19,19,19,19, 18,18,18,18,18,18, 17,17,17,16,17,17, & 20,20,20,20,20,20, 19,19,19,19,19,19, 18,18,18,18,18,18, 17,17,17,16,17,17, &
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,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 & 20,20,20,20,20,20, 19,19,19,19,19,19, 18,18,18,18,18,18, 17,17,17,17,17,16 &
],pInt),shape(HEX_INTERACTIONTWINTWIN),order=[2,1]) !< Twin--twin interaction types for hex (isotropic, 16 in total) ],pInt),shape(HEX_INTERACTIONTWINTWIN),order=[2,1]) !< Twin-twin interaction types for hex
select case(structure) select case(structure)
case('fcc') case('fcc')
@ -2029,30 +2022,26 @@ function lattice_interaction_TwinTwin(Ntwin,interactionValues,structure) result(
interactionTypes = HEX_INTERACTIONTWINTWIN interactionTypes = HEX_INTERACTIONTWINTWIN
NtwinMax = LATTICE_HEX_NTWINSYSTEM NtwinMax = LATTICE_HEX_NTWINSYSTEM
case default case default
call IO_error(132_pInt,ext_msg=trim(structure)//' (twin twin interaction)') call IO_error(137_pInt,ext_msg='lattice_interaction_TwinTwin: '//trim(structure))
end select end select
!if (size(interactionValues) > maxval(interactionTypes)) &
! call IO_error(0_pInt) ! ToDo
interactionMatrix = buildInteraction(Ntwin,Ntwin,NtwinMax,NtwinMax,interactionValues,interactionTypes) interactionMatrix = buildInteraction(Ntwin,Ntwin,NtwinMax,NtwinMax,interactionValues,interactionTypes)
end function lattice_interaction_TwinTwin end function lattice_interaction_TwinTwin
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Populates trans-trans interaction matrix !> @brief Trans-trans interaction matrix
!> details: only active transformation systems are considered !> details only active trans systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_interaction_TransTrans(Ntrans,interactionValues,structure) result(interactionMatrix) function lattice_interaction_TransTrans(Ntrans,interactionValues,structure) result(interactionMatrix)
use IO, only: & use IO, only: &
IO_error IO_error
implicit none implicit none
integer(pInt), dimension(:), intent(in) :: Ntrans !< number of active twin systems per family integer(pInt), dimension(:), intent(in) :: Ntrans !< number of active trans systems per family
real(pReal), dimension(:), intent(in) :: interactionValues !< interaction values twin-twin real(pReal), dimension(:), intent(in) :: interactionValues !< values for trans-trans interaction
character(len=*), intent(in) :: & character(len=*), intent(in) :: structure !< lattice structure (parent crystal)
structure !< lattice structure of parent crystal
real(pReal), dimension(sum(Ntrans),sum(Ntrans)) :: interactionMatrix real(pReal), dimension(sum(Ntrans),sum(Ntrans)) :: interactionMatrix
integer(pInt), dimension(:), allocatable :: NtransMax integer(pInt), dimension(:), allocatable :: NtransMax
@ -2072,24 +2061,23 @@ function lattice_interaction_TransTrans(Ntrans,interactionValues,structure) resu
2,2,2,2,2,2,2,2,2,1,1,1, & 2,2,2,2,2,2,2,2,2,1,1,1, &
2,2,2,2,2,2,2,2,2,1,1,1, & 2,2,2,2,2,2,2,2,2,1,1,1, &
2,2,2,2,2,2,2,2,2,1,1,1 & 2,2,2,2,2,2,2,2,2,1,1,1 &
],pInt),shape(FCC_INTERACTIONTRANSTRANS),order=[2,1]) !< Trans--trans interaction types for fcc ],pInt),shape(FCC_INTERACTIONTRANSTRANS),order=[2,1]) !< Trans-trans interaction types for fcc
if (trim(structure) == 'fcc') then if (trim(structure) == 'fcc') then
interactionTypes = FCC_INTERACTIONTRANSTRANS interactionTypes = FCC_INTERACTIONTRANSTRANS
NtransMax = LATTICE_FCC_NTRANSSYSTEM NtransMax = LATTICE_FCC_NTRANSSYSTEM
else else
call IO_error(132_pInt,ext_msg=trim(structure)//' (trans trans interaction)') call IO_error(137_pInt,ext_msg='lattice_interaction_TransTrans: '//trim(structure))
end if end if
!if (size(interactionValues) > maxval(interactionTypes)) &
! call IO_error(0_pInt) ! ToDo
interactionMatrix = buildInteraction(Ntrans,Ntrans,NtransMax,NtransMax,interactionValues,interactionTypes) interactionMatrix = buildInteraction(Ntrans,Ntrans,NtransMax,NtransMax,interactionValues,interactionTypes)
end function lattice_interaction_TransTrans end function lattice_interaction_TransTrans
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Populates slip-twin interaction matrix !> @brief Slip-twin interaction matrix
!> details: only active slip and twin systems are considered !> details only active slip and twin systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_interaction_SlipTwin(Nslip,Ntwin,interactionValues,structure) result(interactionMatrix) function lattice_interaction_SlipTwin(Nslip,Ntwin,interactionValues,structure) result(interactionMatrix)
use IO, only: & use IO, only: &
@ -2098,7 +2086,7 @@ function lattice_interaction_SlipTwin(Nslip,Ntwin,interactionValues,structure) r
implicit none implicit none
integer(pInt), dimension(:), intent(in) :: Nslip, & !< number of active slip systems per family integer(pInt), dimension(:), intent(in) :: Nslip, & !< number of active slip systems per family
Ntwin !< number of active twin systems per family Ntwin !< number of active twin systems per family
real(pReal), dimension(:), intent(in) :: interactionValues !< interaction values twin-twin real(pReal), dimension(:), intent(in) :: interactionValues !< values for slip-twin interaction
character(len=*), intent(in) :: structure !< lattice structure character(len=*), intent(in) :: structure !< lattice structure
real(pReal), dimension(sum(Nslip),sum(Ntwin)) :: interactionMatrix real(pReal), dimension(sum(Nslip),sum(Ntwin)) :: interactionMatrix
@ -2127,7 +2115,7 @@ function lattice_interaction_SlipTwin(Nslip,Ntwin,interactionValues,structure) r
4,4,4,4,4,4,4,4,4,4,4,4, & 4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4, & 4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4 & 4,4,4,4,4,4,4,4,4,4,4,4 &
],pInt),shape(FCC_INTERACTIONSLIPTWIN),order=[2,1]) !< Slip--twin interaction types for fcc ],pInt),shape(FCC_INTERACTIONSLIPTWIN),order=[2,1]) !< Slip-twin interaction types for fcc
!< 1: coplanar interaction !< 1: coplanar interaction
!< 2: screw trace between slip system and twin habit plane (easy cross slip) !< 2: screw trace between slip system and twin habit plane (easy cross slip)
!< 3: other interaction !< 3: other interaction
@ -2158,7 +2146,7 @@ function lattice_interaction_SlipTwin(Nslip,Ntwin,interactionValues,structure) r
3,3,3,2,2,3,3,3,3,1,3,3, & 3,3,3,2,2,3,3,3,3,1,3,3, &
2,3,3,3,3,3,3,2,3,3,1,3, & 2,3,3,3,3,3,3,2,3,3,1,3, &
3,2,3,3,3,3,2,3,3,3,3,1 & 3,2,3,3,3,3,2,3,3,3,3,1 &
],pInt),shape(BCC_INTERACTIONSLIPTWIN),order=[2,1]) !< Slip--twin interaction types for bcc ],pInt),shape(BCC_INTERACTIONSLIPTWIN),order=[2,1]) !< Slip-twin interaction types for bcc
!< 1: coplanar interaction !< 1: coplanar interaction
!< 2: screw trace between slip system and twin habit plane (easy cross slip) !< 2: screw trace between slip system and twin habit plane (easy cross slip)
!< 3: other interaction !< 3: other interaction
@ -2203,7 +2191,7 @@ function lattice_interaction_SlipTwin(Nslip,Ntwin,interactionValues,structure) r
21,21,21,21,21,21, 22,22,22,22,22,22, 23,23,23,23,23,23, 24,24,24,24,24,24, & 21,21,21,21,21,21, 22,22,22,22,22,22, 23,23,23,23,23,23, 24,24,24,24,24,24, &
21,21,21,21,21,21, 22,22,22,22,22,22, 23,23,23,23,23,23, 24,24,24,24,24,24 & 21,21,21,21,21,21, 22,22,22,22,22,22, 23,23,23,23,23,23, 24,24,24,24,24,24 &
! !
],pInt),shape(HEX_INTERACTIONSLIPTWIN),order=[2,1]) !< Slip--twin interaction types for hex (isotropic, 24 in total) ],pInt),shape(HEX_INTERACTIONSLIPTWIN),order=[2,1]) !< Slip-twin interaction types for hex
select case(structure) select case(structure)
@ -2220,20 +2208,17 @@ function lattice_interaction_SlipTwin(Nslip,Ntwin,interactionValues,structure) r
NslipMax = LATTICE_HEX_NSLIPSYSTEM NslipMax = LATTICE_HEX_NSLIPSYSTEM
NtwinMax = LATTICE_HEX_NTWINSYSTEM NtwinMax = LATTICE_HEX_NTWINSYSTEM
case default case default
call IO_error(132_pInt,ext_msg=trim(structure)//' (slip twin interaction)') call IO_error(137_pInt,ext_msg='lattice_interaction_SlipTwin: '//trim(structure))
end select end select
!if (size(interactionValues) > maxval(interactionTypes)) &
! call IO_error(0_pInt) ! ToDo
interactionMatrix = buildInteraction(Nslip,Ntwin,NslipMax,NtwinMax,interactionValues,interactionTypes) interactionMatrix = buildInteraction(Nslip,Ntwin,NslipMax,NtwinMax,interactionValues,interactionTypes)
end function lattice_interaction_SlipTwin end function lattice_interaction_SlipTwin
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Populates trans-trans interaction matrix !> @brief Slip-trans interaction matrix
!> details: only active transformation systems are considered !> details only active slip and trans systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_interaction_SlipTrans(Nslip,Ntrans,interactionValues,structure) result(interactionMatrix) function lattice_interaction_SlipTrans(Nslip,Ntrans,interactionValues,structure) result(interactionMatrix)
use IO, only: & use IO, only: &
@ -2242,9 +2227,9 @@ function lattice_interaction_SlipTrans(Nslip,Ntrans,interactionValues,structure)
implicit none implicit none
integer(pInt), dimension(:), intent(in) :: Nslip, & !< number of active slip systems per family integer(pInt), dimension(:), intent(in) :: Nslip, & !< number of active slip systems per family
Ntrans !< number of active trans systems per family Ntrans !< number of active trans systems per family
real(pReal), dimension(:), intent(in) :: interactionValues !< interaction values slip--trans real(pReal), dimension(:), intent(in) :: interactionValues !< values for slip-trans interaction
character(len=*), intent(in) :: & character(len=*), intent(in) :: &
structure !< lattice structure of parent crystal structure !< lattice structure (parent crystal)
real(pReal), dimension(sum(Nslip),sum(Ntrans)) :: interactionMatrix real(pReal), dimension(sum(Nslip),sum(Ntrans)) :: interactionMatrix
integer(pInt), dimension(:), allocatable :: NslipMax, & integer(pInt), dimension(:), allocatable :: NslipMax, &
@ -2272,7 +2257,7 @@ function lattice_interaction_SlipTrans(Nslip,Ntrans,interactionValues,structure)
4,4,4,4,4,4,4,4,4,4,4,4, & 4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4, & 4,4,4,4,4,4,4,4,4,4,4,4, &
4,4,4,4,4,4,4,4,4,4,4,4 & 4,4,4,4,4,4,4,4,4,4,4,4 &
],pInt),shape(FCC_INTERACTIONSLIPTRANS),order=[2,1]) !< Slip--trans interaction types for fcc ],pInt),shape(FCC_INTERACTIONSLIPTRANS),order=[2,1]) !< Slip-trans interaction types for fcc
select case(structure) select case(structure)
case('fcc') case('fcc')
@ -2280,19 +2265,17 @@ function lattice_interaction_SlipTrans(Nslip,Ntrans,interactionValues,structure)
NslipMax = LATTICE_FCC_NSLIPSYSTEM NslipMax = LATTICE_FCC_NSLIPSYSTEM
NtransMax = LATTICE_FCC_NTRANSSYSTEM NtransMax = LATTICE_FCC_NTRANSSYSTEM
case default case default
call IO_error(132_pInt,ext_msg=trim(structure)//' (slip trans interaction)') call IO_error(137_pInt,ext_msg='lattice_interaction_SlipTrans: '//trim(structure))
end select end select
!if (size(interactionValues) > maxval(interactionTypes)) &
! call IO_error(0_pInt) ! ToDo
interactionMatrix = buildInteraction(Nslip,Ntrans,NslipMax,NtransMax,interactionValues,interactionTypes) interactionMatrix = buildInteraction(Nslip,Ntrans,NslipMax,NtransMax,interactionValues,interactionTypes)
end function lattice_interaction_SlipTrans end function lattice_interaction_SlipTrans
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Populates twin-slip interaction matrix !> @brief Twin-slip interaction matrix
!> details: only active twin and slip systems are considered !> details only active twin and slip systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_interaction_TwinSlip(Ntwin,Nslip,interactionValues,structure) result(interactionMatrix) function lattice_interaction_TwinSlip(Ntwin,Nslip,interactionValues,structure) result(interactionMatrix)
use IO, only: & use IO, only: &
@ -2301,7 +2284,7 @@ function lattice_interaction_TwinSlip(Ntwin,Nslip,interactionValues,structure) r
implicit none implicit none
integer(pInt), dimension(:), intent(in) :: Ntwin, & !< number of active twin systems per family integer(pInt), dimension(:), intent(in) :: Ntwin, & !< number of active twin systems per family
Nslip !< number of active slip systems per family Nslip !< number of active slip systems per family
real(pReal), dimension(:), intent(in) :: interactionValues !< interaction values twin-twin real(pReal), dimension(:), intent(in) :: interactionValues !< values for twin-twin interaction
character(len=*), intent(in) :: structure !< lattice structure character(len=*), intent(in) :: structure !< lattice structure
real(pReal), dimension(sum(Ntwin),sum(Nslip)) :: interactionMatrix real(pReal), dimension(sum(Ntwin),sum(Nslip)) :: interactionMatrix
@ -2310,10 +2293,10 @@ function lattice_interaction_TwinSlip(Ntwin,Nslip,interactionValues,structure) r
integer(pInt), dimension(:,:), allocatable :: interactionTypes integer(pInt), dimension(:,:), allocatable :: interactionTypes
integer(pInt), dimension(LATTICE_FCC_NTWIN,LATTICE_FCC_NSLIP), parameter :: & integer(pInt), dimension(LATTICE_FCC_NTWIN,LATTICE_FCC_NSLIP), parameter :: &
FCC_INTERACTIONTWINSLIP = 1_pInt !< Twin--Slip interaction types for fcc FCC_INTERACTIONTWINSLIP = 1_pInt !< Twin-Slip interaction types for fcc
integer(pInt), dimension(LATTICE_BCC_NTWIN,LATTICE_BCC_NSLIP), parameter :: & integer(pInt), dimension(LATTICE_BCC_NTWIN,LATTICE_BCC_NSLIP), parameter :: &
BCC_INTERACTIONTWINSLIP = 1_pInt !< Twin--slip interaction types for bcc BCC_INTERACTIONTWINSLIP = 1_pInt !< Twin-slip interaction types for bcc
integer(pInt), dimension(LATTICE_HEX_NTWIN,LATTICE_HEX_NSLIP), parameter :: & integer(pInt), dimension(LATTICE_HEX_NTWIN,LATTICE_HEX_NSLIP), parameter :: &
HEX_INTERACTIONTWINSLIP = reshape(int( [& HEX_INTERACTIONTWINSLIP = reshape(int( [&
@ -2344,7 +2327,7 @@ function lattice_interaction_TwinSlip(Ntwin,Nslip,interactionValues,structure) r
4, 4, 4, 8, 8, 8, 12,12,12, 16,16,16,16,16,16, 20,20,20,20,20,20,20,20,20,20,20,20, 24,24,24,24,24,24, & 4, 4, 4, 8, 8, 8, 12,12,12, 16,16,16,16,16,16, 20,20,20,20,20,20,20,20,20,20,20,20, 24,24,24,24,24,24, &
4, 4, 4, 8, 8, 8, 12,12,12, 16,16,16,16,16,16, 20,20,20,20,20,20,20,20,20,20,20,20, 24,24,24,24,24,24, & 4, 4, 4, 8, 8, 8, 12,12,12, 16,16,16,16,16,16, 20,20,20,20,20,20,20,20,20,20,20,20, 24,24,24,24,24,24, &
4, 4, 4, 8, 8, 8, 12,12,12, 16,16,16,16,16,16, 20,20,20,20,20,20,20,20,20,20,20,20, 24,24,24,24,24,24 & 4, 4, 4, 8, 8, 8, 12,12,12, 16,16,16,16,16,16, 20,20,20,20,20,20,20,20,20,20,20,20, 24,24,24,24,24,24 &
],pInt),shape(HEX_INTERACTIONTWINSLIP),order=[2,1]) !< Twin--twin interaction types for hex (isotropic, 20 in total) ],pInt),shape(HEX_INTERACTIONTWINSLIP),order=[2,1]) !< Twin-twin interaction types for hex
select case(structure) select case(structure)
case('fcc') case('fcc')
@ -2360,19 +2343,17 @@ function lattice_interaction_TwinSlip(Ntwin,Nslip,interactionValues,structure) r
NtwinMax = LATTICE_HEX_NTWINSYSTEM NtwinMax = LATTICE_HEX_NTWINSYSTEM
NslipMax = LATTICE_HEX_NSLIPSYSTEM NslipMax = LATTICE_HEX_NSLIPSYSTEM
case default case default
call IO_error(132_pInt,ext_msg=trim(structure)//' (twin slip interaction)') call IO_error(137_pInt,ext_msg='lattice_interaction_TwinSlip: '//trim(structure))
end select end select
!if (size(interactionValues) > maxval(interactionTypes)) &
! call IO_error(0_pInt) ! ToDo
interactionMatrix = buildInteraction(Ntwin,Nslip,NtwinMax,NslipMax,interactionValues,interactionTypes) interactionMatrix = buildInteraction(Ntwin,Nslip,NtwinMax,NslipMax,interactionValues,interactionTypes)
end function lattice_interaction_TwinSlip end function lattice_interaction_TwinSlip
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Calculates Schmid matrix for active slip systems !> @brief Schmid matrix for slip
!> details only active slip systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_SchmidMatrix_slip(Nslip,structure,cOverA) result(SchmidMatrix) function lattice_SchmidMatrix_slip(Nslip,structure,cOverA) result(SchmidMatrix)
use prec, only: & use prec, only: &
@ -2408,7 +2389,7 @@ function lattice_SchmidMatrix_slip(Nslip,structure,cOverA) result(SchmidMatrix)
NslipMax = LATTICE_BCT_NSLIPSYSTEM NslipMax = LATTICE_BCT_NSLIPSYSTEM
slipSystems = LATTICE_BCT_SYSTEMSLIP slipSystems = LATTICE_BCT_SYSTEMSLIP
case default case default
call IO_error(130_pInt,ext_msg=trim(structure)//' (lattice_SchmidMatrix_slip)') call IO_error(137_pInt,ext_msg='lattice_SchmidMatrix_slip: '//trim(structure))
end select end select
if (any(NslipMax(1:size(Nslip)) - Nslip < 0_pInt)) & if (any(NslipMax(1:size(Nslip)) - Nslip < 0_pInt)) &
@ -2428,7 +2409,8 @@ end function lattice_SchmidMatrix_slip
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Calculates Schmid matrix for active twin systems !> @brief Schmid matrix for twinning
!> details only active twin systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_SchmidMatrix_twin(Ntwin,structure,cOverA) result(SchmidMatrix) function lattice_SchmidMatrix_twin(Ntwin,structure,cOverA) result(SchmidMatrix)
use prec, only: & use prec, only: &
@ -2442,7 +2424,7 @@ function lattice_SchmidMatrix_twin(Ntwin,structure,cOverA) result(SchmidMatrix)
implicit none implicit none
integer(pInt), dimension(:), intent(in) :: Ntwin !< number of active twin systems per family integer(pInt), dimension(:), intent(in) :: Ntwin !< number of active twin systems per family
character(len=*), intent(in) :: structure !< lattice structure character(len=*), intent(in) :: structure !< lattice structure
real(pReal), intent(in) :: cOverA real(pReal), intent(in) :: cOverA !< c/a ratio
real(pReal), dimension(3,3,sum(Ntwin)) :: SchmidMatrix real(pReal), dimension(3,3,sum(Ntwin)) :: SchmidMatrix
real(pReal), dimension(3,3,sum(Ntwin)) :: coordinateSystem real(pReal), dimension(3,3,sum(Ntwin)) :: coordinateSystem
@ -2461,7 +2443,7 @@ function lattice_SchmidMatrix_twin(Ntwin,structure,cOverA) result(SchmidMatrix)
NtwinMax = LATTICE_HEX_NTWINSYSTEM NtwinMax = LATTICE_HEX_NTWINSYSTEM
twinSystems = LATTICE_HEX_SYSTEMTWIN twinSystems = LATTICE_HEX_SYSTEMTWIN
case default case default
call IO_error(130_pInt,ext_msg=trim(structure)//' (lattice_SchmidMatrix_twin)') call IO_error(137_pInt,ext_msg='lattice_SchmidMatrix_twin: '//trim(structure))
end select end select
if (any(NtwinMax(1:size(Ntwin)) - Ntwin < 0_pInt)) & if (any(NtwinMax(1:size(Ntwin)) - Ntwin < 0_pInt)) &
@ -2481,7 +2463,8 @@ end function lattice_SchmidMatrix_twin
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Calculates Schmid matrix for active cleavage systems !> @brief Schmid matrix for cleavage
!> details only active cleavage systems are considered
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_SchmidMatrix_cleavage(Ncleavage,structure,cOverA) result(SchmidMatrix) function lattice_SchmidMatrix_cleavage(Ncleavage,structure,cOverA) result(SchmidMatrix)
use math, only: & use math, only: &
@ -2492,7 +2475,7 @@ function lattice_SchmidMatrix_cleavage(Ncleavage,structure,cOverA) result(Schmid
implicit none implicit none
integer(pInt), dimension(:), intent(in) :: Ncleavage !< number of active cleavage systems per family integer(pInt), dimension(:), intent(in) :: Ncleavage !< number of active cleavage systems per family
character(len=*), intent(in) :: structure !< lattice structure character(len=*), intent(in) :: structure !< lattice structure
real(pReal), intent(in) :: cOverA real(pReal), intent(in) :: cOverA !< c/a ratio
real(pReal), dimension(3,3,3,sum(Ncleavage)) :: SchmidMatrix real(pReal), dimension(3,3,3,sum(Ncleavage)) :: SchmidMatrix
real(pReal), dimension(3,3,sum(Ncleavage)) :: coordinateSystem real(pReal), dimension(3,3,sum(Ncleavage)) :: coordinateSystem
@ -2517,7 +2500,7 @@ function lattice_SchmidMatrix_cleavage(Ncleavage,structure,cOverA) result(Schmid
NcleavageMax = LATTICE_HEX_NCLEAVAGESYSTEM NcleavageMax = LATTICE_HEX_NCLEAVAGESYSTEM
cleavageSystems = LATTICE_HEX_SYSTEMCLEAVAGE cleavageSystems = LATTICE_HEX_SYSTEMCLEAVAGE
case default case default
call IO_error(130_pInt,ext_msg=trim(structure)//' (lattice_SchmidMatrix_cleavage)') call IO_error(137_pInt,ext_msg='lattice_SchmidMatrix_cleavage: '//trim(structure))
end select end select
if (any(NcleavageMax(1:size(Ncleavage)) - Ncleavage < 0_pInt)) & if (any(NcleavageMax(1:size(Ncleavage)) - Ncleavage < 0_pInt)) &
@ -2537,7 +2520,7 @@ end function lattice_SchmidMatrix_cleavage
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Calculates forest projection (for edge dislocations) !> @brief Forest projection (for edge dislocations)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function lattice_forestProjection(Nslip,structure,cOverA) result(projection) function lattice_forestProjection(Nslip,structure,cOverA) result(projection)
use math, only: & use math, only: &
@ -2548,7 +2531,7 @@ function lattice_forestProjection(Nslip,structure,cOverA) result(projection)
implicit none implicit none
integer(pInt), dimension(:), intent(in) :: Nslip !< number of active slip systems per family integer(pInt), dimension(:), intent(in) :: Nslip !< number of active slip systems per family
character(len=*), intent(in) :: structure !< lattice structure character(len=*), intent(in) :: structure !< lattice structure
real(pReal), intent(in) :: cOverA real(pReal), intent(in) :: cOverA !< c/a ratio
real(pReal), dimension(sum(Nslip),sum(Nslip)) :: projection real(pReal), dimension(sum(Nslip),sum(Nslip)) :: projection
real(pReal), dimension(3,3,sum(Nslip)) :: coordinateSystem real(pReal), dimension(3,3,sum(Nslip)) :: coordinateSystem
@ -2570,7 +2553,7 @@ function lattice_forestProjection(Nslip,structure,cOverA) result(projection)
NslipMax = LATTICE_BCT_NSLIPSYSTEM NslipMax = LATTICE_BCT_NSLIPSYSTEM
slipSystems = LATTICE_BCT_SYSTEMSLIP slipSystems = LATTICE_BCT_SYSTEMSLIP
case default case default
call IO_error(130_pInt,ext_msg=trim(structure)//' (lattice_forrestProjection)') call IO_error(137_pInt,ext_msg='lattice_forestProjection: '//trim(structure))
end select end select
if (any(NslipMax(1:size(Nslip)) - Nslip < 0_pInt)) & if (any(NslipMax(1:size(Nslip)) - Nslip < 0_pInt)) &
@ -2590,8 +2573,9 @@ end function lattice_forestProjection
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Populates reduced interaction matrix !> @brief Populates reduced interaction matrix
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
pure function buildInteraction(activeA,activeB,maxA,maxB,values,matrix) function buildInteraction(activeA,activeB,maxA,maxB,values,matrix)
use IO, only: &
IO_error
implicit none implicit none
integer(pInt), dimension(:), intent(in) :: & integer(pInt), dimension(:), intent(in) :: &
activeA, & !< number of active systems as specified in material.config activeA, & !< number of active systems as specified in material.config
@ -2599,7 +2583,7 @@ pure function buildInteraction(activeA,activeB,maxA,maxB,values,matrix)
maxA, & !< number of maximum available systems maxA, & !< number of maximum available systems
maxB !< number of maximum available systems maxB !< number of maximum available systems
real(pReal), dimension(:), intent(in) :: values !< interaction values real(pReal), dimension(:), intent(in) :: values !< interaction values
integer(pInt), dimension(:,:), intent(in) :: matrix !< full interaction matrix integer(pInt), dimension(:,:), intent(in) :: matrix !< complete interaction matrix
real(pReal), dimension(sum(activeA),sum(activeB)) :: buildInteraction real(pReal), dimension(sum(activeA),sum(activeB)) :: buildInteraction
integer(pInt) :: & integer(pInt) :: &
@ -2613,6 +2597,8 @@ pure function buildInteraction(activeA,activeB,maxA,maxB,values,matrix)
otherFamilies: do of = 1_pInt,size(activeB,1) otherFamilies: do of = 1_pInt,size(activeB,1)
index_otherFamily = sum(activeB(1:of-1_pInt)) index_otherFamily = sum(activeB(1:of-1_pInt))
otherSystems: do os = 1_pInt,activeB(of) otherSystems: do os = 1_pInt,activeB(of)
if(matrix(sum(maxA(1:mf-1))+ms, sum(maxB(1:of-1))+os) > size(values)) &
call IO_error(138,ext_msg='buildInteraction')
buildInteraction(index_myFamily+ms,index_otherFamily+os) = & buildInteraction(index_myFamily+ms,index_otherFamily+os) = &
values(matrix(sum(maxA(1:mf-1))+ms, sum(maxB(1:of-1))+os)) values(matrix(sum(maxA(1:mf-1))+ms, sum(maxB(1:of-1))+os))
enddo otherSystems; enddo otherFamilies; enddo otherSystems; enddo otherFamilies;
@ -2624,16 +2610,18 @@ end function buildInteraction
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief build a local coordinate system in a slip, twin, trans, cleavage system !> @brief build a local coordinate system in a slip, twin, trans, cleavage system
!> @details: Order: Direction, plane (normal), and common perpendicular !> @details Order: Direction, plane (normal), and common perpendicular
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function buildCoordinateSystem(active,maximum,system,structure,cOverA) function buildCoordinateSystem(active,complete,system,structure,cOverA)
use IO, only: &
IO_error
use math, only: & use math, only: &
math_crossproduct math_crossproduct
implicit none implicit none
integer(pInt), dimension(:), intent(in) :: & integer(pInt), dimension(:), intent(in) :: &
active, & active, &
maximum complete
real(pReal), dimension(:,:), intent(in) :: & real(pReal), dimension(:,:), intent(in) :: &
system system
character(len=*), intent(in) :: & character(len=*), intent(in) :: &
@ -2646,46 +2634,50 @@ function buildCoordinateSystem(active,maximum,system,structure,cOverA)
real(pReal), dimension(3) :: & real(pReal), dimension(3) :: &
direction, normal direction, normal
integer(pInt) :: & integer(pInt) :: &
i, & !< index in reduced matrix a, & !< index of active system
j, & !< index in full matrix c, & !< index in complete system matrix
f, & !< index of my family f, & !< index of my family
s !< index of my system in current family s !< index of my system in current family
i = 0_pInt a = 0_pInt
activeFamilies: do f = 1_pInt,size(active,1) activeFamilies: do f = 1_pInt,size(active,1)
activeSystems: do s = 1_pInt,active(f) activeSystems: do s = 1_pInt,active(f)
i = i + 1_pInt a = a + 1_pInt
j = sum(maximum(1:f-1))+s c = sum(complete(1:f-1))+s
select case(trim(structure)) select case(trim(structure))
case ('fcc','bcc') case ('fcc','bcc')
direction = system(1:3,j) direction = system(1:3,c)
normal = system(4:6,j) normal = system(4:6,c)
case ('hex') case ('hex')
!ToDo: check if c/a ratio is sensible if (cOverA < 1.0_pReal .or. cOverA > 2.0_pReal) &
! direction [uvtw]->[3u/2 (u+2v)*sqrt(3)/2 w*(c/a)]) call IO_error(131_pInt,ext_msg='buildCoordinateSystem:'//trim(structure))
direction = [ system(1,j)*1.5_pReal, &
(system(1,j)+2.0_pReal*system(2,j))*sqrt(0.75_pReal), &
system(4,j)*CoverA ]
! plane (hkil)->(h (h+2k)/sqrt(3) l/(c/a)) direction = [ system(1,c)*1.5_pReal, &
normal = [ system(5,j), & (system(1,c)+2.0_pReal*system(2,c))*sqrt(0.75_pReal), &
(system(5,j)+2.0_pReal*system(6,j))/ sqrt(3.0_pReal), & system(4,c)*cOverA ] ! direction [uvtw]->[3u/2 (u+2v)*sqrt(3)/2 w*(c/a)])
system(8,j)/CoverA ]
normal = [ system(5,c), &
(system(5,c)+2.0_pReal*system(6,c))/sqrt(3.0_pReal), &
system(8,c)/cOverA ] ! plane (hkil)->(h (h+2k)/sqrt(3) l/(c/a))
case ('bct') case ('bct')
!ToDo: check if c/a ratio is sensible if (cOverA > 2.0_pReal) &
direction = [system(1:2,j),system(3,i)*CoverA] call IO_error(131_pInt,ext_msg='buildCoordinateSystem:'//trim(structure))
normal = [system(4:5,j),system(6,i)/CoverA] direction = [system(1:2,c),system(3,c)*cOverA]
normal = [system(4:5,c),system(6,c)/cOverA]
case default
call IO_error(137_pInt,ext_msg='buildCoordinateSystem: '//trim(structure))
end select end select
buildCoordinateSystem(1:3,1,i) = direction/norm2(direction) buildCoordinateSystem(1:3,1,a) = direction/norm2(direction)
buildCoordinateSystem(1:3,2,i) = normal/norm2(normal) buildCoordinateSystem(1:3,2,a) = normal/norm2(normal)
buildCoordinateSystem(1:3,3,i) = math_crossproduct(buildCoordinateSystem(1:3,1,i),& buildCoordinateSystem(1:3,3,a) = math_crossproduct(buildCoordinateSystem(1:3,1,a),&
buildCoordinateSystem(1:3,2,i)) buildCoordinateSystem(1:3,2,a))
enddo activeSystems enddo activeSystems
enddo activeFamilies enddo activeFamilies
@ -2693,7 +2685,9 @@ function buildCoordinateSystem(active,maximum,system,structure,cOverA)
end function buildCoordinateSystem end function buildCoordinateSystem
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief xxx !> @brief Helper function to define transformation systems
! Needed for Schmid_trans + C66_trans
! ToDo: completely untested and uncommented
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine lattice_Trans(Q,S,Ntrans,cOverA,a_fcc,a_bcc) subroutine lattice_Trans(Q,S,Ntrans,cOverA,a_fcc,a_bcc)
use math, only: & use math, only: &