prepare for varying C66

- check structure centrally
- pure function with guaranteed return/no stop

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(1,1) = T(1,1)
-        T_sym(k,k) = T(1,1)
+      T_sym(2,2) = T(1,1)
-      enddo
+      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

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)
+      if (struct == 'hP' .or. struct == 'tI') then
-      prm%C66(2,3) = elastic%get_asFloat('C_23',defaultVal=0.0_pReal)
+        prm%C66(1,3) = elastic%get_asFloat('C_13')
-      prm%C66(3,3) = elastic%get_asFloat('C_33',defaultVal=0.0_pReal)
+        prm%C66(2,3) = elastic%get_asFloat('C_23')
-      prm%C66(6,6) = elastic%get_asFloat('C_66',defaultVal=0.0_pReal)
+        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