From 299c47fd6f65850f74cef8fbe54179f0d06f7c8e Mon Sep 17 00:00:00 2001 From: Martin Diehl Date: Mon, 24 May 2021 20:33:50 +0200 Subject: [PATCH] prepare for varying C66 - check structure centrally - pure function with guaranteed return/no stop --- src/lattice.f90 | 15 +++++---------- src/phase_mechanical_elastic.f90 | 15 +++++++++++---- 2 files changed, 16 insertions(+), 14 deletions(-) diff --git a/src/lattice.f90 b/src/lattice.f90 index a285d4ffb..405431189 100644 --- a/src/lattice.f90 +++ b/src/lattice.f90 @@ -1664,22 +1664,21 @@ end function lattice_labels_slip !-------------------------------------------------------------------------------------------------- !> @brief Return 3x3 tensor with symmetry according to given crystal structure !-------------------------------------------------------------------------------------------------- -function lattice_applyLatticeSymmetry33(T,structure) result(T_sym) +pure function lattice_applyLatticeSymmetry33(T,structure) result(T_sym) real(pReal), dimension(3,3) :: T_sym real(pReal), dimension(3,3), intent(in) :: T character(len=*), intent(in) :: structure - integer :: k T_sym = 0.0_pReal select case(structure) case('cF','cI') - do k=1,3 - T_sym(k,k) = T(1,1) - enddo + T_sym(1,1) = T(1,1) + T_sym(2,2) = T(1,1) + T_sym(3,3) = T(1,1) case('hP') ! MD TODO: I think that 'tI' has the same symmetry as 'hP' for 2nd order tensors T_sym(1,1) = T(1,1) T_sym(2,2) = T(1,1) @@ -1688,8 +1687,6 @@ function lattice_applyLatticeSymmetry33(T,structure) result(T_sym) T_sym(1,1) = T(1,1) T_sym(2,2) = T(2,2) T_sym(3,3) = T(3,3) - case default - call IO_error(137,ext_msg='lattice_applyLatticeSymmetry33: '//trim(structure)) end select end function lattice_applyLatticeSymmetry33 @@ -1699,7 +1696,7 @@ end function lattice_applyLatticeSymmetry33 !> @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 !-------------------------------------------------------------------------------------------------- -function applyLatticeSymmetryC66(C66,structure) result(C66_sym) +pure function applyLatticeSymmetryC66(C66,structure) result(C66_sym) real(pReal), dimension(6,6) :: C66_sym @@ -1730,8 +1727,6 @@ function applyLatticeSymmetryC66(C66,structure) result(C66_sym) C66_sym(1,3) = C66(1,3); C66_sym(2,3) = C66(1,3) C66_sym(4,4) = C66(4,4); C66_sym(5,5) = C66(4,4) C66_sym(6,6) = C66(6,6) - case default - call IO_error(137,ext_msg='applyLatticeSymmetryC66: '//trim(structure)) end select do i = 1, 6 diff --git a/src/phase_mechanical_elastic.f90 b/src/phase_mechanical_elastic.f90 index a3a204557..3397cc572 100644 --- a/src/phase_mechanical_elastic.f90 +++ b/src/phase_mechanical_elastic.f90 @@ -21,6 +21,7 @@ module subroutine elastic_init(phases) phase, & mech, & elastic + character(len=:), allocatable :: struct print'(/,a)', ' <<<+- phase:mechanical:elastic init -+>>>' @@ -34,14 +35,20 @@ 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)) + struct = phase%get_asString('lattice') + if (struct /= 'cI' .and. struct /= 'cF' .and. struct /= 'hP' .and. struct /= 'tI') & + call IO_error(137,ext_msg=trim(struct)) + 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') - prm%C66(1,3) = elastic%get_asFloat('C_13',defaultVal=0.0_pReal) - prm%C66(2,3) = elastic%get_asFloat('C_23',defaultVal=0.0_pReal) - prm%C66(3,3) = elastic%get_asFloat('C_33',defaultVal=0.0_pReal) - prm%C66(6,6) = elastic%get_asFloat('C_66',defaultVal=0.0_pReal) + if (struct == 'hP' .or. struct == 'tI') then + prm%C66(1,3) = elastic%get_asFloat('C_13') + prm%C66(2,3) = elastic%get_asFloat('C_23') + prm%C66(3,3) = elastic%get_asFloat('C_33') + endif + if (struct == 'tI') prm%C66(6,6) = elastic%get_asFloat('C_66') end associate enddo