diff --git a/src/lattice.f90 b/src/lattice.f90 index d366674c7..3c012ac43 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -365,12 +365,6 @@ module lattice 1, 1, 1, 1,-2, 1 & ],pReal),shape(BCT_SYSTEMSLIP)) !< bct slip systems for c/a = 0.5456 (Sn), sorted by Bieler 2009 (https://doi.org/10.1007/s11664-009-0909-x) -! SHOULD NOT BE PART OF LATTICE BEGIN - real(pReal), dimension(:), allocatable, public, protected :: & - lattice_mu, lattice_nu - real(pReal), dimension(:,:,:), allocatable, public, protected :: & - lattice_C66 -! SHOULD NOT BE PART OF LATTICE END interface lattice_forestProjection_edge module procedure slipProjection_transverse @@ -384,7 +378,8 @@ module lattice lattice_init, & lattice_equivalent_nu, & lattice_equivalent_mu, & - lattice_applyLatticeSymmetry33, & + lattice_symmetrize_33, & + lattice_symmetrize_C66, & lattice_SchmidMatrix_slip, & lattice_SchmidMatrix_twin, & lattice_SchmidMatrix_trans, & @@ -414,53 +409,9 @@ contains !-------------------------------------------------------------------------------------------------- subroutine lattice_init - integer :: Nphases, ph,i - class(tNode), pointer :: & - phases, & - phase, & - mech, & - elasticity - print'(/,a)', ' <<<+- lattice init -+>>>'; flush(IO_STDOUT) -! SHOULD NOT BE PART OF LATTICE BEGIN - - phases => config_material%get('phase') - Nphases = phases%length - - allocate(lattice_C66(6,6,Nphases), source=0.0_pReal) - - allocate(lattice_mu, lattice_nu,& - source=[(0.0_pReal,i=1,Nphases)]) - - do ph = 1, phases%length - phase => phases%get(ph) - mech => phase%get('mechanical') - elasticity => mech%get('elastic') - lattice_C66(1,1,ph) = elasticity%get_asFloat('C_11') - lattice_C66(1,2,ph) = elasticity%get_asFloat('C_12') - lattice_C66(4,4,ph) = elasticity%get_asFloat('C_44') - - lattice_C66(1,3,ph) = elasticity%get_asFloat('C_13',defaultVal=0.0_pReal) - lattice_C66(2,3,ph) = elasticity%get_asFloat('C_23',defaultVal=0.0_pReal) - lattice_C66(3,3,ph) = elasticity%get_asFloat('C_33',defaultVal=0.0_pReal) - lattice_C66(6,6,ph) = elasticity%get_asFloat('C_66',defaultVal=0.0_pReal) - - lattice_C66(1:6,1:6,ph) = applyLatticeSymmetryC66(lattice_C66(1:6,1:6,ph),phase%get_asString('lattice')) - - lattice_nu(ph) = lattice_equivalent_nu(lattice_C66(1:6,1:6,ph),'voigt') - lattice_mu(ph) = lattice_equivalent_mu(lattice_C66(1:6,1:6,ph),'voigt') - - lattice_C66(1:6,1:6,ph) = math_sym3333to66(math_Voigt66to3333(lattice_C66(1:6,1:6,ph))) ! Literature data is in Voigt notation - do i = 1, 6 - if (abs(lattice_C66(i,i,ph)) @brief Return 3x3 tensor with symmetry according to given crystal structure !-------------------------------------------------------------------------------------------------- -pure function lattice_applyLatticeSymmetry33(T,structure) result(T_sym) +pure function lattice_symmetrize_33(T,lattice) result(T_sym) real(pReal), dimension(3,3) :: T_sym real(pReal), dimension(3,3), intent(in) :: T - character(len=*), intent(in) :: structure + character(len=2), intent(in) :: lattice T_sym = 0.0_pReal - select case(structure) + select case(lattice) case('cF','cI') T_sym(1,1) = T(1,1) T_sym(2,2) = T(1,1) @@ -1649,26 +1600,26 @@ pure function lattice_applyLatticeSymmetry33(T,structure) result(T_sym) T_sym(3,3) = T(3,3) end select -end function lattice_applyLatticeSymmetry33 +end function lattice_symmetrize_33 !-------------------------------------------------------------------------------------------------- !> @brief Return stiffness matrix in 6x6 notation with symmetry according to given crystal structure !> @details J. A. Rayne and B. S. Chandrasekhar Phys. Rev. 120, 1658 Erratum Phys. Rev. 122, 1962 !-------------------------------------------------------------------------------------------------- -pure function applyLatticeSymmetryC66(C66,structure) result(C66_sym) +pure function lattice_symmetrize_C66(C66,lattice) result(C66_sym) real(pReal), dimension(6,6) :: C66_sym real(pReal), dimension(6,6), intent(in) :: C66 - character(len=*), intent(in) :: structure + character(len=*), intent(in) :: lattice integer :: i,j C66_sym = 0.0_pReal - select case(structure) + select case(lattice) case ('cF','cI') C66_sym(1,1) = C66(1,1); C66_sym(2,2) = C66(1,1); C66_sym(3,3) = C66(1,1) C66_sym(1,2) = C66(1,2); C66_sym(1,3) = C66(1,2); C66_sym(2,3) = C66(1,2) @@ -1695,7 +1646,7 @@ pure function applyLatticeSymmetryC66(C66,structure) result(C66_sym) enddo enddo -end function applyLatticeSymmetryC66 +end function lattice_symmetrize_C66 !-------------------------------------------------------------------------------------------------- @@ -2203,10 +2154,10 @@ subroutine selfTest do i = 1, 10 call random_number(C) - C_cF = applyLatticeSymmetryC66(C,'cI') - C_cI = applyLatticeSymmetryC66(C,'cF') - C_hP = applyLatticeSymmetryC66(C,'hP') - C_tI = applyLatticeSymmetryC66(C,'tI') + C_cF = lattice_symmetrize_C66(C,'cI') + C_cI = lattice_symmetrize_C66(C,'cF') + C_hP = lattice_symmetrize_C66(C,'hP') + C_tI = lattice_symmetrize_C66(C,'tI') if (any(dNeq(C_cI,transpose(C_cF)))) error stop 'SymmetryC66/cI-cF' if (any(dNeq(C_cF,transpose(C_cI)))) error stop 'SymmetryC66/cF-cI' @@ -2226,10 +2177,10 @@ subroutine selfTest if (any(dNeq(C(4,4),[C_tI(4,4),C_tI(5,5)]))) error stop 'SymmetryC_44-55/tI' call random_number(T) - T_cF = lattice_applyLatticeSymmetry33(T,'cI') - T_cI = lattice_applyLatticeSymmetry33(T,'cF') - T_hP = lattice_applyLatticeSymmetry33(T,'hP') - T_tI = lattice_applyLatticeSymmetry33(T,'tI') + T_cF = lattice_symmetrize_33(T,'cI') + T_cI = lattice_symmetrize_33(T,'cF') + T_hP = lattice_symmetrize_33(T,'hP') + T_tI = lattice_symmetrize_33(T,'tI') if (any(dNeq0(T_cF) .and. math_I3<1.0_pReal)) error stop 'Symmetry33/c' if (any(dNeq0(T_hP) .and. math_I3<1.0_pReal)) error stop 'Symmetry33/hP' @@ -2244,7 +2195,7 @@ subroutine selfTest call random_number(C) C(1,1) = C(1,1) + C(1,2) + 0.1_pReal C(4,4) = 0.5_pReal * (C(1,1) - C(1,2)) - C = applyLatticeSymmetryC66(C,'cI') + C = lattice_symmetrize_C66(C,'cI') if(dNeq(C(4,4),lattice_equivalent_mu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_mu/voigt' if(dNeq(C(4,4),lattice_equivalent_mu(C,'reuss'),1.0e-12_pReal)) error stop 'equivalent_mu/reuss' diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index e749a7fa0..b69969898 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -107,7 +107,7 @@ module subroutine damage_init param(ph)%mu = source%get_asFloat('mu') param(ph)%D(1,1) = source%get_asFloat('D_11') if (any(phase_lattice(ph) == ['hP','tI'])) param(ph)%D(3,3) = source%get_asFloat('D_33') - param(ph)%D = lattice_applyLatticeSymmetry33(param(ph)%D,phase_lattice(ph)) + param(ph)%D = lattice_symmetrize_33(param(ph)%D,phase_lattice(ph)) endif enddo diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 8b04a0009..c972ed3cd 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -168,6 +168,20 @@ submodule(phase) mechanical integer, intent(in) :: ph,en end function plastic_dislotwin_homogenizedC + module function elastic_C66(ph) result(C66) + real(pReal), dimension(6,6) :: C66 + integer, intent(in) :: ph + end function elastic_C66 + + module function elastic_mu(ph) result(mu) + real(pReal) :: mu + integer, intent(in) :: ph + end function elastic_mu + + module function elastic_nu(ph) result(nu) + real(pReal) :: nu + integer, intent(in) :: ph + end function elastic_nu end interface type :: tOutput !< new requested output (per phase) diff --git a/src/phase_mechanical_eigen_thermalexpansion.f90 b/src/phase_mechanical_eigen_thermalexpansion.f90 index dba02d70e..0976c1dfa 100644 --- a/src/phase_mechanical_eigen_thermalexpansion.f90 +++ b/src/phase_mechanical_eigen_thermalexpansion.f90 @@ -69,7 +69,7 @@ module function thermalexpansion_init(kinematics_length) result(myKinematics) prm%A(3,3,3) = kinematic_type%get_asFloat('A_33,T^2',defaultVal=0.0_pReal) endif do i=1, size(prm%A,3) - prm%A(1:3,1:3,i) = lattice_applyLatticeSymmetry33(prm%A(1:3,1:3,i),phase_lattice(p)) + prm%A(1:3,1:3,i) = lattice_symmetrize_33(prm%A(1:3,1:3,i),phase_lattice(p)) enddo end associate endif diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90 index 3249b2921..154ad4640 100644 --- a/src/phase_mechanical_elastic.f90 +++ b/src/phase_mechanical_elastic.f90 @@ -2,7 +2,10 @@ submodule(phase:mechanical) elastic type :: tParameters real(pReal), dimension(6,6) :: & - C66 !< Elastic constants in Voig notation + C66 !< Elastic constants in Voigt notation + real(pReal) :: & + mu, & + nu end type tParameters type(tParameters), allocatable, dimension(:) :: param @@ -37,6 +40,7 @@ module subroutine elastic_init(phases) if (elastic%get_asString('type') /= 'Hooke') call IO_error(200,ext_msg=elastic%get_asString('type')) associate(prm => param(ph)) + prm%C66(1,1) = elastic%get_asFloat('C_11') prm%C66(1,2) = elastic%get_asFloat('C_12') prm%C66(4,4) = elastic%get_asFloat('C_44') @@ -46,6 +50,12 @@ module subroutine elastic_init(phases) prm%C66(3,3) = elastic%get_asFloat('C_33') endif if (phase_lattice(ph) == 'tI') prm%C66(6,6) = elastic%get_asFloat('C_66') + + prm%C66 = lattice_symmetrize_C66(prm%C66,phase_lattice(ph)) + + prm%nu = lattice_equivalent_nu(prm%C66,'voigt') + prm%mu = lattice_equivalent_mu(prm%C66,'voigt') + end associate enddo @@ -104,10 +114,38 @@ module function phase_homogenizedC(ph,en) result(C) case (PLASTICITY_DISLOTWIN_ID) plasticType C = plastic_dislotwin_homogenizedC(ph,en) case default plasticType - C = lattice_C66(1:6,1:6,ph) + C = param(ph)%C66 end select plasticType end function phase_homogenizedC +module function elastic_C66(ph) result(C66) + real(pReal), dimension(6,6) :: C66 + integer, intent(in) :: ph + + + C66 = param(ph)%C66 + +end function elastic_C66 + +module function elastic_mu(ph) result(mu) + + real(pReal) :: mu + integer, intent(in) :: ph + + + mu = param(ph)%mu + +end function elastic_mu + +module function elastic_nu(ph) result(nu) + + real(pReal) :: nu + integer, intent(in) :: ph + + + nu = param(ph)%mu + +end function elastic_nu end submodule elastic diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index 8b50c6845..bef7b5f6a 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -129,8 +129,7 @@ module function plastic_dislotungsten_init() result(myPlasticity) prm%output = pl%get_as1dString('output',defaultVal=emptyStringArray) #endif - ! This data is read in already in lattice - prm%mu = lattice_mu(ph) + prm%mu = elastic_mu(ph) !-------------------------------------------------------------------------------------------------- ! slip related parameters diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index 36b8cfe32..33f148036 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -184,9 +184,9 @@ module function plastic_dislotwin_init() result(myPlasticity) #endif ! This data is read in already in lattice - prm%mu = lattice_mu(ph) - prm%nu = lattice_nu(ph) - prm%C66 = lattice_C66(1:6,1:6,ph) + prm%mu = elastic_mu(ph) + prm%nu = elastic_nu(ph) + prm%C66 = elastic_C66(ph) !-------------------------------------------------------------------------------------------------- ! slip related parameters @@ -481,22 +481,21 @@ module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC) real(pReal) :: f_unrotated - associate(prm => param(ph),& - stt => state(ph)) + associate(prm => param(ph), stt => state(ph)) - f_unrotated = 1.0_pReal & - - sum(stt%f_tw(1:prm%sum_N_tw,en)) & - - sum(stt%f_tr(1:prm%sum_N_tr,en)) + f_unrotated = 1.0_pReal & + - sum(stt%f_tw(1:prm%sum_N_tw,en)) & + - sum(stt%f_tr(1:prm%sum_N_tr,en)) - homogenizedC = f_unrotated * prm%C66 - do i=1,prm%sum_N_tw - homogenizedC = homogenizedC & - + stt%f_tw(i,en)*prm%C66_tw(1:6,1:6,i) - enddo - do i=1,prm%sum_N_tr - homogenizedC = homogenizedC & - + stt%f_tr(i,en)*prm%C66_tr(1:6,1:6,i) - enddo + homogenizedC = f_unrotated * prm%C66 + do i=1,prm%sum_N_tw + homogenizedC = homogenizedC & + + stt%f_tw(i,en)*prm%C66_tw(1:6,1:6,i) + enddo + do i=1,prm%sum_N_tr + homogenizedC = homogenizedC & + + stt%f_tr(i,en)*prm%C66_tr(1:6,1:6,i) + enddo end associate diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index 4ff60c67e..101c0d299 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -240,9 +240,8 @@ module function plastic_nonlocal_init() result(myPlasticity) prm%atol_rho = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal) - ! This data is read in already in lattice - prm%mu = lattice_mu(ph) - prm%nu = lattice_nu(ph) + prm%mu = elastic_mu(ph) + prm%nu = elastic_nu(ph) ini%N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray) prm%sum_N_sl = sum(abs(ini%N_sl)) diff --git a/src/phase_thermal.f90 b/src/phase_thermal.f90 index 895246dfa..9a75263ca 100644 --- a/src/phase_thermal.f90 +++ b/src/phase_thermal.f90 @@ -103,7 +103,7 @@ module subroutine thermal_init(phases) param(ph)%C_p = thermal%get_asFloat('C_p',defaultVal=0.0_pReal) ! ToDo: make mandatory? param(ph)%K(1,1) = thermal%get_asFloat('K_11',defaultVal=0.0_pReal) ! ToDo: make mandatory? param(ph)%K(3,3) = thermal%get_asFloat('K_33',defaultVal=0.0_pReal) ! ToDo: depends on symmtery - param(ph)%K = lattice_applyLatticeSymmetry33(param(ph)%K,phase_lattice(ph)) + param(ph)%K = lattice_symmetrize_33(param(ph)%K,phase_lattice(ph)) sources => thermal%get('source',defaultVal=emptyList) thermal_Nsources(ph) = sources%length