Merge remote-tracking branch 'origin/development' into restructure-Orientation-2

PRIVATE

@@ -1 +1 @@
-Subproject commit 8fec909d1931b092b223b0560dd30c3339c6e5a7
+Subproject commit 174ecac2d3ab7596bdb60184d6bb9e1a52cb7378

examples/config/numerics.yaml

@@ -79,6 +79,5 @@ commercialFEM:
   unitlength:                  1            # physical length of one computational length unit
 
 generic:
-  charLength:                  1.0          # characteristic length scale for gradient problems.
   random_seed:                 0            # fixed seeding for pseudo-random number generator, Default 0: use random seed.
   residualStiffness:           1.0e-6       # non-zero residual damage.

examples/config/phase/damage/anisobrittle_cubic.yaml

@@ -7,5 +7,5 @@ q: 20
 
 output: [f_phi]
 
-K_11: 1.0
+D_11: 1.0
 mu: 0.001

examples/config/phase/damage/isobrittle_generic.yaml

@@ -1,9 +1,7 @@
 type: isobrittle
 W_crit: 1400000.0
-m: 1.0
-isoBrittle_atol: 0.01
 
 output: [f_phi]
 
-K_11: 1.0
+D_11: 1.0
 mu: 0.001

python/damask/VERSION

@@ -1 +1 @@
-v3.0.0-alpha4-94-g63fee141b
+v3.0.0-alpha4-137-gb69b85754

src/CPFEM.f90

@@ -189,9 +189,11 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
       CPFEM_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_eye(6)
 
     else validCalculation
-      if (debugCPFEM%extensive) &
+      if (debugCPFEM%extensive)  print'(a,i8,1x,i2)', '<< CPFEM >> calculation for elFE ip ',elFE,ip
-        print'(a,i8,1x,i2)', '<< CPFEM >> calculation for elFE ip ',elFE,ip
+      call homogenization_mechanical_response(dt,[ip,ip],[elCP,elCP])
-      call materialpoint_stressAndItsTangent(dt,[ip,ip],[elCP,elCP])
+      if (.not. terminallyIll) &
+        call homogenization_mechanical_response2(dt,[ip,ip],[elCP,elCP])
+
 
       terminalIllness: if (terminallyIll) then
 

src/IO.f90

@@ -81,12 +81,7 @@ function IO_readlines(fileName) result(fileContent)
 
   rawData = IO_read(fileName)
 
-!--------------------------------------------------------------------------------------------------
+  N_lines = count([(rawData(l:l) == IO_EOL,l=1,len(rawData))])
-! count lines to allocate string array
-  N_lines = 0
-  do l=1, len(rawData)
-    if (rawData(l:l) == IO_EOL) N_lines = N_lines+1
-  enddo
   allocate(fileContent(N_lines))
 
 !--------------------------------------------------------------------------------------------------
@@ -261,7 +256,7 @@ pure function IO_lc(string)
 
   do i=1,len(string)
     n = index(UPPER,string(i:i))
-    if(n/=0) then
+    if (n/=0) then
       IO_lc(i:i) = LOWER(n:n)
     else
       IO_lc(i:i) = string(i:i)

src/grid/DAMASK_grid.f90

@@ -482,20 +482,13 @@ subroutine getMaskedTensor(values,mask,tensor)
 
 
   values = 0.0
-  if (tensor%length == 9) then ! temporary support for deprecated 1D tensor
+  do i = 1,3
-    do i = 1,9
+    row => tensor%get(i)
-      mask((i-1)/3+1,mod(i-1,3)+1) = tensor%get_asString(i) /= 'x'
+    do j = 1,3
-      if (mask((i-1)/3+1,mod(i-1,3)+1)) values((i-1)/3+1,mod(i-1,3)+1) = tensor%get_asFloat(i)
+      mask(i,j) = row%get_asString(j) /= 'x'
+      if (mask(i,j)) values(i,j) = row%get_asFloat(j)
     enddo
-  else
+  enddo
-    do i = 1,3
-      row => tensor%get(i)
-      do j = 1,3
-        mask(i,j) = row%get_asString(j) /= 'x'                      ! ToDo change to np.masked behavior
-        if (mask(i,j)) values(i,j) = row%get_asFloat(j)
-      enddo
-    enddo
-  endif
 
 end subroutine
 

src/grid/grid_mech_FEM.f90

@@ -261,7 +261,7 @@ subroutine grid_mechanical_FEM_init
     F         = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3)
   endif restartRead
 
-  homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3])                            ! set starting condition for materialpoint_stressAndItsTangent
+  homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3])                            ! set starting condition for homogenization_mechanical_response
   call utilities_updateCoords(F)
   call utilities_constitutiveResponse(P_current,P_av,C_volAvg,devNull, &                            ! stress field, stress avg, global average of stiffness and (min+max)/2
                                       F, &                                                          ! target F
@@ -490,8 +490,8 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,i
     BCTol
 
   err_div = fnorm*sqrt(wgt)*geomSize(1)/scaledGeomSize(1)/detJ
-  divTol = max(maxval(abs(P_av))*num%eps_div_rtol   ,num%eps_div_atol)
+  divTol = max(maxval(abs(P_av))*num%eps_div_rtol, num%eps_div_atol)
-  BCTol  = max(maxval(abs(P_av))*num%eps_stress_rtol,num%eps_stress_atol)
+  BCTol  = max(maxval(abs(P_av))*num%eps_stress_rtol, num%eps_stress_atol)
 
   if ((totalIter >= num%itmin .and. all([err_div/divTol, err_BC/BCTol] < 1.0_pReal)) &
        .or. terminallyIll) then
@@ -502,8 +502,6 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,fnorm,reason,dummy,i
     reason = 0
   endif
 
-!--------------------------------------------------------------------------------------------------
-! report
   print'(1/,a)', ' ... reporting .............................................................'
   print'(1/,a,f12.2,a,es8.2,a,es9.2,a)', ' error divergence = ', &
           err_div/divTol,  ' (',err_div,' / m, tol = ',divTol,')'

src/grid/grid_mech_spectral_basic.f90

@@ -211,7 +211,7 @@ subroutine grid_mechanical_spectral_basic_init
     F = reshape(F_lastInc,[9,grid(1),grid(2),grid3])
   endif restartRead
 
-  homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3])                            ! set starting condition for materialpoint_stressAndItsTangent
+  homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3])                            ! set starting condition for homogenization_mechanical_response
   call utilities_updateCoords(reshape(F,shape(F_lastInc)))
   call utilities_constitutiveResponse(P,P_av,C_volAvg,C_minMaxAvg, &                                ! stress field, stress avg, global average of stiffness and (min+max)/2
                                       reshape(F,shape(F_lastInc)), &                                ! target F
@@ -429,8 +429,8 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm
     divTol, &
     BCTol
 
-  divTol = max(maxval(abs(P_av))*num%eps_div_rtol   ,num%eps_div_atol)
+  divTol = max(maxval(abs(P_av))*num%eps_div_rtol, num%eps_div_atol)
-  BCTol  = max(maxval(abs(P_av))*num%eps_stress_rtol,num%eps_stress_atol)
+  BCTol = max(maxval(abs(P_av))*num%eps_stress_rtol, num%eps_stress_atol)
 
   if ((totalIter >= num%itmin .and. all([err_div/divTol, err_BC/BCTol] < 1.0_pReal)) &
        .or. terminallyIll) then
@@ -441,8 +441,6 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm
     reason = 0
   endif
 
-!--------------------------------------------------------------------------------------------------
-! report
   print'(1/,a)', ' ... reporting .............................................................'
   print'(1/,a,f12.2,a,es8.2,a,es9.2,a)', ' error divergence = ', &
           err_div/divTol,  ' (',err_div,' / m, tol = ',divTol,')'

src/grid/grid_mech_spectral_polarisation.f90

@@ -237,7 +237,7 @@ subroutine grid_mechanical_spectral_polarisation_init
     F_tau_lastInc = 2.0_pReal*F_lastInc
   endif restartRead
 
-  homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3])                            ! set starting condition for materialpoint_stressAndItsTangent
+  homogenization_F0 = reshape(F_lastInc, [3,3,product(grid(1:2))*grid3])                            ! set starting condition for homogenization_mechanical_response
   call utilities_updateCoords(reshape(F,shape(F_lastInc)))
   call utilities_constitutiveResponse(P,P_av,C_volAvg,C_minMaxAvg, &                                ! stress field, stress avg, global average of stiffness and (min+max)/2
                                       reshape(F,shape(F_lastInc)), &                                ! target F
@@ -487,9 +487,9 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm
     divTol, &
     BCTol
 
-  curlTol    = max(maxval(abs(F_aim-math_I3))*num%eps_curl_rtol  ,num%eps_curl_atol)
+  curlTol = max(maxval(abs(F_aim-math_I3))*num%eps_curl_rtol, num%eps_curl_atol)
-  divTol     = max(maxval(abs(P_av))         *num%eps_div_rtol   ,num%eps_div_atol)
+  divTol = max(maxval(abs(P_av))*num%eps_div_rtol, num%eps_div_atol)
-  BCTol      = max(maxval(abs(P_av))         *num%eps_stress_rtol,num%eps_stress_atol)
+  BCTol = max(maxval(abs(P_av))*num%eps_stress_rtol, num%eps_stress_atol)
 
   if ((totalIter >= num%itmin .and. all([err_div/divTol, err_curl/curlTol, err_BC/BCTol] < 1.0_pReal)) &
        .or. terminallyIll) then
@@ -500,8 +500,6 @@ subroutine converged(snes_local,PETScIter,devNull1,devNull2,devNull3,reason,dumm
     reason = 0
   endif
 
-!--------------------------------------------------------------------------------------------------
-! report
   print'(1/,a)', ' ... reporting .............................................................'
   print'(1/,a,f12.2,a,es8.2,a,es9.2,a)', ' error divergence = ', &
             err_div/divTol,  ' (',err_div, ' / m, tol = ',divTol,')'
@@ -542,13 +540,13 @@ subroutine formResidual(in, FandF_tau, &
 !---------------------------------------------------------------------------------------------------
 
   F              => FandF_tau(1:3,1:3,1,&
-                                 XG_RANGE,YG_RANGE,ZG_RANGE)
+                              XG_RANGE,YG_RANGE,ZG_RANGE)
   F_tau          => FandF_tau(1:3,1:3,2,&
-                                 XG_RANGE,YG_RANGE,ZG_RANGE)
+                              XG_RANGE,YG_RANGE,ZG_RANGE)
   residual_F     => residuum(1:3,1:3,1,&
-                                 X_RANGE, Y_RANGE, Z_RANGE)
+                             X_RANGE, Y_RANGE, Z_RANGE)
   residual_F_tau => residuum(1:3,1:3,2,&
-                                 X_RANGE, Y_RANGE, Z_RANGE)
+                             X_RANGE, Y_RANGE, Z_RANGE)
 
   F_av = sum(sum(sum(F,dim=5),dim=4),dim=3) * wgt
   call MPI_Allreduce(MPI_IN_PLACE,F_av,9,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr)

src/grid/spectral_utilities.f90

@@ -815,10 +815,14 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
 
   homogenization_F  = reshape(F,[3,3,product(grid(1:2))*grid3])                                     ! set materialpoint target F to estimated field
 
-  call materialpoint_stressAndItsTangent(timeinc,[1,1],[1,product(grid(1:2))*grid3])                ! calculate P field
+  call homogenization_mechanical_response(timeinc,[1,1],[1,product(grid(1:2))*grid3])               ! calculate P field
+  if (.not. terminallyIll) &
+    call homogenization_thermal_response(timeinc,[1,1],[1,product(grid(1:2))*grid3])
+  if (.not. terminallyIll) &
+    call homogenization_mechanical_response2(timeinc,[1,1],[1,product(grid(1:2))*grid3])
 
   P = reshape(homogenization_P, [3,3,grid(1),grid(2),grid3])
-  P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt                                                   ! average of P
+  P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt
   call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr)
   if (debugRotation) print'(/,a,/,2(3(2x,f12.4,1x)/),3(2x,f12.4,1x))', &
     ' Piola--Kirchhoff stress (lab) / MPa =', transpose(P_av)*1.e-6_pReal

src/homogenization.f90

@@ -41,15 +41,6 @@ module homogenization
   integer(kind(DAMAGE_none_ID)),              dimension(:),   allocatable :: &
     damage_type                                                                                     !< nonlocal damage model
 
-  type, private :: tNumerics_damage
-    real(pReal) :: &
-    charLength                                                                                      !< characteristic length scale for gradient problems
-  end type tNumerics_damage
-
-  type(tNumerics_damage), private :: &
-    num_damage
-
-
   logical, public :: &
     terminallyIll = .false.                                                                         !< at least one material point is terminally ill
 
@@ -101,8 +92,8 @@ module homogenization
       integer,     intent(in) :: ce
     end subroutine damage_partition
 
-    module subroutine mechanical_homogenize(dt,ce)
+    module subroutine mechanical_homogenize(Delta_t,ce)
-     real(pReal), intent(in) :: dt
+     real(pReal), intent(in) :: Delta_t
      integer, intent(in) :: &
        ce                                                                                           !< cell
     end subroutine mechanical_homogenize
@@ -178,7 +169,9 @@ module homogenization
 
   public ::  &
     homogenization_init, &
-    materialpoint_stressAndItsTangent, &
+    homogenization_mechanical_response, &
+    homogenization_mechanical_response2, &
+    homogenization_thermal_response, &
     homogenization_mu_T, &
     homogenization_K_T, &
     homogenization_f_T, &
@@ -227,24 +220,24 @@ end subroutine homogenization_init
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief  parallelized calculation of stress and corresponding tangent at material points
+!> @brief
 !--------------------------------------------------------------------------------------------------
-subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execElem)
+subroutine homogenization_mechanical_response(Delta_t,FEsolving_execIP,FEsolving_execElem)
 
-  real(pReal), intent(in) :: dt                                                                     !< time increment
+  real(pReal), intent(in) :: Delta_t                                                                !< time increment
   integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP
   integer :: &
     NiterationMPstate, &
     ip, &                                                                                            !< integration point number
     el, &                                                                                            !< element number
-    co, ce, ho, en, ph
+    co, ce, ho, en
   logical :: &
     converged
   logical, dimension(2) :: &
     doneAndHappy
 
-  !$OMP PARALLEL
+
-  !$OMP DO PRIVATE(ce,en,ho,NiterationMPstate,converged,doneAndHappy)
+  !$OMP PARALLEL DO PRIVATE(ce,en,ho,NiterationMPstate,converged,doneAndHappy)
   do el = FEsolving_execElem(1),FEsolving_execElem(2)
 
     do ip = FEsolving_execIP(1),FEsolving_execIP(2)
@@ -266,65 +259,89 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
         NiterationMPstate = NiterationMPstate + 1
 
         call mechanical_partition(homogenization_F(1:3,1:3,ce),ce)
-        converged = .true.
+        converged = all([(phase_mechanical_constitutive(Delta_t,co,ip,el),co=1,homogenization_Nconstituents(ho))])
-        do co = 1, homogenization_Nconstituents(ho)
-          converged = converged .and. crystallite_stress(dt,co,ip,el)
-        enddo
-
         if (converged) then
-          doneAndHappy = mechanical_updateState(dt,homogenization_F(1:3,1:3,ce),ce)
+          doneAndHappy = mechanical_updateState(Delta_t,homogenization_F(1:3,1:3,ce),ce)
           converged = all(doneAndHappy)
         else
           doneAndHappy = [.true.,.false.]
         endif
-
       enddo convergenceLooping
+
+      converged = converged .and. all([(phase_damage_constitutive(Delta_t,co,ip,el),co=1,homogenization_Nconstituents(ho))])
+
       if (.not. converged) then
-        if (.not. terminallyIll) print*, ' Integration point ', ip,' at element ', el, ' terminally ill'
+        if (.not. terminallyIll) print*, ' Cell ', ce, ' terminally ill'
         terminallyIll = .true.
       endif
     enddo
   enddo
-  !$OMP END DO
+  !$OMP END PARALLEL DO
 
-  if (.not. terminallyIll) then
+end subroutine homogenization_mechanical_response
-    !$OMP DO PRIVATE(ho,ph,ce)
+
-    do el = FEsolving_execElem(1),FEsolving_execElem(2)
+
-      if (terminallyIll) continue
+!--------------------------------------------------------------------------------------------------
-      do ip = FEsolving_execIP(1),FEsolving_execIP(2)
+!> @brief
-        ce = (el-1)*discretization_nIPs + ip
+!--------------------------------------------------------------------------------------------------
-        ho = material_homogenizationID(ce)
+subroutine homogenization_thermal_response(Delta_t,FEsolving_execIP,FEsolving_execElem)
-        call thermal_partition(ce)
+
-        do co = 1, homogenization_Nconstituents(ho)
+  real(pReal), intent(in) :: Delta_t                                                                !< time increment
-          ph = material_phaseID(co,ce)
+  integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP
-          if (.not. thermal_stress(dt,ph,material_phaseMemberAt(co,ip,el))) then
+  integer :: &
-            if (.not. terminallyIll) &                                                              ! so first signals terminally ill...
+    ip, &                                                                                           !< integration point number
-              print*, ' Integration point ', ip,' at element ', el, ' terminally ill'
+    el, &                                                                                           !< element number
-            terminallyIll = .true.                                                                  ! ...and kills all others
+    co, ce, ho
-          endif
+
-        enddo
+
+  !$OMP PARALLEL DO PRIVATE(ho,ce)
+  do el = FEsolving_execElem(1),FEsolving_execElem(2)
+    if (terminallyIll) continue
+    do ip = FEsolving_execIP(1),FEsolving_execIP(2)
+      ce = (el-1)*discretization_nIPs + ip
+      ho = material_homogenizationID(ce)
+      call thermal_partition(ce)
+      do co = 1, homogenization_Nconstituents(ho)
+        if (.not. phase_thermal_constitutive(Delta_t,material_phaseID(co,ce),material_phaseEntry(co,ce))) then
+          if (.not. terminallyIll) print*, ' Cell ', ce, ' terminally ill'
+          terminallyIll = .true.
+        endif
       enddo
     enddo
-    !$OMP END DO
+  enddo
+  !$OMP END PARALLEL DO
 
-    !$OMP DO PRIVATE(ho,ce)
+end subroutine homogenization_thermal_response
-    elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2)
-      IpLooping3: do ip = FEsolving_execIP(1),FEsolving_execIP(2)
-        ce = (el-1)*discretization_nIPs + ip
-        ho = material_homogenizationID(ce)
-        do co = 1, homogenization_Nconstituents(ho)
-          call crystallite_orientations(co,ip,el)
-        enddo
-        call mechanical_homogenize(dt,ce)
-      enddo IpLooping3
-    enddo elementLooping3
-    !$OMP END DO
-  else
-    print'(/,a,/)', ' << HOMOG >> Material Point terminally ill'
-  endif
-  !$OMP END PARALLEL
 
-end subroutine materialpoint_stressAndItsTangent
+
+!--------------------------------------------------------------------------------------------------
+!> @brief
+!--------------------------------------------------------------------------------------------------
+subroutine homogenization_mechanical_response2(Delta_t,FEsolving_execIP,FEsolving_execElem)
+
+  real(pReal), intent(in) :: Delta_t                                                                !< time increment
+  integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP
+  integer :: &
+    ip, &                                                                                           !< integration point number
+    el, &                                                                                           !< element number
+    co, ce, ho
+
+
+  !$OMP PARALLEL DO PRIVATE(ho,ce)
+  elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2)
+    IpLooping3: do ip = FEsolving_execIP(1),FEsolving_execIP(2)
+      ce = (el-1)*discretization_nIPs + ip
+      ho = material_homogenizationID(ce)
+      do co = 1, homogenization_Nconstituents(ho)
+        call crystallite_orientations(co,ip,el)
+      enddo
+      call mechanical_homogenize(Delta_t,ce)
+    enddo IpLooping3
+  enddo elementLooping3
+  !$OMP END PARALLEL DO
+
+
+end subroutine homogenization_mechanical_response2
 
 
 !--------------------------------------------------------------------------------------------------

src/homogenization_damage.f90

@@ -37,8 +37,7 @@ module subroutine damage_init()
   class(tNode), pointer :: &
     configHomogenizations, &
     configHomogenization, &
-    configHomogenizationDamage, &
+    configHomogenizationDamage
-    num_generic
   integer :: ho,Nmembers
 
 
@@ -70,11 +69,6 @@ module subroutine damage_init()
     end associate
   enddo
 
-!------------------------------------------------------------------------------------
-! read numerics parameter
-  num_generic => config_numerics%get('generic',defaultVal= emptyDict)
-  num_damage%charLength = num_generic%get_asFloat('charLength',defaultVal=1.0_pReal)
-
   call pass_init()
 
 end subroutine damage_init
@@ -119,8 +113,7 @@ module function homogenization_K_phi(ce) result(K)
   real(pReal), dimension(3,3) :: K
 
 
-  K = phase_K_phi(1,ce) &
+  K = phase_K_phi(1,ce)
-    * num_damage%charLength**2
 
 end function homogenization_K_phi
 

src/homogenization_mechanical.f90

@@ -123,21 +123,21 @@ end subroutine mechanical_partition
 !--------------------------------------------------------------------------------------------------
 !> @brief Average P and dPdF from the individual constituents.
 !--------------------------------------------------------------------------------------------------
-module subroutine mechanical_homogenize(dt,ce)
+module subroutine mechanical_homogenize(Delta_t,ce)
 
-  real(pReal), intent(in) :: dt
+  real(pReal), intent(in) :: Delta_t
   integer, intent(in) :: ce
 
   integer :: co
 
 
   homogenization_P(1:3,1:3,ce)            = phase_P(1,ce)
-  homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = phase_mechanical_dPdF(dt,1,ce)
+  homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = phase_mechanical_dPdF(Delta_t,1,ce)
   do co = 2, homogenization_Nconstituents(material_homogenizationID(ce))
     homogenization_P(1:3,1:3,ce)            = homogenization_P(1:3,1:3,ce) &
                                             + phase_P(co,ce)
     homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = homogenization_dPdF(1:3,1:3,1:3,1:3,ce) &
-                                            + phase_mechanical_dPdF(dt,co,ce)
+                                            + phase_mechanical_dPdF(Delta_t,co,ce)
   enddo
 
   homogenization_P(1:3,1:3,ce)            = homogenization_P(1:3,1:3,ce) &

src/lattice.f90

@@ -3,7 +3,7 @@
 !> @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  contains lattice structure definitions including Schmid matrices for slip, twin, trans,
+!> @brief  contains lattice definitions including Schmid matrices for slip, twin, trans,
 !          and cleavage as well as interaction among the various systems
 !--------------------------------------------------------------------------------------------------
 module lattice
@@ -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
+  call selfTest
-
-  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))<tol_math_check) &
-        call IO_error(135,el=i,ip=ph,ext_msg='matrix diagonal "el"ement of phase "ip"')
-    enddo
-! SHOULD NOT BE PART OF LATTICE END
-
-    call selfTest
-
-  enddo
 
 end subroutine lattice_init
 
@@ -468,10 +419,10 @@ end subroutine lattice_init
 !--------------------------------------------------------------------------------------------------
 !> @brief Characteristic shear for twinning
 !--------------------------------------------------------------------------------------------------
-function lattice_characteristicShear_Twin(Ntwin,structure,CoverA) result(characteristicShear)
+function lattice_characteristicShear_Twin(Ntwin,lattice,CoverA) result(characteristicShear)
 
   integer,     dimension(:),            intent(in) :: Ntwin                                         !< number of active twin systems per family
-  character(len=*),                     intent(in) :: structure                                     !< lattice structure
+  character(len=2),                     intent(in) :: lattice                                       !< Bravais lattice (Pearson symbol)
   real(pReal),                          intent(in) :: cOverA                                        !< c/a ratio
   real(pReal), dimension(sum(Ntwin))               :: characteristicShear
 
@@ -513,7 +464,7 @@ function lattice_characteristicShear_Twin(Ntwin,structure,CoverA) result(charact
   myFamilies: do f = 1,size(Ntwin,1)
     mySystems: do s = 1,Ntwin(f)
       a = a + 1
-      select case(structure)
+      select case(lattice)
         case('cF','cI')
           characteristicShear(a) = 0.5_pReal*sqrt(2.0_pReal)
         case('hP')
@@ -531,7 +482,7 @@ function lattice_characteristicShear_Twin(Ntwin,structure,CoverA) result(charact
               characteristicShear(a) = 2.0_pReal*(cOverA**2.0_pReal-2.0_pReal)/3.0_pReal/cOverA
           end select
         case default
-          call IO_error(137,ext_msg='lattice_characteristicShear_Twin: '//trim(structure))
+          call IO_error(137,ext_msg='lattice_characteristicShear_Twin: '//trim(lattice))
       end select
     enddo mySystems
   enddo myFamilies
@@ -542,10 +493,10 @@ end function lattice_characteristicShear_Twin
 !--------------------------------------------------------------------------------------------------
 !> @brief Rotated elasticity matrices for twinning in 66-vector notation
 !--------------------------------------------------------------------------------------------------
-function lattice_C66_twin(Ntwin,C66,structure,CoverA)
+function lattice_C66_twin(Ntwin,C66,lattice,CoverA)
 
   integer,     dimension(:),            intent(in) :: Ntwin                                         !< number of active twin systems per family
-  character(len=*),                     intent(in) :: structure                                     !< lattice structure
+  character(len=2),                     intent(in) :: lattice                                       !< Bravais lattice (Pearson symbol)
   real(pReal), dimension(6,6),          intent(in) :: C66                                           !< unrotated parent stiffness matrix
   real(pReal),                          intent(in) :: cOverA                                        !< c/a ratio
   real(pReal), dimension(6,6,sum(Ntwin))           :: lattice_C66_twin
@@ -554,18 +505,18 @@ function lattice_C66_twin(Ntwin,C66,structure,CoverA)
   type(rotation)                        :: R
   integer                               :: i
 
-  select case(structure)
+  select case(lattice)
     case('cF')
       coordinateSystem = buildCoordinateSystem(Ntwin,FCC_NSLIPSYSTEM,FCC_SYSTEMTWIN,&
-                                               structure,0.0_pReal)
+                                               lattice,0.0_pReal)
     case('cI')
       coordinateSystem = buildCoordinateSystem(Ntwin,BCC_NSLIPSYSTEM,BCC_SYSTEMTWIN,&
-                                               structure,0.0_pReal)
+                                               lattice,0.0_pReal)
     case('hP')
       coordinateSystem = buildCoordinateSystem(Ntwin,HEX_NSLIPSYSTEM,HEX_SYSTEMTWIN,&
-                                               structure,cOverA)
+                                               lattice,cOverA)
     case default
-      call IO_error(137,ext_msg='lattice_C66_twin: '//trim(structure))
+      call IO_error(137,ext_msg='lattice_C66_twin: '//trim(lattice))
   end select
 
   do i = 1, sum(Ntwin)
@@ -579,11 +530,11 @@ end function lattice_C66_twin
 !--------------------------------------------------------------------------------------------------
 !> @brief Rotated elasticity matrices for transformation in 66-vector notation
 !--------------------------------------------------------------------------------------------------
-function lattice_C66_trans(Ntrans,C_parent66,structure_target, &
+function lattice_C66_trans(Ntrans,C_parent66,lattice_target, &
                            cOverA_trans,a_bcc,a_fcc)
 
   integer,     dimension(:),             intent(in) :: Ntrans                                       !< number of active twin systems per family
-  character(len=*),                      intent(in) :: structure_target                             !< lattice structure
+  character(len=2),                      intent(in) :: lattice_target                               !< Bravais lattice (Pearson symbol)
   real(pReal), dimension(6,6),           intent(in) :: C_parent66
   real(pReal), dimension(6,6,sum(Ntrans))           :: lattice_C66_trans
 
@@ -595,9 +546,9 @@ function lattice_C66_trans(Ntrans,C_parent66,structure_target, &
 
  !--------------------------------------------------------------------------------------------------
  ! elasticity matrix of the target phase in cube orientation
-  if (structure_target == 'hP') then
+  if (lattice_target == 'hP') then
     if (cOverA_trans < 1.0_pReal .or. cOverA_trans > 2.0_pReal) &
-      call IO_error(131,ext_msg='lattice_C66_trans: '//trim(structure_target))
+      call IO_error(131,ext_msg='lattice_C66_trans: '//trim(lattice_target))
     C_bar66(1,1) = (C_parent66(1,1) + C_parent66(1,2) + 2.0_pReal*C_parent66(4,4))/2.0_pReal
     C_bar66(1,2) = (C_parent66(1,1) + 5.0_pReal*C_parent66(1,2) - 2.0_pReal*C_parent66(4,4))/6.0_pReal
     C_bar66(3,3) = (C_parent66(1,1) + 2.0_pReal*C_parent66(1,2) + 4.0_pReal*C_parent66(4,4))/3.0_pReal
@@ -611,13 +562,13 @@ function lattice_C66_trans(Ntrans,C_parent66,structure_target, &
     C_target_unrotated66(1,3) = C_bar66(1,3)
     C_target_unrotated66(3,3) = C_bar66(3,3)
     C_target_unrotated66(4,4) = C_bar66(4,4) - C_bar66(1,4)**2.0_pReal/(0.5_pReal*(C_bar66(1,1) - C_bar66(1,2)))
-    C_target_unrotated66 = applyLatticeSymmetryC66(C_target_unrotated66,'hP')
+    C_target_unrotated66 = lattice_symmetrize_C66(C_target_unrotated66,'hP')
-  elseif (structure_target  == 'cI') then
+  elseif (lattice_target  == 'cI') then
     if (a_bcc <= 0.0_pReal .or. a_fcc <= 0.0_pReal) &
-      call IO_error(134,ext_msg='lattice_C66_trans: '//trim(structure_target))
+      call IO_error(134,ext_msg='lattice_C66_trans: '//trim(lattice_target))
     C_target_unrotated66 = C_parent66
   else
-    call IO_error(137,ext_msg='lattice_C66_trans : '//trim(structure_target))
+    call IO_error(137,ext_msg='lattice_C66_trans : '//trim(lattice_target))
   endif
 
   do i = 1, 6
@@ -686,11 +637,11 @@ end function lattice_nonSchmidMatrix
 !> @brief Slip-slip interaction matrix
 !> details only active slip systems are considered
 !--------------------------------------------------------------------------------------------------
-function lattice_interaction_SlipBySlip(Nslip,interactionValues,structure) result(interactionMatrix)
+function lattice_interaction_SlipBySlip(Nslip,interactionValues,lattice) result(interactionMatrix)
 
   integer,         dimension(:),                   intent(in) :: Nslip                              !< number of active slip systems per family
   real(pReal),     dimension(:),                   intent(in) :: interactionValues                  !< values for slip-slip interaction
-  character(len=*),                                intent(in) :: structure                          !< lattice structure
+  character(len=2),                                intent(in) :: lattice                            !< Bravais lattice (Pearson symbol)
   real(pReal),     dimension(sum(Nslip),sum(Nslip))           :: interactionMatrix
 
   integer, dimension(:),   allocatable :: NslipMax
@@ -908,7 +859,7 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,structure) resul
       ],shape(BCT_INTERACTIONSLIPSLIP))
 
 
-  select case(structure)
+  select case(lattice)
     case('cF')
       interactionTypes = FCC_INTERACTIONSLIPSLIP
       NslipMax         = FCC_NSLIPSYSTEM
@@ -922,7 +873,7 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,structure) resul
       interactionTypes = BCT_INTERACTIONSLIPSLIP
       NslipMax         = BCT_NSLIPSYSTEM
     case default
-      call IO_error(137,ext_msg='lattice_interaction_SlipBySlip: '//trim(structure))
+      call IO_error(137,ext_msg='lattice_interaction_SlipBySlip: '//trim(lattice))
   end select
 
   interactionMatrix = buildInteraction(Nslip,Nslip,NslipMax,NslipMax,interactionValues,interactionTypes)
@@ -934,11 +885,11 @@ end function lattice_interaction_SlipBySlip
 !> @brief Twin-twin interaction matrix
 !> details only active twin systems are considered
 !--------------------------------------------------------------------------------------------------
-function lattice_interaction_TwinByTwin(Ntwin,interactionValues,structure) result(interactionMatrix)
+function lattice_interaction_TwinByTwin(Ntwin,interactionValues,lattice) result(interactionMatrix)
 
   integer,         dimension(:),                   intent(in) :: Ntwin                              !< number of active twin systems per family
   real(pReal),     dimension(:),                   intent(in) :: interactionValues                  !< values for twin-twin interaction
-  character(len=*),                                intent(in) :: structure                          !< lattice structure
+  character(len=2),                                intent(in) :: lattice                            !< Bravais lattice (Pearson symbol)
   real(pReal),     dimension(sum(Ntwin),sum(Ntwin))           :: interactionMatrix
 
   integer, dimension(:),   allocatable :: NtwinMax
@@ -1009,7 +960,7 @@ function lattice_interaction_TwinByTwin(Ntwin,interactionValues,structure) resul
       20,20,20,20,20,20,  19,19,19,19,19,19,  18,18,18,18,18,18,  17,17,17,17,17,16  &
       ],shape(HEX_INTERACTIONTWINTWIN))                                                             !< Twin-twin interaction types for hex
 
-  select case(structure)
+  select case(lattice)
     case('cF')
       interactionTypes = FCC_INTERACTIONTWINTWIN
       NtwinMax         = FCC_NTWINSYSTEM
@@ -1020,7 +971,7 @@ function lattice_interaction_TwinByTwin(Ntwin,interactionValues,structure) resul
       interactionTypes = HEX_INTERACTIONTWINTWIN
       NtwinMax         = HEX_NTWINSYSTEM
     case default
-      call IO_error(137,ext_msg='lattice_interaction_TwinByTwin: '//trim(structure))
+      call IO_error(137,ext_msg='lattice_interaction_TwinByTwin: '//trim(lattice))
   end select
 
   interactionMatrix = buildInteraction(Ntwin,Ntwin,NtwinMax,NtwinMax,interactionValues,interactionTypes)
@@ -1032,11 +983,11 @@ end function lattice_interaction_TwinByTwin
 !> @brief Trans-trans interaction matrix
 !> details only active trans systems are considered
 !--------------------------------------------------------------------------------------------------
-function lattice_interaction_TransByTrans(Ntrans,interactionValues,structure) result(interactionMatrix)
+function lattice_interaction_TransByTrans(Ntrans,interactionValues,lattice) result(interactionMatrix)
 
   integer,         dimension(:),                     intent(in) :: Ntrans                           !< number of active trans systems per family
   real(pReal),     dimension(:),                     intent(in) :: interactionValues                !< values for trans-trans interaction
-  character(len=*),                                  intent(in) :: structure                        !< lattice structure (parent crystal)
+  character(len=2),                                  intent(in) :: lattice                          !<Bravais lattice (Pearson symbol) (parent crystal)
   real(pReal),     dimension(sum(Ntrans),sum(Ntrans))           :: interactionMatrix
 
   integer, dimension(:),   allocatable :: NtransMax
@@ -1058,11 +1009,11 @@ function lattice_interaction_TransByTrans(Ntrans,interactionValues,structure) re
       2,2,2,2,2,2,2,2,2,1,1,1  &
       ],shape(FCC_INTERACTIONTRANSTRANS))                                                           !< Trans-trans interaction types for fcc
 
-  if(structure == 'cF') then
+  if(lattice == 'cF') then
     interactionTypes = FCC_INTERACTIONTRANSTRANS
     NtransMax        = FCC_NTRANSSYSTEM
   else
-    call IO_error(137,ext_msg='lattice_interaction_TransByTrans: '//trim(structure))
+    call IO_error(137,ext_msg='lattice_interaction_TransByTrans: '//trim(lattice))
   end if
 
   interactionMatrix = buildInteraction(Ntrans,Ntrans,NtransMax,NtransMax,interactionValues,interactionTypes)
@@ -1074,12 +1025,12 @@ end function lattice_interaction_TransByTrans
 !> @brief Slip-twin interaction matrix
 !> details only active slip and twin systems are considered
 !--------------------------------------------------------------------------------------------------
-function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,structure) result(interactionMatrix)
+function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,lattice) result(interactionMatrix)
 
   integer,         dimension(:),                   intent(in) :: Nslip, &                           !< number of active slip systems per family
                                                                  Ntwin                              !< number of active twin systems per family
   real(pReal),     dimension(:),                   intent(in) :: interactionValues                  !< values for slip-twin interaction
-  character(len=*),                                intent(in) :: structure                          !< lattice structure
+  character(len=2),                                intent(in) :: lattice                            !< Bravais lattice (Pearson symbol)
   real(pReal),     dimension(sum(Nslip),sum(Ntwin))           :: interactionMatrix
 
   integer, dimension(:),   allocatable :: NslipMax, &
@@ -1211,7 +1162,7 @@ function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,structure)
       21,21,21,21,21,21,  22,22,22,22,22,22,  23,23,23,23,23,23,  24,24,24,24,24,24  &
       ],shape(HEX_INTERACTIONSLIPTWIN))                                                             !< Slip-twin interaction types for hex
 
-  select case(structure)
+  select case(lattice)
     case('cF')
       interactionTypes = FCC_INTERACTIONSLIPTWIN
       NslipMax         = FCC_NSLIPSYSTEM
@@ -1225,7 +1176,7 @@ function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,structure)
       NslipMax         = HEX_NSLIPSYSTEM
       NtwinMax         = HEX_NTWINSYSTEM
     case default
-      call IO_error(137,ext_msg='lattice_interaction_SlipByTwin: '//trim(structure))
+      call IO_error(137,ext_msg='lattice_interaction_SlipByTwin: '//trim(lattice))
   end select
 
   interactionMatrix = buildInteraction(Nslip,Ntwin,NslipMax,NtwinMax,interactionValues,interactionTypes)
@@ -1237,12 +1188,12 @@ end function lattice_interaction_SlipByTwin
 !> @brief Slip-trans interaction matrix
 !> details only active slip and trans systems are considered
 !--------------------------------------------------------------------------------------------------
-function lattice_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,structure) result(interactionMatrix)
+function lattice_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,lattice) result(interactionMatrix)
 
   integer,         dimension(:),                    intent(in) :: Nslip, &                          !< number of active slip systems per family
                                                                   Ntrans                            !< number of active trans systems per family
   real(pReal),     dimension(:),                    intent(in) :: interactionValues                 !< values for slip-trans interaction
-  character(len=*),                                 intent(in) :: structure                         !< lattice structure (parent crystal)
+  character(len=2),                                 intent(in) :: lattice                           !< Bravais lattice (Pearson symbol) (parent crystal)
   real(pReal),     dimension(sum(Nslip),sum(Ntrans))           :: interactionMatrix
 
   integer, dimension(:),   allocatable :: NslipMax, &
@@ -1272,13 +1223,13 @@ function lattice_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,structur
       4,4,4,4,4,4,4,4,4,4,4,4  &
       ],shape(FCC_INTERACTIONSLIPTRANS))                                                            !< Slip-trans interaction types for fcc
 
-  select case(structure)
+  select case(lattice)
     case('cF')
       interactionTypes = FCC_INTERACTIONSLIPTRANS
       NslipMax         = FCC_NSLIPSYSTEM
       NtransMax        = FCC_NTRANSSYSTEM
     case default
-      call IO_error(137,ext_msg='lattice_interaction_SlipByTrans: '//trim(structure))
+      call IO_error(137,ext_msg='lattice_interaction_SlipByTrans: '//trim(lattice))
   end select
 
   interactionMatrix = buildInteraction(Nslip,Ntrans,NslipMax,NtransMax,interactionValues,interactionTypes)
@@ -1290,12 +1241,12 @@ function lattice_interaction_SlipByTrans(Nslip,Ntrans,interactionValues,structur
 !> @brief Twin-slip interaction matrix
 !> details only active twin and slip systems are considered
 !--------------------------------------------------------------------------------------------------
-function lattice_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,structure) result(interactionMatrix)
+function lattice_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,lattice) result(interactionMatrix)
 
   integer,         dimension(:),                   intent(in) :: Ntwin, &                           !< number of active twin systems per family
                                                                  Nslip                              !< number of active slip systems per family
   real(pReal),     dimension(:),                   intent(in) :: interactionValues                  !< values for twin-twin interaction
-  character(len=*),                                intent(in) :: structure                          !< lattice structure
+  character(len=2),                                intent(in) :: lattice                            !< Bravais lattice (Pearson symbol)
   real(pReal),     dimension(sum(Ntwin),sum(Nslip))           :: interactionMatrix
 
   integer, dimension(:),   allocatable :: NtwinMax, &
@@ -1339,7 +1290,7 @@ function lattice_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,structure)
       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  &
       ],shape(HEX_INTERACTIONTWINSLIP))                                                             !< Twin-slip interaction types for hex
 
-  select case(structure)
+  select case(lattice)
     case('cF')
       interactionTypes = FCC_INTERACTIONTWINSLIP
       NtwinMax         = FCC_NTWINSYSTEM
@@ -1353,7 +1304,7 @@ function lattice_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,structure)
       NtwinMax         = HEX_NTWINSYSTEM
       NslipMax         = HEX_NSLIPSYSTEM
     case default
-      call IO_error(137,ext_msg='lattice_interaction_TwinBySlip: '//trim(structure))
+      call IO_error(137,ext_msg='lattice_interaction_TwinBySlip: '//trim(lattice))
   end select
 
   interactionMatrix = buildInteraction(Ntwin,Nslip,NtwinMax,NslipMax,interactionValues,interactionTypes)
@@ -1365,10 +1316,10 @@ end function lattice_interaction_TwinBySlip
 !> @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,lattice,cOverA) result(SchmidMatrix)
 
   integer,         dimension(:),            intent(in) :: Nslip                                     !< number of active slip systems per family
-  character(len=*),                         intent(in) :: structure                                 !< lattice structure
+  character(len=2),                         intent(in) :: lattice                                   !< Bravais lattice (Pearson symbol)
   real(pReal),                              intent(in) :: cOverA
   real(pReal),     dimension(3,3,sum(Nslip))           :: SchmidMatrix
 
@@ -1377,7 +1328,7 @@ function lattice_SchmidMatrix_slip(Nslip,structure,cOverA) result(SchmidMatrix)
   integer,     dimension(:),             allocatable :: NslipMax
   integer                                            :: i
 
-  select case(structure)
+  select case(lattice)
     case('cF')
       NslipMax    = FCC_NSLIPSYSTEM
       slipSystems = FCC_SYSTEMSLIP
@@ -1392,15 +1343,15 @@ function lattice_SchmidMatrix_slip(Nslip,structure,cOverA) result(SchmidMatrix)
       slipSystems = BCT_SYSTEMSLIP
     case default
       allocate(NslipMax(0))
-      call IO_error(137,ext_msg='lattice_SchmidMatrix_slip: '//trim(structure))
+      call IO_error(137,ext_msg='lattice_SchmidMatrix_slip: '//trim(lattice))
   end select
 
   if (any(NslipMax(1:size(Nslip)) - Nslip < 0)) &
-    call IO_error(145,ext_msg='Nslip '//trim(structure))
+    call IO_error(145,ext_msg='Nslip '//trim(lattice))
   if (any(Nslip < 0)) &
-    call IO_error(144,ext_msg='Nslip '//trim(structure))
+    call IO_error(144,ext_msg='Nslip '//trim(lattice))
 
-  coordinateSystem = buildCoordinateSystem(Nslip,NslipMax,slipSystems,structure,cOverA)
+  coordinateSystem = buildCoordinateSystem(Nslip,NslipMax,slipSystems,lattice,cOverA)
 
   do i = 1, sum(Nslip)
     SchmidMatrix(1:3,1:3,i) = math_outer(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i))
@@ -1415,10 +1366,10 @@ end function lattice_SchmidMatrix_slip
 !> @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,lattice,cOverA) result(SchmidMatrix)
 
   integer,         dimension(:),            intent(in) :: Ntwin                                     !< number of active twin systems per family
-  character(len=*),                         intent(in) :: structure                                 !< lattice structure
+  character(len=2),                         intent(in) :: lattice                                   !< Bravais lattice (Pearson symbol)
   real(pReal),                              intent(in) :: cOverA                                    !< c/a ratio
   real(pReal),     dimension(3,3,sum(Ntwin))           :: SchmidMatrix
 
@@ -1427,7 +1378,7 @@ function lattice_SchmidMatrix_twin(Ntwin,structure,cOverA) result(SchmidMatrix)
   integer,     dimension(:),             allocatable :: NtwinMax
   integer                                            :: i
 
-  select case(structure)
+  select case(lattice)
     case('cF')
       NtwinMax    = FCC_NTWINSYSTEM
       twinSystems = FCC_SYSTEMTWIN
@@ -1439,15 +1390,15 @@ function lattice_SchmidMatrix_twin(Ntwin,structure,cOverA) result(SchmidMatrix)
       twinSystems = HEX_SYSTEMTWIN
     case default
       allocate(NtwinMax(0))
-      call IO_error(137,ext_msg='lattice_SchmidMatrix_twin: '//trim(structure))
+      call IO_error(137,ext_msg='lattice_SchmidMatrix_twin: '//trim(lattice))
   end select
 
   if (any(NtwinMax(1:size(Ntwin)) - Ntwin < 0)) &
-    call IO_error(145,ext_msg='Ntwin '//trim(structure))
+    call IO_error(145,ext_msg='Ntwin '//trim(lattice))
   if (any(Ntwin < 0)) &
-    call IO_error(144,ext_msg='Ntwin '//trim(structure))
+    call IO_error(144,ext_msg='Ntwin '//trim(lattice))
 
-  coordinateSystem = buildCoordinateSystem(Ntwin,NtwinMax,twinSystems,structure,cOverA)
+  coordinateSystem = buildCoordinateSystem(Ntwin,NtwinMax,twinSystems,lattice,cOverA)
 
   do i = 1, sum(Ntwin)
     SchmidMatrix(1:3,1:3,i) = math_outer(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i))
@@ -1462,24 +1413,24 @@ end function lattice_SchmidMatrix_twin
 !> @brief Schmid matrix for twinning
 !> details only active twin systems are considered
 !--------------------------------------------------------------------------------------------------
-function lattice_SchmidMatrix_trans(Ntrans,structure_target,cOverA,a_bcc,a_fcc) result(SchmidMatrix)
+function lattice_SchmidMatrix_trans(Ntrans,lattice_target,cOverA,a_bcc,a_fcc) result(SchmidMatrix)
 
   integer,         dimension(:),             intent(in) :: Ntrans                                   !< number of active twin systems per family
-  character(len=*),                          intent(in) :: structure_target                         !< lattice structure
+  character(len=2),                          intent(in) :: lattice_target                           !< Bravais lattice (Pearson symbol)
   real(pReal),                               intent(in) :: cOverA                                   !< c/a ratio
   real(pReal),     dimension(3,3,sum(Ntrans))           :: SchmidMatrix
 
   real(pReal), dimension(3,3,sum(Ntrans)) :: devNull
   real(pReal)                             :: a_bcc, a_fcc
 
-  if (structure_target /= 'cI' .and. structure_target /= 'hP') &
+  if (lattice_target /= 'cI' .and. lattice_target /= 'hP') &
-    call IO_error(137,ext_msg='lattice_SchmidMatrix_trans: '//trim(structure_target))
+    call IO_error(137,ext_msg='lattice_SchmidMatrix_trans: '//trim(lattice_target))
 
-  if (structure_target == 'hP' .and. (cOverA < 1.0_pReal .or. cOverA > 2.0_pReal)) &
+  if (lattice_target == 'hP' .and. (cOverA < 1.0_pReal .or. cOverA > 2.0_pReal)) &
-    call IO_error(131,ext_msg='lattice_SchmidMatrix_trans: '//trim(structure_target))
+    call IO_error(131,ext_msg='lattice_SchmidMatrix_trans: '//trim(lattice_target))
 
-  if (structure_target == 'cI' .and. (a_bcc <= 0.0_pReal .or. a_fcc <= 0.0_pReal)) &
+  if (lattice_target == 'cI' .and. (a_bcc <= 0.0_pReal .or. a_fcc <= 0.0_pReal)) &
-    call IO_error(134,ext_msg='lattice_SchmidMatrix_trans: '//trim(structure_target))
+    call IO_error(134,ext_msg='lattice_SchmidMatrix_trans: '//trim(lattice_target))
 
   call buildTransformationSystem(devNull,SchmidMatrix,Ntrans,cOverA,a_fcc,a_bcc)
 
@@ -1490,10 +1441,10 @@ end function lattice_SchmidMatrix_trans
 !> @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,lattice,cOverA) result(SchmidMatrix)
 
   integer,         dimension(:),                  intent(in) :: Ncleavage                           !< number of active cleavage systems per family
-  character(len=*),                               intent(in) :: structure                           !< lattice structure
+  character(len=2),                               intent(in) :: lattice                             !< Bravais lattice (Pearson symbol)
   real(pReal),                                    intent(in) :: cOverA                              !< c/a ratio
   real(pReal),     dimension(3,3,3,sum(Ncleavage))           :: SchmidMatrix
 
@@ -1502,7 +1453,7 @@ function lattice_SchmidMatrix_cleavage(Ncleavage,structure,cOverA) result(Schmid
   integer,     dimension(:),                 allocatable :: NcleavageMax
   integer                                                :: i
 
-  select case(structure)
+  select case(lattice)
     case('cF')
       NcleavageMax    = FCC_NCLEAVAGESYSTEM
       cleavageSystems = FCC_SYSTEMCLEAVAGE
@@ -1511,15 +1462,15 @@ function lattice_SchmidMatrix_cleavage(Ncleavage,structure,cOverA) result(Schmid
       cleavageSystems = BCC_SYSTEMCLEAVAGE
     case default
       allocate(NcleavageMax(0))
-      call IO_error(137,ext_msg='lattice_SchmidMatrix_cleavage: '//trim(structure))
+      call IO_error(137,ext_msg='lattice_SchmidMatrix_cleavage: '//trim(lattice))
   end select
 
   if (any(NcleavageMax(1:size(Ncleavage)) - Ncleavage < 0)) &
-    call IO_error(145,ext_msg='Ncleavage '//trim(structure))
+    call IO_error(145,ext_msg='Ncleavage '//trim(lattice))
   if (any(Ncleavage < 0)) &
-    call IO_error(144,ext_msg='Ncleavage '//trim(structure))
+    call IO_error(144,ext_msg='Ncleavage '//trim(lattice))
 
-  coordinateSystem = buildCoordinateSystem(Ncleavage,NcleavageMax,cleavageSystems,structure,cOverA)
+  coordinateSystem = buildCoordinateSystem(Ncleavage,NcleavageMax,cleavageSystems,lattice,cOverA)
 
   do i = 1, sum(Ncleavage)
     SchmidMatrix(1:3,1:3,1,i) = math_outer(coordinateSystem(1:3,1,i),coordinateSystem(1:3,2,i))
@@ -1533,16 +1484,16 @@ end function lattice_SchmidMatrix_cleavage
 !--------------------------------------------------------------------------------------------------
 !> @brief Slip direction of slip systems (|| b)
 !--------------------------------------------------------------------------------------------------
-function lattice_slip_direction(Nslip,structure,cOverA) result(d)
+function lattice_slip_direction(Nslip,lattice,cOverA) result(d)
 
   integer,         dimension(:),           intent(in) :: Nslip                                      !< number of active slip systems per family
-  character(len=*),                        intent(in) :: structure                                  !< lattice structure
+  character(len=2),                        intent(in) :: lattice                                    !< Bravais lattice (Pearson symbol)
   real(pReal),                             intent(in) :: cOverA                                     !< c/a ratio
   real(pReal),     dimension(3,sum(Nslip))            :: d
 
   real(pReal), dimension(3,3,sum(Nslip)) :: coordinateSystem
 
-  coordinateSystem = coordinateSystem_slip(Nslip,structure,cOverA)
+  coordinateSystem = coordinateSystem_slip(Nslip,lattice,cOverA)
   d = coordinateSystem(1:3,1,1:sum(Nslip))
 
 end function lattice_slip_direction
@@ -1551,16 +1502,16 @@ end function lattice_slip_direction
 !--------------------------------------------------------------------------------------------------
 !> @brief Normal direction of slip systems (|| n)
 !--------------------------------------------------------------------------------------------------
-function lattice_slip_normal(Nslip,structure,cOverA) result(n)
+function lattice_slip_normal(Nslip,lattice,cOverA) result(n)
 
   integer,         dimension(:),           intent(in) :: Nslip                                      !< number of active slip systems per family
-  character(len=*),                        intent(in) :: structure                                  !< lattice structure
+  character(len=2),                        intent(in) :: lattice                                    !< Bravais lattice (Pearson symbol)
   real(pReal),                             intent(in) :: cOverA                                     !< c/a ratio
   real(pReal),     dimension(3,sum(Nslip))            :: n
 
   real(pReal), dimension(3,3,sum(Nslip)) :: coordinateSystem
 
-  coordinateSystem = coordinateSystem_slip(Nslip,structure,cOverA)
+  coordinateSystem = coordinateSystem_slip(Nslip,lattice,cOverA)
   n = coordinateSystem(1:3,2,1:sum(Nslip))
 
 end function lattice_slip_normal
@@ -1569,16 +1520,16 @@ end function lattice_slip_normal
 !--------------------------------------------------------------------------------------------------
 !> @brief Transverse direction of slip systems (|| t = b x n)
 !--------------------------------------------------------------------------------------------------
-function lattice_slip_transverse(Nslip,structure,cOverA) result(t)
+function lattice_slip_transverse(Nslip,lattice,cOverA) result(t)
 
   integer,         dimension(:),           intent(in) :: Nslip                                      !< number of active slip systems per family
-  character(len=*),                        intent(in) :: structure                                  !< lattice structure
+  character(len=2),                        intent(in) :: lattice                                    !< Bravais lattice (Pearson symbol)
   real(pReal),                             intent(in) :: cOverA                                     !< c/a ratio
   real(pReal),     dimension(3,sum(Nslip))            :: t
 
   real(pReal), dimension(3,3,sum(Nslip)) :: coordinateSystem
 
-  coordinateSystem = coordinateSystem_slip(Nslip,structure,cOverA)
+  coordinateSystem = coordinateSystem_slip(Nslip,lattice,cOverA)
   t = coordinateSystem(1:3,3,1:sum(Nslip))
 
 end function lattice_slip_transverse
@@ -1588,17 +1539,17 @@ end function lattice_slip_transverse
 !> @brief Labels for slip systems
 !> details only active slip systems are considered
 !--------------------------------------------------------------------------------------------------
-function lattice_labels_slip(Nslip,structure) result(labels)
+function lattice_labels_slip(Nslip,lattice) result(labels)
 
   integer,         dimension(:),  intent(in)  :: Nslip                                              !< number of active slip systems per family
-  character(len=*),               intent(in)  :: structure                                          !< lattice structure
+  character(len=2),               intent(in)  :: lattice                                            !< Bravais lattice (Pearson symbol)
 
   character(len=:), dimension(:),   allocatable :: labels
 
   real(pReal),      dimension(:,:), allocatable :: slipSystems
   integer,          dimension(:),   allocatable :: NslipMax
 
-  select case(structure)
+  select case(lattice)
     case('cF')
       NslipMax    = FCC_NSLIPSYSTEM
       slipSystems = FCC_SYSTEMSLIP
@@ -1612,13 +1563,13 @@ function lattice_labels_slip(Nslip,structure) result(labels)
       NslipMax    = BCT_NSLIPSYSTEM
       slipSystems = BCT_SYSTEMSLIP
     case default
-      call IO_error(137,ext_msg='lattice_labels_slip: '//trim(structure))
+      call IO_error(137,ext_msg='lattice_labels_slip: '//trim(lattice))
   end select
 
   if (any(NslipMax(1:size(Nslip)) - Nslip < 0)) &
-    call IO_error(145,ext_msg='Nslip '//trim(structure))
+    call IO_error(145,ext_msg='Nslip '//trim(lattice))
   if (any(Nslip < 0)) &
-    call IO_error(144,ext_msg='Nslip '//trim(structure))
+    call IO_error(144,ext_msg='Nslip '//trim(lattice))
 
   labels = getLabels(Nslip,NslipMax,slipSystems)
 
@@ -1626,19 +1577,19 @@ end function lattice_labels_slip
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Return 3x3 tensor with symmetry according to given crystal structure
+!> @brief Return 3x3 tensor with symmetry according to given Bravais lattice
 !--------------------------------------------------------------------------------------------------
-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                                                !< Bravais lattice (Pearson symbol)
 
 
   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
+!> @brief Return stiffness matrix in 6x6 notation with symmetry according to given Bravais lattice
 !> @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=2),            intent(in) :: lattice                                                !< Bravais lattice (Pearson symbol)
 
   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,24 +1646,24 @@ pure function applyLatticeSymmetryC66(C66,structure) result(C66_sym)
      enddo
    enddo
 
-end function applyLatticeSymmetryC66
+end function lattice_symmetrize_C66
 
 
 !--------------------------------------------------------------------------------------------------
 !> @brief Labels for twin systems
 !> details only active twin systems are considered
 !--------------------------------------------------------------------------------------------------
-function lattice_labels_twin(Ntwin,structure) result(labels)
+function lattice_labels_twin(Ntwin,lattice) result(labels)
 
   integer,         dimension(:),  intent(in)  :: Ntwin                                              !< number of active slip systems per family
-  character(len=*),               intent(in)  :: structure                                          !< lattice structure
+  character(len=2),               intent(in)  :: lattice                                            !< Bravais lattice (Pearson symbol)
 
   character(len=:), dimension(:),   allocatable :: labels
 
   real(pReal),      dimension(:,:), allocatable :: twinSystems
   integer,          dimension(:),   allocatable :: NtwinMax
 
-  select case(structure)
+  select case(lattice)
     case('cF')
       NtwinMax    = FCC_NTWINSYSTEM
       twinSystems = FCC_SYSTEMTWIN
@@ -1723,13 +1674,13 @@ function lattice_labels_twin(Ntwin,structure) result(labels)
       NtwinMax    = HEX_NTWINSYSTEM
       twinSystems = HEX_SYSTEMTWIN
     case default
-      call IO_error(137,ext_msg='lattice_labels_twin: '//trim(structure))
+      call IO_error(137,ext_msg='lattice_labels_twin: '//trim(lattice))
   end select
 
   if (any(NtwinMax(1:size(Ntwin)) - Ntwin < 0)) &
-    call IO_error(145,ext_msg='Ntwin '//trim(structure))
+    call IO_error(145,ext_msg='Ntwin '//trim(lattice))
   if (any(Ntwin < 0)) &
-    call IO_error(144,ext_msg='Ntwin '//trim(structure))
+    call IO_error(144,ext_msg='Ntwin '//trim(lattice))
 
   labels = getLabels(Ntwin,NtwinMax,twinSystems)
 
@@ -1740,18 +1691,18 @@ end function lattice_labels_twin
 !> @brief Projection of the transverse direction onto the slip plane
 !> @details: This projection is used to calculate forest hardening for edge dislocations
 !--------------------------------------------------------------------------------------------------
-function slipProjection_transverse(Nslip,structure,cOverA) result(projection)
+function slipProjection_transverse(Nslip,lattice,cOverA) result(projection)
 
   integer,         dimension(:),                   intent(in) :: Nslip                              !< number of active slip systems per family
-  character(len=*),                                intent(in) :: structure                          !< lattice structure
+  character(len=2),                                intent(in) :: lattice                            !< Bravais lattice (Pearson symbol)
   real(pReal),                                     intent(in) :: cOverA                             !< c/a ratio
   real(pReal),     dimension(sum(Nslip),sum(Nslip))           :: projection
 
   real(pReal), dimension(3,sum(Nslip)) :: n, t
   integer                              :: i, j
 
-  n = lattice_slip_normal    (Nslip,structure,cOverA)
+  n = lattice_slip_normal    (Nslip,lattice,cOverA)
-  t = lattice_slip_transverse(Nslip,structure,cOverA)
+  t = lattice_slip_transverse(Nslip,lattice,cOverA)
 
   do i=1, sum(Nslip); do j=1, sum(Nslip)
     projection(i,j) = abs(math_inner(n(:,i),t(:,j)))
@@ -1764,18 +1715,18 @@ end function slipProjection_transverse
 !> @brief Projection of the slip direction onto the slip plane
 !> @details: This projection is used to calculate forest hardening for screw dislocations
 !--------------------------------------------------------------------------------------------------
-function slipProjection_direction(Nslip,structure,cOverA) result(projection)
+function slipProjection_direction(Nslip,lattice,cOverA) result(projection)
 
   integer,         dimension(:),                   intent(in) :: Nslip                              !< number of active slip systems per family
-  character(len=*),                                intent(in) :: structure                          !< lattice structure
+  character(len=2),                                intent(in) :: lattice                            !< Bravais lattice (Pearson symbol)
   real(pReal),                                     intent(in) :: cOverA                             !< c/a ratio
   real(pReal),     dimension(sum(Nslip),sum(Nslip))           :: projection
 
   real(pReal), dimension(3,sum(Nslip)) :: n, d
   integer                              :: i, j
 
-  n = lattice_slip_normal   (Nslip,structure,cOverA)
+  n = lattice_slip_normal   (Nslip,lattice,cOverA)
-  d = lattice_slip_direction(Nslip,structure,cOverA)
+  d = lattice_slip_direction(Nslip,lattice,cOverA)
 
   do i=1, sum(Nslip); do j=1, sum(Nslip)
     projection(i,j) = abs(math_inner(n(:,i),d(:,j)))
@@ -1788,17 +1739,17 @@ end function slipProjection_direction
 !> @brief build a local coordinate system on slip systems
 !> @details Order: Direction, plane (normal), and common perpendicular
 !--------------------------------------------------------------------------------------------------
-function coordinateSystem_slip(Nslip,structure,cOverA) result(coordinateSystem)
+function coordinateSystem_slip(Nslip,lattice,cOverA) result(coordinateSystem)
 
   integer,          dimension(:),            intent(in) :: Nslip                                    !< number of active slip systems per family
-  character(len=*),                          intent(in) :: structure                                !< lattice structure
+  character(len=2),                          intent(in) :: lattice                                  !< Bravais lattice (Pearson symbol)
   real(pReal),                               intent(in) :: cOverA                                   !< c/a ratio
   real(pReal),     dimension(3,3,sum(Nslip))            :: coordinateSystem
 
   real(pReal), dimension(:,:), allocatable :: slipSystems
   integer,     dimension(:),   allocatable :: NslipMax
 
-  select case(structure)
+  select case(lattice)
     case('cF')
       NslipMax    = FCC_NSLIPSYSTEM
       slipSystems = FCC_SYSTEMSLIP
@@ -1813,15 +1764,15 @@ function coordinateSystem_slip(Nslip,structure,cOverA) result(coordinateSystem)
       slipSystems = BCT_SYSTEMSLIP
     case default
       allocate(NslipMax(0))
-      call IO_error(137,ext_msg='coordinateSystem_slip: '//trim(structure))
+      call IO_error(137,ext_msg='coordinateSystem_slip: '//trim(lattice))
   end select
 
   if (any(NslipMax(1:size(Nslip)) - Nslip < 0)) &
-    call IO_error(145,ext_msg='Nslip '//trim(structure))
+    call IO_error(145,ext_msg='Nslip '//trim(lattice))
   if (any(Nslip < 0)) &
-    call IO_error(144,ext_msg='Nslip '//trim(structure))
+    call IO_error(144,ext_msg='Nslip '//trim(lattice))
 
-  coordinateSystem = buildCoordinateSystem(Nslip,NslipMax,slipSystems,structure,cOverA)
+  coordinateSystem = buildCoordinateSystem(Nslip,NslipMax,slipSystems,lattice,cOverA)
 
 end function coordinateSystem_slip
 
@@ -1873,15 +1824,15 @@ end function buildInteraction
 !> @brief Build a local coordinate system on slip, twin, trans, cleavage systems
 !> @details Order: Direction, plane (normal), and common perpendicular
 !--------------------------------------------------------------------------------------------------
-function buildCoordinateSystem(active,potential,system,structure,cOverA)
+function buildCoordinateSystem(active,potential,system,lattice,cOverA)
 
   integer, dimension(:), intent(in) :: &
     active, &                                                                                       !< # of active systems per family
     potential                                                                                       !< # of potential systems per family
   real(pReal), dimension(:,:), intent(in) :: &
     system
-  character(len=*),            intent(in) :: &
+  character(len=2),            intent(in) :: &
-    structure                                                                                       !< lattice structure
+    lattice                                                                                         !< Bravais lattice (Pearson symbol)
   real(pReal),                 intent(in) :: &
     cOverA
   real(pReal), dimension(3,3,sum(active)) :: &
@@ -1895,10 +1846,10 @@ function buildCoordinateSystem(active,potential,system,structure,cOverA)
     f, &                                                                                            !< index of my family
     s                                                                                               !< index of my system in current family
 
-  if (structure == 'tI' .and. cOverA > 2.0_pReal) &
+  if (lattice == 'tI' .and. cOverA > 2.0_pReal) &
-    call IO_error(131,ext_msg='buildCoordinateSystem:'//trim(structure))
+    call IO_error(131,ext_msg='buildCoordinateSystem:'//trim(lattice))
-  if (structure == 'hP' .and. (cOverA < 1.0_pReal .or. cOverA > 2.0_pReal)) &
+  if (lattice == 'hP' .and. (cOverA < 1.0_pReal .or. cOverA > 2.0_pReal)) &
-    call IO_error(131,ext_msg='buildCoordinateSystem:'//trim(structure))
+    call IO_error(131,ext_msg='buildCoordinateSystem:'//trim(lattice))
 
   a = 0
   activeFamilies: do f = 1,size(active,1)
@@ -1906,7 +1857,7 @@ function buildCoordinateSystem(active,potential,system,structure,cOverA)
       a = a + 1
       p = sum(potential(1:f-1))+s
 
-      select case(structure)
+      select case(lattice)
 
         case ('cF','cI','tI')
           direction = system(1:3,p)
@@ -1921,7 +1872,7 @@ function buildCoordinateSystem(active,potential,system,structure,cOverA)
                         system(8,p)/cOverA ]                                                        ! plane (hkil)->(h (h+2k)/sqrt(3) l/(p/a))
 
         case default
-          call IO_error(137,ext_msg='buildCoordinateSystem: '//trim(structure))
+          call IO_error(137,ext_msg='buildCoordinateSystem: '//trim(lattice))
 
       end select
 
@@ -1950,9 +1901,9 @@ subroutine buildTransformationSystem(Q,S,Ntrans,cOverA,a_fcc,a_bcc)
     Q, &                                                                                            !< Total rotation: Q = R*B
     S                                                                                               !< Eigendeformation tensor for phase transformation
   real(pReal),                 intent(in) :: &
-    cOverA, &                                                                                       !< c/a for target hex structure
+    cOverA, &                                                                                       !< c/a for target hex lattice
-    a_bcc, &                                                                                        !< lattice parameter a for target bcc structure
+    a_bcc, &                                                                                        !< lattice parameter a for bcc target lattice
-    a_fcc                                                                                           !< lattice parameter a for parent fcc structure
+    a_fcc                                                                                           !< lattice parameter a for fcc parent lattice
 
   type(rotation) :: &
     R, &                                                                                            !< Pitsch rotation
@@ -2126,9 +2077,10 @@ end function getlabels
 function lattice_equivalent_nu(C,assumption) result(nu)
 
   real(pReal), dimension(6,6), intent(in) :: C                                                      !< Stiffness tensor (Voigt notation)
-  character(len=*),            intent(in) :: assumption                                             !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress)
+  character(len=5),            intent(in) :: assumption                                             !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress)
-  real(pReal)                 :: K, mu, nu
+  real(pReal) :: nu
 
+  real(pReal) :: K, mu
   logical                     :: error
   real(pReal), dimension(6,6) :: S
 
@@ -2158,8 +2110,8 @@ end function lattice_equivalent_nu
 function lattice_equivalent_mu(C,assumption) result(mu)
 
   real(pReal), dimension(6,6), intent(in) :: C                                                      !< Stiffness tensor (Voigt notation)
-  character(len=*),            intent(in) :: assumption                                             !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress)
+  character(len=5),            intent(in) :: assumption                                             !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress)
-  real(pReal)                 :: mu
+  real(pReal) :: mu
 
   logical                     :: error
   real(pReal), dimension(6,6) :: S
@@ -2203,10 +2155,10 @@ subroutine selfTest
 
   do i = 1, 10
     call random_number(C)
-    C_cF = applyLatticeSymmetryC66(C,'cI')
+    C_cF = lattice_symmetrize_C66(C,'cI')
-    C_cI = applyLatticeSymmetryC66(C,'cF')
+    C_cI = lattice_symmetrize_C66(C,'cF')
-    C_hP = applyLatticeSymmetryC66(C,'hP')
+    C_hP = lattice_symmetrize_C66(C,'hP')
-    C_tI = applyLatticeSymmetryC66(C,'tI')
+    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 +2178,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_cF = lattice_symmetrize_33(T,'cI')
-    T_cI = lattice_applyLatticeSymmetry33(T,'cF')
+    T_cI = lattice_symmetrize_33(T,'cF')
-    T_hP = lattice_applyLatticeSymmetry33(T,'hP')
+    T_hP = lattice_symmetrize_33(T,'hP')
-    T_tI = lattice_applyLatticeSymmetry33(T,'tI')
+    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 +2196,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'
 

src/math.f90

@@ -1068,6 +1068,7 @@ integer pure function math_factorial(n)
 
   integer, intent(in) :: n
 
+
   math_factorial = product(math_range(n))
 
 end function math_factorial
@@ -1081,6 +1082,7 @@ integer pure function math_binomial(n,k)
   integer, intent(in) :: n, k
   integer :: i, k_, n_
 
+
   k_ = min(k,n-k)
   n_ = n
   math_binomial = merge(1,0,k_>-1)                                                                  ! handling special cases k < 0 or k > n
@@ -1095,15 +1097,13 @@ end function math_binomial
 !--------------------------------------------------------------------------------------------------
 !> @brief multinomial coefficient
 !--------------------------------------------------------------------------------------------------
-integer pure function math_multinomial(alpha)
+integer pure function math_multinomial(k)
 
-  integer, intent(in), dimension(:) :: alpha
+  integer, intent(in), dimension(:) :: k
   integer :: i
 
-  math_multinomial = 1
+
-  do i = 1, size(alpha)
+  math_multinomial = product([(math_binomial(sum(k(1:i)),k(i)),i=1,size(k))])
-    math_multinomial = math_multinomial*math_binomial(sum(alpha(1:i)),alpha(i))
-  enddo
 
 end function math_multinomial
 
@@ -1116,6 +1116,7 @@ real(pReal) pure function math_volTetrahedron(v1,v2,v3,v4)
   real(pReal), dimension (3), intent(in) :: v1,v2,v3,v4
   real(pReal), dimension (3,3) :: m
 
+
   m(1:3,1) = v1-v2
   m(1:3,2) = v1-v3
   m(1:3,3) = v1-v4
@@ -1289,6 +1290,9 @@ subroutine selfTest
   if(math_binomial(49,6) /= 13983816) &
     error stop 'math_binomial'
 
+  if(math_multinomial([1,2,3,4]) /= 12600) &
+    error stop 'math_multinomial'
+
   ijk = cshift([1,2,3],int(r*1.0e2_pReal))
   if(dNeq(math_LeviCivita(ijk(1),ijk(2),ijk(3)),+1.0_pReal)) &
     error stop 'math_LeviCivita(even)'

src/mesh/FEM_utilities.f90

@@ -27,17 +27,17 @@ module FEM_utilities
   logical, public             :: cutBack = .false.                                                  !< cut back of BVP solver in case convergence is not achieved or a material point is terminally ill
   integer, public, parameter  :: maxFields = 6
   integer, public             :: nActiveFields = 0
- 
+
 !--------------------------------------------------------------------------------------------------
 ! grid related information information
   real(pReal),   public       :: wgt                                                                !< weighting factor 1/Nelems
-  
+
 
 !--------------------------------------------------------------------------------------------------
 ! field labels information
   character(len=*),                         parameter,            public :: &
     FIELD_MECH_label     = 'mechanical'
- 
+
   enum, bind(c); enumerator :: &
     FIELD_UNDEFINED_ID, &
     FIELD_MECH_ID
@@ -48,7 +48,7 @@ module FEM_utilities
     COMPONENT_MECH_Y_ID, &
     COMPONENT_MECH_Z_ID
   end enum
- 
+
 !--------------------------------------------------------------------------------------------------
 ! variables controlling debugging
  logical :: &
@@ -57,23 +57,23 @@ module FEM_utilities
 !--------------------------------------------------------------------------------------------------
 ! derived types
   type, public :: tSolutionState                                                                    !< return type of solution from FEM solver variants
-    logical :: converged        = .true.   
+    logical :: converged        = .true.
-    logical :: stagConverged    = .true.   
+    logical :: stagConverged    = .true.
     integer :: iterationsNeeded = 0
   end type tSolutionState
- 
+
   type, public :: tComponentBC
     integer(kind(COMPONENT_UNDEFINED_ID))          :: ID
     real(pReal),                       allocatable, dimension(:) :: Value
-    logical,                           allocatable, dimension(:) :: Mask 
+    logical,                           allocatable, dimension(:) :: Mask
   end type tComponentBC
- 
+
   type, public :: tFieldBC
     integer(kind(FIELD_UNDEFINED_ID))  :: ID
     integer                            :: nComponents = 0
     type(tComponentBC),    allocatable :: componentBC(:)
   end type tFieldBC
- 
+
   type, public :: tLoadCase
     real(pReal)  :: time                   = 0.0_pReal                                              !< length of increment
     integer      :: incs                   = 0, &                                                   !< number of increments
@@ -83,7 +83,7 @@ module FEM_utilities
     integer,        allocatable, dimension(:) :: faceID
     type(tFieldBC), allocatable, dimension(:) :: fieldBC
   end type tLoadCase
-  
+
   public :: &
     FEM_utilities_init, &
     utilities_constitutiveResponse, &
@@ -94,14 +94,14 @@ module FEM_utilities
     COMPONENT_MECH_Y_ID, &
     COMPONENT_MECH_Z_ID
 
-contains 
+contains
 
 !ToDo: use functions in variable call
 !--------------------------------------------------------------------------------------------------
 !> @brief allocates all neccessary fields, sets debug flags
 !--------------------------------------------------------------------------------------------------
 subroutine FEM_utilities_init
-   
+
   character(len=pStringLen) :: petsc_optionsOrder
   class(tNode), pointer :: &
     num_mesh, &
@@ -113,7 +113,7 @@ subroutine FEM_utilities_init
   PetscErrorCode            :: ierr
 
   print'(/,a)',   ' <<<+-  FEM_utilities init  -+>>>'
- 
+
   num_mesh    => config_numerics%get('mesh',defaultVal=emptyDict)
   structOrder =  num_mesh%get_asInt('structOrder', defaultVal = 2)
 
@@ -141,7 +141,7 @@ subroutine FEM_utilities_init
   write(petsc_optionsOrder,'(a,i0)') '-mechFE_petscspace_degree ', structOrder
   call PetscOptionsInsertString(PETSC_NULL_OPTIONS,trim(petsc_optionsOrder),ierr)
   CHKERRQ(ierr)
-  
+
   wgt = 1.0/real(mesh_maxNips*mesh_NcpElemsGlobal,pReal)
 
 
@@ -152,20 +152,21 @@ end subroutine FEM_utilities_init
 !> @brief calculates constitutive response
 !--------------------------------------------------------------------------------------------------
 subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
- 
+
   real(pReal), intent(in)                 :: timeinc                                                !< loading time
   logical,     intent(in)                 :: forwardData                                            !< age results
-  
+
   real(pReal),intent(out), dimension(3,3) :: P_av                                                   !< average PK stress
-  
+
   PetscErrorCode :: ierr
 
   print'(/,a)', ' ... evaluating constitutive response ......................................'
 
-  call materialpoint_stressAndItsTangent(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems])                ! calculate P field
+  call homogenization_mechanical_response(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems])                ! calculate P field
+  if (.not. terminallyIll) &
+    call homogenization_mechanical_response2(timeinc,[1,mesh_maxNips],[1,mesh_NcpElems])
+  cutBack = .false.
 
-  cutBack = .false.                                                                                 ! reset cutBack status
-  
   P_av = sum(homogenization_P,dim=3) * wgt
   call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD,ierr)
 
@@ -198,12 +199,12 @@ subroutine utilities_projectBCValues(localVec,section,field,comp,bcPointsIS,BCVa
     do dof = offset+comp+1, offset+numDof, numComp
       localArray(dof) = localArray(dof) + BCValue + BCDotValue*timeinc
     enddo
-  enddo    
+  enddo
   call VecRestoreArrayF90(localVec,localArray,ierr); CHKERRQ(ierr)
   call VecAssemblyBegin(localVec, ierr); CHKERRQ(ierr)
   call VecAssemblyEnd  (localVec, ierr); CHKERRQ(ierr)
   if (nBcPoints > 0) call ISRestoreIndicesF90(bcPointsIS,bcPoints,ierr)
- 
+
 end subroutine utilities_projectBCValues
 
 end module FEM_utilities

src/phase.f90

@@ -108,9 +108,13 @@ module phase
       logical, intent(in) :: includeL
     end subroutine mechanical_restore
 
+    module subroutine damage_restore(ce)
+      integer, intent(in) :: ce
+    end subroutine damage_restore
 
-    module function phase_mechanical_dPdF(dt,co,ce) result(dPdF)
+
-      real(pReal), intent(in) :: dt
+    module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF)
+      real(pReal), intent(in) :: Delta_t
       integer, intent(in) :: &
         co, &                                                                                       !< counter in constituent loop
         ce
@@ -164,8 +168,8 @@ module phase
       real(pReal) :: dot_T
     end function thermal_dot_T
 
-    module function damage_phi(ph,me) result(phi)
+    module function damage_phi(ph,en) result(phi)
-      integer, intent(in) :: ph,me
+      integer, intent(in) :: ph,en
       real(pReal) :: phi
     end function damage_phi
 
@@ -209,29 +213,27 @@ module phase
 
 ! == cleaned:end ===================================================================================
 
-    module function thermal_stress(Delta_t,ph,en) result(converged_)
+    module function phase_thermal_constitutive(Delta_t,ph,en) result(converged_)
 
       real(pReal), intent(in) :: Delta_t
       integer, intent(in) :: ph, en
       logical :: converged_
 
-    end function thermal_stress
+    end function phase_thermal_constitutive
 
-    module function integrateDamageState(dt,co,ce) result(broken)
+    module function phase_damage_constitutive(Delta_t,co,ip,el) result(converged_)
-      real(pReal), intent(in) :: dt
+      real(pReal), intent(in) :: Delta_t
-      integer, intent(in) :: &
-        ce, &
-        co
-      logical :: broken
-    end function integrateDamageState
-
-    module function crystallite_stress(dt,co,ip,el) result(converged_)
-      real(pReal), intent(in) :: dt
       integer, intent(in) :: co, ip, el
       logical :: converged_
-    end function crystallite_stress
+    end function phase_damage_constitutive
 
-    !ToDo: Try to merge the all stiffness functions
+    module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged_)
+      real(pReal), intent(in) :: Delta_t
+      integer, intent(in) :: co, ip, el
+      logical :: converged_
+    end function phase_mechanical_constitutive
+
+    !ToDo: Merge all the stiffness functions
     module function phase_homogenizedC(ph,en) result(C)
       integer, intent(in) :: ph, en
       real(pReal), dimension(6,6) :: C
@@ -271,8 +273,8 @@ module phase
     end subroutine plastic_dependentState
 
 
-    module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me)
+    module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,en)
-      integer, intent(in) :: ph, me
+      integer, intent(in) :: ph, en
       real(pReal),   intent(in),  dimension(3,3) :: &
         S
       real(pReal),   intent(out), dimension(3,3) :: &
@@ -315,14 +317,14 @@ module phase
     plastic_nonlocal_updateCompatibility, &
     converged, &
     crystallite_init, &
-    crystallite_stress, &
+    phase_mechanical_constitutive, &
-    thermal_stress, &
+    phase_thermal_constitutive, &
+    phase_damage_constitutive, &
     phase_mechanical_dPdF, &
     crystallite_orientations, &
     crystallite_push33ToRef, &
     phase_restartWrite, &
     phase_restartRead, &
-    integrateDamageState, &
     phase_thermal_setField, &
     phase_set_phi, &
     phase_P, &
@@ -435,17 +437,9 @@ subroutine phase_restore(ce,includeL)
   logical, intent(in) :: includeL
   integer, intent(in) :: ce
 
-  integer :: &
-    co
-
-
-  do co = 1,homogenization_Nconstituents(material_homogenizationID(ce))
-    if (damageState(material_phaseID(co,ce))%sizeState > 0) &
-    damageState(material_phaseID(co,ce))%state( :,material_phaseEntry(co,ce)) = &
-      damageState(material_phaseID(co,ce))%state0(:,material_phaseEntry(co,ce))
-  enddo
 
   call mechanical_restore(ce,includeL)
+  call damage_restore(ce)
 
 end subroutine phase_restore
 
@@ -545,10 +539,6 @@ subroutine crystallite_init()
 
   phases => config_material%get('phase')
 
-  do ph = 1, phases%length
-    if (damageState(ph)%sizeState > 0) allocate(damageState(ph)%subState0,source=damageState(ph)%state0)  ! ToDo: hack
-  enddo
-
   print'(a42,1x,i10)', '    # of elements:                       ', eMax
   print'(a42,1x,i10)', '    # of integration points/element:     ', iMax
   print'(a42,1x,i10)', 'max # of constituents/integration point: ', cMax

src/phase_damage.f90

@@ -5,7 +5,7 @@ submodule(phase) damage
 
   type :: tDamageParameters
     real(pReal) ::                 mu = 0.0_pReal                                                   !< viscosity
-    real(pReal), dimension(3,3) :: K  = 0.0_pReal                                                   !< conductivity/diffusivity
+    real(pReal), dimension(3,3) :: D  = 0.0_pReal                                                   !< conductivity/diffusivity
   end type tDamageParameters
 
   enum, bind(c); enumerator :: &
@@ -18,7 +18,7 @@ submodule(phase) damage
 
 
   type :: tDataContainer
-    real(pReal), dimension(:), allocatable :: phi, d_phi_d_dot_phi
+    real(pReal), dimension(:), allocatable :: phi
   end type tDataContainer
 
   integer(kind(DAMAGE_UNDEFINED_ID)),     dimension(:), allocatable :: &
@@ -39,8 +39,8 @@ submodule(phase) damage
     end function isobrittle_init
 
 
-    module subroutine isobrittle_deltaState(C, Fe, ph, me)
+    module subroutine isobrittle_deltaState(C, Fe, ph, en)
-      integer, intent(in) :: ph,me
+      integer, intent(in) :: ph,en
       real(pReal),  intent(in), dimension(3,3) :: &
         Fe
       real(pReal),  intent(in), dimension(6,6) :: &
@@ -48,8 +48,8 @@ submodule(phase) damage
     end subroutine isobrittle_deltaState
 
 
-    module subroutine anisobrittle_dotState(S, ph, me)
+    module subroutine anisobrittle_dotState(S, ph, en)
-      integer, intent(in) :: ph,me
+      integer, intent(in) :: ph,en
       real(pReal),  intent(in), dimension(3,3) :: &
         S
     end subroutine anisobrittle_dotState
@@ -97,7 +97,6 @@ module subroutine damage_init
     Nmembers = count(material_phaseID == ph)
 
     allocate(current(ph)%phi(Nmembers),source=1.0_pReal)
-    allocate(current(ph)%d_phi_d_dot_phi(Nmembers),source=0.0_pReal)
 
     phase => phases%get(ph)
     sources => phase%get('damage',defaultVal=emptyList)
@@ -105,10 +104,10 @@ module subroutine damage_init
     if (sources%length == 1) then
       damage_active = .true.
       source => sources%get(1)
-      param(ph)%mu     = source%get_asFloat('mu',defaultVal=0.0_pReal)                              ! ToDo: make mandatory?
+      param(ph)%mu     = source%get_asFloat('mu')
-      param(ph)%K(1,1) = source%get_asFloat('K_11',defaultVal=0.0_pReal)                            ! ToDo: make mandatory?
+      param(ph)%D(1,1) = source%get_asFloat('D_11')
-      param(ph)%K(3,3) = source%get_asFloat('K_33',defaultVal=0.0_pReal)                            ! ToDo: depends on symmetry
+      if (any(phase_lattice(ph) == ['hP','tI'])) param(ph)%D(3,3) = source%get_asFloat('D_33')
-      param(ph)%K = lattice_applyLatticeSymmetry33(param(ph)%K,phase_lattice(ph))
+      param(ph)%D = lattice_symmetrize_33(param(ph)%D,phase_lattice(ph))
     endif
 
   enddo
@@ -125,6 +124,29 @@ module subroutine damage_init
 end subroutine damage_init
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief calculate stress (P)
+!--------------------------------------------------------------------------------------------------
+module function phase_damage_constitutive(Delta_t,co,ip,el) result(converged_)
+
+  real(pReal), intent(in) :: Delta_t
+  integer, intent(in) :: &
+    co, &
+    ip, &
+    el
+  logical :: converged_
+
+  integer :: &
+    ph, en
+
+  ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
+  en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
+
+  converged_ = .not. integrateDamageState(Delta_t,ph,en)
+
+end function phase_damage_constitutive
+
+
 !--------------------------------------------------------------------------------------------------
 !> @brief returns the degraded/modified elasticity matrix
 !--------------------------------------------------------------------------------------------------
@@ -144,6 +166,26 @@ module function phase_damage_C(C_homogenized,ph,en) result(C)
 end function phase_damage_C
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief Restore data after homog cutback.
+!--------------------------------------------------------------------------------------------------
+module subroutine damage_restore(ce)
+
+  integer, intent(in) :: ce
+
+  integer :: &
+    co
+
+
+  do co = 1,homogenization_Nconstituents(material_homogenizationID(ce))
+    if (damageState(material_phaseID(co,ce))%sizeState > 0) &
+    damageState(material_phaseID(co,ce))%state( :,material_phaseEntry(co,ce)) = &
+      damageState(material_phaseID(co,ce))%state0(:,material_phaseEntry(co,ce))
+  enddo
+
+end subroutine damage_restore
+
+
 !----------------------------------------------------------------------------------------------
 !< @brief returns local part of nonlocal damage driving force
 !----------------------------------------------------------------------------------------------
@@ -173,23 +215,20 @@ module function phase_f_phi(phi,co,ce) result(f)
 end function phase_f_phi
 
 
-
 !--------------------------------------------------------------------------------------------------
 !> @brief integrate stress, state with adaptive 1st order explicit Euler method
 !> using Fixed Point Iteration to adapt the stepsize
 !--------------------------------------------------------------------------------------------------
-module function integrateDamageState(dt,co,ce) result(broken)
+function integrateDamageState(Delta_t,ph,en) result(broken)
 
-  real(pReal), intent(in) :: dt
+  real(pReal), intent(in) :: Delta_t
   integer, intent(in) :: &
-    ce, &
+    ph, &
-    co
+    en
   logical :: broken
 
   integer :: &
     NiterationState, &                                                                              !< number of iterations in state loop
-    ph, &
-    me, &
     size_so
   real(pReal) :: &
     zeta
@@ -199,8 +238,6 @@ module function integrateDamageState(dt,co,ce) result(broken)
   logical :: &
     converged_
 
-  ph = material_phaseID(co,ce)
-  me = material_phaseEntry(co,ce)
 
   if (damageState(ph)%sizeState == 0) then
     broken = .false.
@@ -208,37 +245,37 @@ module function integrateDamageState(dt,co,ce) result(broken)
   endif
 
   converged_ = .true.
-  broken = phase_damage_collectDotState(ph,me)
+  broken = phase_damage_collectDotState(ph,en)
   if(broken) return
 
-    size_so = damageState(ph)%sizeDotState
+  size_so = damageState(ph)%sizeDotState
-    damageState(ph)%state(1:size_so,me) = damageState(ph)%subState0(1:size_so,me) &
+  damageState(ph)%state(1:size_so,en) = damageState(ph)%state0  (1:size_so,en) &
-                                        + damageState(ph)%dotState (1:size_so,me) * dt
+                                      + damageState(ph)%dotState(1:size_so,en) * Delta_t
-    source_dotState(1:size_so,2) = 0.0_pReal
+  source_dotState(1:size_so,2) = 0.0_pReal
 
   iteration: do NiterationState = 1, num%nState
 
     if(nIterationState > 1) source_dotState(1:size_so,2) = source_dotState(1:size_so,1)
-    source_dotState(1:size_so,1) = damageState(ph)%dotState(:,me)
+    source_dotState(1:size_so,1) = damageState(ph)%dotState(:,en)
 
-    broken = phase_damage_collectDotState(ph,me)
+    broken = phase_damage_collectDotState(ph,en)
     if(broken) exit iteration
 
 
-      zeta = damper(damageState(ph)%dotState(:,me),source_dotState(1:size_so,1),source_dotState(1:size_so,2))
+      zeta = damper(damageState(ph)%dotState(:,en),source_dotState(1:size_so,1),source_dotState(1:size_so,2))
-      damageState(ph)%dotState(:,me) = damageState(ph)%dotState(:,me) * zeta &
+      damageState(ph)%dotState(:,en) = damageState(ph)%dotState(:,en) * zeta &
                                      + source_dotState(1:size_so,1)* (1.0_pReal - zeta)
-      r(1:size_so) = damageState(ph)%state    (1:size_so,me)  &
+      r(1:size_so) = damageState(ph)%state   (1:size_so,en)  &
-                   - damageState(ph)%subState0(1:size_so,me)  &
+                   - damageState(ph)%State0  (1:size_so,en)  &
-                   - damageState(ph)%dotState (1:size_so,me) * dt
+                   - damageState(ph)%dotState(1:size_so,en) * Delta_t
-      damageState(ph)%state(1:size_so,me) = damageState(ph)%state(1:size_so,me) - r(1:size_so)
+      damageState(ph)%state(1:size_so,en) = damageState(ph)%state(1:size_so,en) - r(1:size_so)
       converged_ = converged_  .and. converged(r(1:size_so), &
-                                               damageState(ph)%state(1:size_so,me), &
+                                               damageState(ph)%state(1:size_so,en), &
                                                damageState(ph)%atol(1:size_so))
 
 
     if(converged_) then
-      broken = phase_damage_deltaState(mechanical_F_e(ph,me),ph,me)
+      broken = phase_damage_deltaState(mechanical_F_e(ph,en),ph,en)
       exit iteration
     endif
 
@@ -300,11 +337,11 @@ end subroutine damage_results
 !--------------------------------------------------------------------------------------------------
 !> @brief contains the constitutive equation for calculating the rate of change of microstructure
 !--------------------------------------------------------------------------------------------------
-function phase_damage_collectDotState(ph,me) result(broken)
+function phase_damage_collectDotState(ph,en) result(broken)
 
   integer, intent(in) :: &
     ph, &
-    me                                                                                         !< counter in source loop
+    en                                                                                         !< counter in source loop
   logical :: broken
 
 
@@ -315,11 +352,11 @@ function phase_damage_collectDotState(ph,me) result(broken)
     sourceType: select case (phase_damage(ph))
 
       case (DAMAGE_ANISOBRITTLE_ID) sourceType
-        call anisobrittle_dotState(mechanical_S(ph,me), ph,me) ! correct stress?
+        call anisobrittle_dotState(mechanical_S(ph,en), ph,en) ! correct stress?
 
     end select sourceType
 
-    broken = broken .or. any(IEEE_is_NaN(damageState(ph)%dotState(:,me)))
+    broken = broken .or. any(IEEE_is_NaN(damageState(ph)%dotState(:,en)))
 
   endif
 
@@ -347,9 +384,10 @@ module function phase_K_phi(co,ce) result(K)
 
   integer, intent(in) :: co, ce
   real(pReal), dimension(3,3) :: K
-
+  real(pReal), parameter :: l = 1.0_pReal
-
+  
-  K = crystallite_push33ToRef(co,ce,param(material_phaseID(co,ce))%K)
+  K = crystallite_push33ToRef(co,ce,param(material_phaseID(co,ce))%D) \
+    * l**2.0_pReal
 
 end function phase_K_phi
 
@@ -358,11 +396,11 @@ end function phase_K_phi
 !> @brief for constitutive models having an instantaneous change of state
 !> will return false if delta state is not needed/supported by the constitutive model
 !--------------------------------------------------------------------------------------------------
-function phase_damage_deltaState(Fe, ph, me) result(broken)
+function phase_damage_deltaState(Fe, ph, en) result(broken)
 
   integer, intent(in) :: &
     ph, &
-    me
+    en
   real(pReal),   intent(in), dimension(3,3) :: &
     Fe                                                                                              !< elastic deformation gradient
   integer :: &
@@ -379,13 +417,13 @@ function phase_damage_deltaState(Fe, ph, me) result(broken)
    sourceType: select case (phase_damage(ph))
 
     case (DAMAGE_ISOBRITTLE_ID) sourceType
-      call isobrittle_deltaState(phase_homogenizedC(ph,me), Fe, ph,me)
+      call isobrittle_deltaState(phase_homogenizedC(ph,en), Fe, ph,en)
-      broken = any(IEEE_is_NaN(damageState(ph)%deltaState(:,me)))
+      broken = any(IEEE_is_NaN(damageState(ph)%deltaState(:,en)))
       if(.not. broken) then
         myOffset = damageState(ph)%offsetDeltaState
         mySize   = damageState(ph)%sizeDeltaState
-        damageState(ph)%state(myOffset + 1: myOffset + mySize,me) = &
+        damageState(ph)%state(myOffset + 1: myOffset + mySize,en) = &
-        damageState(ph)%state(myOffset + 1: myOffset + mySize,me) + damageState(ph)%deltaState(1:mySize,me)
+        damageState(ph)%state(myOffset + 1: myOffset + mySize,en) + damageState(ph)%deltaState(1:mySize,en)
       endif
 
   end select sourceType
@@ -436,13 +474,13 @@ module subroutine phase_set_phi(phi,co,ce)
 end subroutine phase_set_phi
 
 
-module function damage_phi(ph,me) result(phi)
+module function damage_phi(ph,en) result(phi)
 
-  integer, intent(in) :: ph, me
+  integer, intent(in) :: ph, en
   real(pReal) :: phi
 
 
-  phi = current(ph)%phi(me)
+  phi = current(ph)%phi(en)
 
 end function damage_phi
 

src/phase_damage_anisobrittle.f90

@@ -40,7 +40,7 @@ module function anisobrittle_init() result(mySources)
     phase, &
     sources, &
     src
-  integer :: Nmembers,p
+  integer :: Nmembers,ph
   integer, dimension(:), allocatable :: N_cl
   character(len=pStringLen) :: extmsg = ''
 
@@ -56,12 +56,12 @@ module function anisobrittle_init() result(mySources)
   allocate(param(phases%length))
 
 
-  do p = 1, phases%length
+  do ph = 1, phases%length
-    if(mySources(p)) then
+    if(mySources(ph)) then
-    phase => phases%get(p)
+      phase => phases%get(ph)
-    sources => phase%get('damage')
+      sources => phase%get('damage')
 
-        associate(prm  => param(p))
+      associate(prm  => param(ph))
         src => sources%get(1)
 
         N_cl = src%get_as1dInt('N_cl',defaultVal=emptyIntArray)
@@ -73,8 +73,7 @@ module function anisobrittle_init() result(mySources)
         prm%s_crit  = src%get_as1dFloat('s_crit',  requiredSize=size(N_cl))
         prm%g_crit  = src%get_as1dFloat('g_crit',  requiredSize=size(N_cl))
 
-        prm%cleavage_systems = lattice_SchmidMatrix_cleavage(N_cl,phase%get_asString('lattice'),&
+        prm%cleavage_systems = lattice_SchmidMatrix_cleavage(N_cl,phase_lattice(ph),phase_cOverA(ph))
-                                                             phase%get_asFloat('c/a',defaultVal=0.0_pReal))
 
         ! expand: family => system
         prm%s_crit = math_expand(prm%s_crit,N_cl)
@@ -92,17 +91,15 @@ module function anisobrittle_init() result(mySources)
         if (any(prm%g_crit <  0.0_pReal)) extmsg = trim(extmsg)//' g_crit'
         if (any(prm%s_crit <  0.0_pReal)) extmsg = trim(extmsg)//' s_crit'
 
-        Nmembers = count(material_phaseID==p)
+        Nmembers = count(material_phaseID==ph)
-        call phase_allocateState(damageState(p),Nmembers,1,1,0)
+        call phase_allocateState(damageState(ph),Nmembers,1,1,0)
-        damageState(p)%atol = src%get_asFloat('atol_phi',defaultVal=1.0e-9_pReal)
+        damageState(ph)%atol = src%get_asFloat('atol_phi',defaultVal=1.0e-9_pReal)
-        if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' atol_phi'
+        if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' atol_phi'
 
-        end associate
+      end associate
 
-!--------------------------------------------------------------------------------------------------
+      if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_anisoBrittle)')
-!  exit if any parameter is out of range
+    endif
-        if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_anisoBrittle)')
-      endif
 
   enddo
 
@@ -112,10 +109,10 @@ end function anisobrittle_init
 !--------------------------------------------------------------------------------------------------
 !> @brief
 !--------------------------------------------------------------------------------------------------
-module subroutine anisobrittle_dotState(S, ph,me)
+module subroutine anisobrittle_dotState(S, ph,en)
 
   integer, intent(in) :: &
-    ph,me
+    ph,en
   real(pReal),  intent(in), dimension(3,3) :: &
     S
 
@@ -126,15 +123,15 @@ module subroutine anisobrittle_dotState(S, ph,me)
 
 
   associate(prm => param(ph))
-    damageState(ph)%dotState(1,me) = 0.0_pReal
+    damageState(ph)%dotState(1,en) = 0.0_pReal
     do i = 1, prm%sum_N_cl
       traction_d = math_tensordot(S,prm%cleavage_systems(1:3,1:3,1,i))
       traction_t = math_tensordot(S,prm%cleavage_systems(1:3,1:3,2,i))
       traction_n = math_tensordot(S,prm%cleavage_systems(1:3,1:3,3,i))
 
-      traction_crit = prm%g_crit(i)*damage_phi(ph,me)**2.0_pReal
+      traction_crit = prm%g_crit(i)*damage_phi(ph,en)**2.0_pReal
 
-      damageState(ph)%dotState(1,me) = damageState(ph)%dotState(1,me) &
+      damageState(ph)%dotState(1,en) = damageState(ph)%dotState(1,en) &
           + prm%dot_o / prm%s_crit(i) &
             * ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**prm%q + &
                (max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**prm%q + &
@@ -171,10 +168,10 @@ end subroutine anisobrittle_results
 !--------------------------------------------------------------------------------------------------
 !> @brief  contains the constitutive equation for calculating the velocity gradient
 !--------------------------------------------------------------------------------------------------
-module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me)
+module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,en)
 
   integer, intent(in) :: &
-    ph,me
+    ph,en
   real(pReal),   intent(in),  dimension(3,3) :: &
     S
   real(pReal),   intent(out), dimension(3,3) :: &
@@ -193,7 +190,7 @@ module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me)
   dLd_dTstar = 0.0_pReal
   associate(prm => param(ph))
   do i = 1,prm%sum_N_cl
-    traction_crit = prm%g_crit(i)*damage_phi(ph,me)**2.0_pReal
+    traction_crit = prm%g_crit(i)*damage_phi(ph,en)**2.0_pReal
 
     traction_d = math_tensordot(S,prm%cleavage_systems(1:3,1:3,1,i))
     if (abs(traction_d) > traction_crit + tol_math_check) then

src/phase_damage_isobrittle.f90

@@ -13,7 +13,15 @@ submodule(phase:damage) isobrittle
       output
   end type tParameters
 
-  type(tParameters), dimension(:), allocatable :: param                                             !< containers of constitutive parameters (len Ninstances)
+  type :: tIsobrittleState
+    real(pReal), pointer, dimension(:) :: &                                                         !< vectors along Nmembers
+      r_W                                                                                           !< ratio between actual and critical strain energy density
+  end type tIsobrittleState
+
+  type(tParameters),      allocatable, dimension(:) :: param                                        !< containers of constitutive parameters (len Ninstances)
+  type(tIsobrittleState), allocatable, dimension(:) :: &
+    deltaState, &
+    state
 
 contains
 
@@ -44,13 +52,15 @@ module function isobrittle_init() result(mySources)
 
   phases => config_material%get('phase')
   allocate(param(phases%length))
+  allocate(state(phases%length))
+  allocate(deltaState(phases%length))
 
   do ph = 1, phases%length
     if(mySources(ph)) then
-    phase => phases%get(ph)
+      phase => phases%get(ph)
-    sources => phase%get('damage')
+      sources => phase%get('damage')
 
-        associate(prm  => param(ph))
+      associate(prm => param(ph), dlt => deltaState(ph), stt => state(ph))
         src => sources%get(1)
 
         prm%W_crit = src%get_asFloat('W_crit')
@@ -69,12 +79,14 @@ module function isobrittle_init() result(mySources)
         damageState(ph)%atol = src%get_asFloat('atol_phi',defaultVal=1.0e-9_pReal)
         if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' atol_phi'
 
-        end associate
+        stt%r_W => damageState(ph)%state(1,:)
+        dlt%r_W => damageState(ph)%deltaState(1,:)
 
-!--------------------------------------------------------------------------------------------------
+      end associate
-!  exit if any parameter is out of range
+
-        if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_isobrittle)')
+
-      endif
+      if (extmsg /= '') call IO_error(211,ext_msg=trim(extmsg)//'(damage_isobrittle)')
+    endif
 
   enddo
 
@@ -85,29 +97,27 @@ end function isobrittle_init
 !--------------------------------------------------------------------------------------------------
 !> @brief calculates derived quantities from state
 !--------------------------------------------------------------------------------------------------
-module subroutine isobrittle_deltaState(C, Fe, ph,me)
+module subroutine isobrittle_deltaState(C, Fe, ph,en)
 
-  integer, intent(in) :: ph,me
+  integer, intent(in) :: ph,en
   real(pReal),  intent(in), dimension(3,3) :: &
     Fe
   real(pReal),  intent(in), dimension(6,6) :: &
     C
 
   real(pReal), dimension(6) :: &
-    strain
+    epsilon
   real(pReal) :: &
-    strainenergy
+    r_W
 
 
-  strain = 0.5_pReal*math_sym33to6(matmul(transpose(Fe),Fe)-math_I3)
+  epsilon = 0.5_pReal*math_sym33to6(matmul(transpose(Fe),Fe)-math_I3)
 
-  associate(prm => param(ph))
+  associate(prm => param(ph), stt => state(ph), dlt => deltaState(ph))
-    strainenergy = 2.0_pReal*sum(strain*matmul(C,strain))/prm%W_crit
+
-    ! ToDo: check strainenergy = 2.0_pReal*dot_product(strain,matmul(C,strain))/prm%W_crit
+    r_W = 2.0_pReal*dot_product(epsilon,matmul(C,epsilon))/prm%W_crit
+    dlt%r_W(en) = merge(r_W - stt%r_W(en), 0.0_pReal, r_W > stt%r_W(en))
 
-    damageState(ph)%deltaState(1,me) = merge(strainenergy - damageState(ph)%state(1,me), &
-                                             damageState(ph)%subState0(1,me) - damageState(ph)%state(1,me), &
-                                             strainenergy > damageState(ph)%subState0(1,me))
   end associate
 
 end subroutine isobrittle_deltaState
@@ -124,14 +134,15 @@ module subroutine isobrittle_results(phase,group)
   integer :: o
 
 
-  associate(prm => param(phase), &
+  associate(prm => param(phase), stt => damageState(phase)%state) ! point to state and output r_W (is scalar, not 1D vector)
-            stt => damageState(phase)%state)
+
     outputsLoop: do o = 1,size(prm%output)
       select case(trim(prm%output(o)))
         case ('f_phi')
-          call results_writeDataset(stt,group,trim(prm%output(o)),'driving force','J/m³')
+          call results_writeDataset(stt,group,trim(prm%output(o)),'driving force','J/m³') ! Wrong, this is dimensionless
       end select
     enddo outputsLoop
+
   end associate
 
 end subroutine isobrittle_results

src/phase_mechanical.f90

@@ -44,9 +44,9 @@ submodule(phase) mechanical
 
   interface
 
-    module subroutine eigendeformation_init(phases)
+    module subroutine eigen_init(phases)
       class(tNode), pointer :: phases
-    end subroutine eigendeformation_init
+    end subroutine eigen_init
 
     module subroutine elastic_init(phases)
       class(tNode), pointer :: phases
@@ -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)
@@ -197,8 +211,7 @@ module subroutine mechanical_init(materials,phases)
     Nmembers
   class(tNode), pointer :: &
     num_crystallite, &
-    material, &
+
-    constituents, &
     phase, &
     mech
 
@@ -303,7 +316,7 @@ module subroutine mechanical_init(materials,phases)
   end select
 
 
-  call eigendeformation_init(phases)
+  call eigen_init(phases)
 
 
 end subroutine mechanical_init
@@ -977,9 +990,9 @@ end subroutine mechanical_forward
 !--------------------------------------------------------------------------------------------------
 !> @brief calculate stress (P)
 !--------------------------------------------------------------------------------------------------
-module function crystallite_stress(dt,co,ip,el) result(converged_)
+module function phase_mechanical_constitutive(Delta_t,co,ip,el) result(converged_)
 
-  real(pReal), intent(in) :: dt
+  real(pReal), intent(in) :: Delta_t
   integer, intent(in) :: &
     co, &
     ip, &
@@ -1009,17 +1022,13 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
   subLi0 = phase_mechanical_Li0(ph)%data(1:3,1:3,en)
   subLp0 = phase_mechanical_Lp0(ph)%data(1:3,1:3,en)
   subState0 = plasticState(ph)%State0(:,en)
-
-  if (damageState(ph)%sizeState > 0) &
-    damageState(ph)%subState0(:,en) = damageState(ph)%state0(:,en)
-
   subFp0 = phase_mechanical_Fp0(ph)%data(1:3,1:3,en)
   subFi0 = phase_mechanical_Fi0(ph)%data(1:3,1:3,en)
   subF0  = phase_mechanical_F0(ph)%data(1:3,1:3,en)
   subFrac = 0.0_pReal
-  subStep = 1.0_pReal/num%subStepSizeCryst
   todo = .true.
-  converged_ = .false.                                                      ! pretend failed step of 1/subStepSizeCryst
+  subStep = 1.0_pReal/num%subStepSizeCryst
+  converged_ = .false.                                                                              ! pretend failed step of 1/subStepSizeCryst
 
   todo = .true.
   cutbackLooping: do while (todo)
@@ -1038,9 +1047,6 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
         subFp0 = phase_mechanical_Fp(ph)%data(1:3,1:3,en)
         subFi0 = phase_mechanical_Fi(ph)%data(1:3,1:3,en)
         subState0 = plasticState(ph)%state(:,en)
-        if (damageState(ph)%sizeState > 0) &
-          damageState(ph)%subState0(:,en) = damageState(ph)%state(:,en)
-
       endif
 !--------------------------------------------------------------------------------------------------
 !  cut back (reduced time and restore)
@@ -1048,16 +1054,13 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
       subStep       = num%subStepSizeCryst * subStep
       phase_mechanical_Fp(ph)%data(1:3,1:3,en) = subFp0
       phase_mechanical_Fi(ph)%data(1:3,1:3,en) = subFi0
-      phase_mechanical_S(ph)%data(1:3,1:3,en) = phase_mechanical_S0(ph)%data(1:3,1:3,en)          ! why no subS0 ? is S0 of any use?
+      phase_mechanical_S(ph)%data(1:3,1:3,en) = phase_mechanical_S0(ph)%data(1:3,1:3,en)            ! why no subS0 ? is S0 of any use?
       if (subStep < 1.0_pReal) then                                                                 ! actual (not initial) cutback
         phase_mechanical_Lp(ph)%data(1:3,1:3,en) = subLp0
         phase_mechanical_Li(ph)%data(1:3,1:3,en) = subLi0
       endif
       plasticState(ph)%state(:,en) = subState0
-      if (damageState(ph)%sizeState > 0) &
+      todo = subStep > num%subStepMinCryst                                                          ! still on track or already done (beyond repair)
-        damageState(ph)%state(:,en) = damageState(ph)%subState0(:,en)
-
-      todo = subStep > num%subStepMinCryst                          ! still on track or already done (beyond repair)
     endif
 
 !--------------------------------------------------------------------------------------------------
@@ -1065,13 +1068,12 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
     if (todo) then
       subF = subF0 &
            + subStep * (phase_mechanical_F(ph)%data(1:3,1:3,en) - phase_mechanical_F0(ph)%data(1:3,1:3,en))
-      converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * dt,co,ip,el)
+      converged_ = .not. integrateState(subF0,subF,subFp0,subFi0,subState0(1:sizeDotState),subStep * Delta_t,co,ip,el)
-      converged_ = converged_ .and. .not. integrateDamageState(subStep * dt,co,(el-1)*discretization_nIPs + ip)
     endif
 
   enddo cutbackLooping
 
-end function crystallite_stress
+end function phase_mechanical_constitutive
 
 
 !--------------------------------------------------------------------------------------------------
@@ -1104,12 +1106,13 @@ module subroutine mechanical_restore(ce,includeL)
 
 end subroutine mechanical_restore
 
+
 !--------------------------------------------------------------------------------------------------
 !> @brief Calculate tangent (dPdF).
 !--------------------------------------------------------------------------------------------------
-module function phase_mechanical_dPdF(dt,co,ce) result(dPdF)
+module function phase_mechanical_dPdF(Delta_t,co,ce) result(dPdF)
 
-  real(pReal), intent(in) :: dt
+  real(pReal), intent(in) :: Delta_t
   integer, intent(in) :: &
     co, &                                                                                            !< counter in constituent loop
     ce
@@ -1159,11 +1162,11 @@ module function phase_mechanical_dPdF(dt,co,ce) result(dPdF)
     lhs_3333 = 0.0_pReal; rhs_3333 = 0.0_pReal
     do o=1,3; do p=1,3
       lhs_3333(1:3,1:3,o,p) = lhs_3333(1:3,1:3,o,p) &
-                            + matmul(invSubFi0,dLidFi(1:3,1:3,o,p)) * dt
+                            + matmul(invSubFi0,dLidFi(1:3,1:3,o,p)) * Delta_t
       lhs_3333(1:3,o,1:3,p) = lhs_3333(1:3,o,1:3,p) &
                             + invFi*invFi(p,o)
       rhs_3333(1:3,1:3,o,p) = rhs_3333(1:3,1:3,o,p) &
-                            - matmul(invSubFi0,dLidS(1:3,1:3,o,p)) * dt
+                            - matmul(invSubFi0,dLidS(1:3,1:3,o,p)) * Delta_t
     enddo; enddo
     call math_invert(temp_99,error,math_3333to99(lhs_3333))
     if (error) then
@@ -1192,7 +1195,7 @@ module function phase_mechanical_dPdF(dt,co,ce) result(dPdF)
     temp_3333(1:3,1:3,p,o) = matmul(matmul(temp_33_2,dLpdS(1:3,1:3,p,o)), invFi) &
                            + matmul(temp_33_3,dLidS(1:3,1:3,p,o))
   enddo; enddo
-  lhs_3333 = math_mul3333xx3333(dSdFe,temp_3333) * dt &
+  lhs_3333 = math_mul3333xx3333(dSdFe,temp_3333) * Delta_t &
            + math_mul3333xx3333(dSdFi,dFidS)
 
   call math_invert(temp_99,error,math_eye(9)+math_3333to99(lhs_3333))
@@ -1208,7 +1211,7 @@ module function phase_mechanical_dPdF(dt,co,ce) result(dPdF)
 ! calculate dFpinvdF
   temp_3333 = math_mul3333xx3333(dLpdS,dSdF)
   do o=1,3; do p=1,3
-    dFpinvdF(1:3,1:3,p,o) = - matmul(invSubFp0, matmul(temp_3333(1:3,1:3,p,o),invFi)) * dt
+    dFpinvdF(1:3,1:3,p,o) = - matmul(invSubFp0, matmul(temp_3333(1:3,1:3,p,o),invFi)) * Delta_t
   enddo; enddo
 
 !--------------------------------------------------------------------------------------------------

src/phase_mechanical_eigen.f90

@@ -32,7 +32,7 @@ submodule(phase:mechanical) eigen
 contains
 
 
-module subroutine eigendeformation_init(phases)
+module subroutine eigen_init(phases)
 
   class(tNode), pointer :: &
     phases
@@ -68,7 +68,7 @@ module subroutine eigendeformation_init(phases)
   where(damage_anisobrittle_init())  model_damage = KINEMATICS_cleavage_opening_ID
 
 
-end subroutine eigendeformation_init
+end subroutine eigen_init
 
 
 !--------------------------------------------------------------------------------------------------
@@ -86,17 +86,17 @@ function kinematics_active(kinematics_label,kinematics_length)  result(active_ki
     kinematics, &
     kinematics_type, &
     mechanics
-  integer :: p,k
+  integer :: ph,k
 
   phases => config_material%get('phase')
   allocate(active_kinematics(kinematics_length,phases%length), source = .false. )
-  do p = 1, phases%length
+  do ph = 1, phases%length
-    phase => phases%get(p)
+    phase => phases%get(ph)
     mechanics => phase%get('mechanical')
     kinematics => mechanics%get('eigen',defaultVal=emptyList)
     do k = 1, kinematics%length
       kinematics_type => kinematics%get(k)
-      active_kinematics(k,p) = kinematics_type%get_asString('type') == kinematics_label
+      active_kinematics(k,ph) = kinematics_type%get_asString('type') == kinematics_label
     enddo
   enddo
 
@@ -118,17 +118,17 @@ function kinematics_active2(kinematics_label)  result(active_kinematics)
     phase, &
     kinematics, &
     kinematics_type
-  integer :: p
+  integer :: ph
 
   phases => config_material%get('phase')
-  allocate(active_kinematics(phases%length), source = .false. )
+  allocate(active_kinematics(phases%length), source = .false.)
-  do p = 1, phases%length
+  do ph = 1, phases%length
-    phase => phases%get(p)
+    phase => phases%get(ph)
     kinematics => phase%get('damage',defaultVal=emptyList)
     if(kinematics%length < 1) return
     kinematics_type => kinematics%get(1)
     if (.not. kinematics_type%contains('type')) continue
-    active_kinematics(p) = kinematics_type%get_asString('type',defaultVal='n/a') == kinematics_label
+    active_kinematics(ph) = kinematics_type%get_asString('type',defaultVal='n/a') == kinematics_label
   enddo
 
 

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

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 = 0.0_pReal                                                                               !< 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,14 @@ 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')
+      
+      prm%C66 = math_sym3333to66(math_Voigt66to3333(prm%C66))                                       ! Literature data is in Voigt notation
+
     end associate
   enddo
 
@@ -104,10 +116,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

src/phase_mechanical_plastic_dislotungsten.f90

@@ -36,8 +36,8 @@ submodule(phase:plastic) dislotungsten
       forestProjection
     real(pReal),               allocatable, dimension(:,:,:) :: &
       P_sl, &
-      nonSchmid_pos, &
+      P_nS_pos, &
-      nonSchmid_neg
+      P_nS_neg
     integer :: &
       sum_N_sl                                                                                      !< total number of active slip system
     character(len=pStringLen), allocatable, dimension(:) :: &
@@ -129,30 +129,28 @@ 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 = elastic_mu(ph)
-    prm%mu = lattice_mu(ph)
 
 !--------------------------------------------------------------------------------------------------
 ! slip related parameters
     N_sl         = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray)
     prm%sum_N_sl = sum(abs(N_sl))
     slipActive: if (prm%sum_N_sl > 0) then
-      prm%P_sl = lattice_SchmidMatrix_slip(N_sl,phase%get_asString('lattice'),&
+      prm%P_sl = lattice_SchmidMatrix_slip(N_sl,phase_lattice(ph),phase_cOverA(ph))
-                                           phase%get_asFloat('c/a',defaultVal=0.0_pReal))
 
-      if(trim(phase%get_asString('lattice')) == 'cI') then
+      if (phase_lattice(ph) == 'cI') then
         a = pl%get_as1dFloat('a_nonSchmid',defaultVal = emptyRealArray)
-        prm%nonSchmid_pos = lattice_nonSchmidMatrix(N_sl,a,+1)
+        prm%P_nS_pos = lattice_nonSchmidMatrix(N_sl,a,+1)
-        prm%nonSchmid_neg = lattice_nonSchmidMatrix(N_sl,a,-1)
+        prm%P_nS_neg = lattice_nonSchmidMatrix(N_sl,a,-1)
       else
-        prm%nonSchmid_pos = prm%P_sl
+        prm%P_nS_pos = prm%P_sl
-        prm%nonSchmid_neg = prm%P_sl
+        prm%P_nS_neg = prm%P_sl
       endif
 
       prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dFloat('h_sl-sl'), &
-                                                   phase%get_asString('lattice'))
+                                                   phase_lattice(ph))
-      prm%forestProjection = lattice_forestProjection_edge(N_sl,phase%get_asString('lattice'),&
+      prm%forestProjection = lattice_forestProjection_edge(N_sl,phase_lattice(ph),&
-                                                           phase%get_asFloat('c/a',defaultVal=0.0_pReal))
+                                                           phase_cOverA(ph))
       prm%forestProjection = transpose(prm%forestProjection)
 
       rho_mob_0       = pl%get_as1dFloat('rho_mob_0',     requiredSize=size(N_sl))
@@ -163,39 +161,37 @@ module function plastic_dislotungsten_init() result(myPlasticity)
 
       prm%i_sl        = pl%get_as1dFloat('i_sl',          requiredSize=size(N_sl))
       prm%tau_Peierls = pl%get_as1dFloat('tau_Peierls',   requiredSize=size(N_sl))
-      prm%p           = pl%get_as1dFloat('p_sl',          requiredSize=size(N_sl), &
+      prm%p           = pl%get_as1dFloat('p_sl',          requiredSize=size(N_sl))
-                                         defaultVal=[(1.0_pReal,i=1,size(N_sl))])
+      prm%q           = pl%get_as1dFloat('q_sl',          requiredSize=size(N_sl))
-      prm%q           = pl%get_as1dFloat('q_sl',          requiredSize=size(N_sl), &
-                                         defaultVal=[(1.0_pReal,i=1,size(N_sl))])
       prm%h           = pl%get_as1dFloat('h',             requiredSize=size(N_sl))
       prm%w           = pl%get_as1dFloat('w',             requiredSize=size(N_sl))
       prm%omega       = pl%get_as1dFloat('omega',         requiredSize=size(N_sl))
       prm%B           = pl%get_as1dFloat('B',             requiredSize=size(N_sl))
 
-      prm%D               = pl%get_asFloat('D')
+      prm%D    = pl%get_asFloat('D')
-      prm%D_0             = pl%get_asFloat('D_0')
+      prm%D_0  = pl%get_asFloat('D_0')
-      prm%Q_cl            = pl%get_asFloat('Q_cl')
+      prm%Q_cl = pl%get_asFloat('Q_cl')
-      prm%f_at            = pl%get_asFloat('f_at')       * prm%b_sl**3.0_pReal
+      prm%f_at = pl%get_asFloat('f_at') * prm%b_sl**3.0_pReal
-      prm%D_a             = pl%get_asFloat('D_a')        * prm%b_sl
+      prm%D_a  = pl%get_asFloat('D_a') * prm%b_sl
 
       prm%dipoleformation = .not. pl%get_asBool('no_dipole_formation', defaultVal = .false.)
 
       ! expand: family => system
-      rho_mob_0          = math_expand(rho_mob_0,          N_sl)
+      rho_mob_0          = math_expand(rho_mob_0,       N_sl)
-      rho_dip_0          = math_expand(rho_dip_0,          N_sl)
+      rho_dip_0          = math_expand(rho_dip_0,       N_sl)
-      prm%q              = math_expand(prm%q,              N_sl)
+      prm%q              = math_expand(prm%q,           N_sl)
-      prm%p              = math_expand(prm%p,              N_sl)
+      prm%p              = math_expand(prm%p,           N_sl)
-      prm%Q_s            = math_expand(prm%Q_s,            N_sl)
+      prm%Q_s            = math_expand(prm%Q_s,         N_sl)
-      prm%b_sl           = math_expand(prm%b_sl,           N_sl)
+      prm%b_sl           = math_expand(prm%b_sl,        N_sl)
-      prm%h              = math_expand(prm%h,              N_sl)
+      prm%h              = math_expand(prm%h,           N_sl)
-      prm%w              = math_expand(prm%w,              N_sl)
+      prm%w              = math_expand(prm%w,           N_sl)
-      prm%omega          = math_expand(prm%omega,          N_sl)
+      prm%omega          = math_expand(prm%omega,       N_sl)
-      prm%tau_Peierls    = math_expand(prm%tau_Peierls,    N_sl)
+      prm%tau_Peierls    = math_expand(prm%tau_Peierls, N_sl)
-      prm%v_0            = math_expand(prm%v_0,            N_sl)
+      prm%v_0            = math_expand(prm%v_0,         N_sl)
-      prm%B              = math_expand(prm%B,              N_sl)
+      prm%B              = math_expand(prm%B,           N_sl)
-      prm%i_sl           = math_expand(prm%i_sl,           N_sl)
+      prm%i_sl           = math_expand(prm%i_sl,        N_sl)
-      prm%f_at           = math_expand(prm%f_at,           N_sl)
+      prm%f_at           = math_expand(prm%f_at,        N_sl)
-      prm%D_a            = math_expand(prm%D_a,            N_sl)
+      prm%D_a            = math_expand(prm%D_a,         N_sl)
 
       ! sanity checks
       if (    prm%D_0          <= 0.0_pReal)  extmsg = trim(extmsg)//' D_0'
@@ -298,8 +294,8 @@ pure module subroutine dislotungsten_LpAndItsTangent(Lp,dLp_dMp, &
     Lp = Lp + (dot_gamma_pos(i)+dot_gamma_neg(i))*prm%P_sl(1:3,1:3,i)
     forall (k=1:3,l=1:3,m=1:3,n=1:3) &
       dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
-                       + ddot_gamma_dtau_pos(i) * prm%P_sl(k,l,i) * prm%nonSchmid_pos(m,n,i) &
+                       + ddot_gamma_dtau_pos(i) * prm%P_sl(k,l,i) * prm%P_nS_pos(m,n,i) &
-                       + ddot_gamma_dtau_neg(i) * prm%P_sl(k,l,i) * prm%nonSchmid_neg(m,n,i)
+                       + ddot_gamma_dtau_neg(i) * prm%P_sl(k,l,i) * prm%P_nS_neg(m,n,i)
   enddo
 
   end associate
@@ -323,7 +319,7 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en)
   real(pReal) :: &
     VacancyDiffusion
   real(pReal), dimension(param(ph)%sum_N_sl) :: &
-    gdot_pos, gdot_neg,&
+    dot_gamma_pos, dot_gamma_neg,&
     tau_pos,&
     tau_neg, &
     v_cl, &
@@ -335,10 +331,10 @@ module subroutine dislotungsten_dotState(Mp,T,ph,en)
             dot => dotState(ph), dst => dependentState(ph))
 
   call kinetics(Mp,T,ph,en,&
-                gdot_pos,gdot_neg, &
+                dot_gamma_pos,dot_gamma_neg, &
                 tau_pos_out = tau_pos,tau_neg_out = tau_neg)
 
-  dot%gamma_sl(:,en) = (gdot_pos+gdot_neg)                                                          ! ToDo: needs to be abs
+  dot%gamma_sl(:,en) = (dot_gamma_pos+dot_gamma_neg)                                                          ! ToDo: needs to be abs
   VacancyDiffusion = prm%D_0*exp(-prm%Q_cl/(kB*T))
 
   where(dEq0(tau_pos))                                                                              ! ToDo: use avg of pos and neg
@@ -380,13 +376,14 @@ module subroutine dislotungsten_dependentState(ph,en)
   real(pReal), dimension(param(ph)%sum_N_sl) :: &
     dislocationSpacing
 
+
   associate(prm => param(ph), stt => state(ph),dst => dependentState(ph))
 
-  dislocationSpacing = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,en)+stt%rho_dip(:,en)))
+    dislocationSpacing = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,en)+stt%rho_dip(:,en)))
-  dst%threshold_stress(:,en) = prm%mu*prm%b_sl &
+    dst%threshold_stress(:,en) = prm%mu*prm%b_sl &
-                             * sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en)))
+                               * sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en)))
 
-  dst%Lambda_sl(:,en) = prm%D/(1.0_pReal+prm%D*dislocationSpacing/prm%i_sl)
+    dst%Lambda_sl(:,en) = prm%D/(1.0_pReal+prm%D*dislocationSpacing/prm%i_sl)
 
   end associate
 
@@ -466,8 +463,8 @@ pure subroutine kinetics(Mp,T,ph,en, &
   associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
 
   do j = 1, prm%sum_N_sl
-    tau_pos(j) = math_tensordot(Mp,prm%nonSchmid_pos(1:3,1:3,j))
+    tau_pos(j) = math_tensordot(Mp,prm%P_nS_pos(1:3,1:3,j))
-    tau_neg(j) = math_tensordot(Mp,prm%nonSchmid_neg(1:3,1:3,j))
+    tau_neg(j) = math_tensordot(Mp,prm%P_nS_neg(1:3,1:3,j))
   enddo
 
 

src/phase_mechanical_plastic_dislotwin.f90

@@ -184,27 +184,23 @@ module function plastic_dislotwin_init() result(myPlasticity)
 #endif
 
     ! This data is read in already in lattice
-    prm%mu  = lattice_mu(ph)
+    prm%mu  = elastic_mu(ph)
-    prm%nu  = lattice_nu(ph)
+    prm%nu  = elastic_nu(ph)
-    prm%C66 = lattice_C66(1:6,1:6,ph)
+    prm%C66 = elastic_C66(ph)
 
 !--------------------------------------------------------------------------------------------------
 ! slip related parameters
     N_sl         = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray)
     prm%sum_N_sl = sum(abs(N_sl))
     slipActive: if (prm%sum_N_sl > 0) then
-      prm%P_sl    = lattice_SchmidMatrix_slip(N_sl,phase%get_asString('lattice'),&
+      prm%P_sl    = lattice_SchmidMatrix_slip(N_sl,phase_lattice(ph),phase_cOverA(ph))
-                                              phase%get_asFloat('c/a',defaultVal=0.0_pReal))
+      prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dFloat('h_sl-sl'),phase_lattice(ph))
-      prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dFloat('h_sl-sl'), &
+      prm%forestProjection = lattice_forestProjection_edge(N_sl,phase_lattice(ph),phase_cOverA(ph))
-                                                   phase%get_asString('lattice'))
-      prm%forestProjection = lattice_forestProjection_edge(N_sl,phase%get_asString('lattice'),&
-                                                           phase%get_asFloat('c/a',defaultVal=0.0_pReal))
       prm%forestProjection = transpose(prm%forestProjection)
 
-      prm%n0_sl            = lattice_slip_normal(N_sl,phase%get_asString('lattice'),&
+      prm%n0_sl = lattice_slip_normal(N_sl,phase_lattice(ph),phase_cOverA(ph))
-                                                 phase%get_asFloat('c/a',defaultVal=0.0_pReal))
       prm%fccTwinTransNucleation = phase_lattice(ph) == 'cF' .and. (N_sl(1) == 12)
-      if(prm%fccTwinTransNucleation) prm%fcc_twinNucleationSlipPair = lattice_FCC_TWINNUCLEATIONSLIPPAIR
+      if (prm%fccTwinTransNucleation) prm%fcc_twinNucleationSlipPair = lattice_FCC_TWINNUCLEATIONSLIPPAIR
 
       rho_mob_0                = pl%get_as1dFloat('rho_mob_0',   requiredSize=size(N_sl))
       rho_dip_0                = pl%get_as1dFloat('rho_dip_0',   requiredSize=size(N_sl))
@@ -265,11 +261,9 @@ module function plastic_dislotwin_init() result(myPlasticity)
     N_tw         = pl%get_as1dInt('N_tw', defaultVal=emptyIntArray)
     prm%sum_N_tw = sum(abs(N_tw))
     twinActive: if (prm%sum_N_tw > 0) then
-      prm%P_tw  = lattice_SchmidMatrix_twin(N_tw,phase%get_asString('lattice'),&
+      prm%P_tw  = lattice_SchmidMatrix_twin(N_tw,phase_lattice(ph),phase_cOverA(ph))
-                                                   phase%get_asFloat('c/a',defaultVal=0.0_pReal))
+      prm%h_tw_tw   = lattice_interaction_TwinByTwin(N_tw,pl%get_as1dFloat('h_tw-tw'), &
-      prm%h_tw_tw   = lattice_interaction_TwinByTwin(N_tw,&
+                                                     phase_lattice(ph))
-                                                     pl%get_as1dFloat('h_tw-tw'), &
-                                                     phase%get_asString('lattice'))
 
       prm%b_tw      = pl%get_as1dFloat('b_tw',     requiredSize=size(N_tw))
       prm%t_tw      = pl%get_as1dFloat('t_tw',     requiredSize=size(N_tw))
@@ -279,11 +273,9 @@ module function plastic_dislotwin_init() result(myPlasticity)
       prm%L_tw      = pl%get_asFloat('L_tw')
       prm%i_tw      = pl%get_asFloat('i_tw')
 
-      prm%gamma_char= lattice_characteristicShear_Twin(N_tw,phase%get_asString('lattice'),&
+      prm%gamma_char= lattice_characteristicShear_Twin(N_tw,phase_lattice(ph),phase_cOverA(ph))
-                                                       phase%get_asFloat('c/a',defaultVal=0.0_pReal))
 
-      prm%C66_tw    = lattice_C66_twin(N_tw,prm%C66,phase%get_asString('lattice'),&
+      prm%C66_tw    = lattice_C66_twin(N_tw,prm%C66,phase_lattice(ph),phase_cOverA(ph))
-                                       phase%get_asFloat('c/a',defaultVal=0.0_pReal))
 
       if (.not. prm%fccTwinTransNucleation) then
         prm%dot_N_0_tw = pl%get_as1dFloat('dot_N_0_tw')
@@ -324,8 +316,8 @@ module function plastic_dislotwin_init() result(myPlasticity)
       prm%x_c_tr        = pl%get_asFloat('x_c_tr',  defaultVal=0.0_pReal) ! ToDo: How to handle that???
       prm%L_tr          = pl%get_asFloat('L_tr')
 
-      prm%h_tr_tr = lattice_interaction_TransByTrans(N_tr,pl%get_as1dFloat('h_tr-tr'), &
+      prm%h_tr_tr = lattice_interaction_TransByTrans(N_tr,pl%get_as1dFloat('h_tr-tr'),&
-                                                     phase%get_asString('lattice'))
+                                                     phase_lattice(ph))
 
       prm%C66_tr  = lattice_C66_trans(N_tr,prm%C66,pl%get_asString('lattice_tr'), &
                                       0.0_pReal, &
@@ -391,17 +383,15 @@ module function plastic_dislotwin_init() result(myPlasticity)
     endif
 
     slipAndTwinActive: if (prm%sum_N_sl * prm%sum_N_tw > 0) then
-      prm%h_sl_tw = lattice_interaction_SlipByTwin(N_sl,N_tw,&
+      prm%h_sl_tw = lattice_interaction_SlipByTwin(N_sl,N_tw,pl%get_as1dFloat('h_sl-tw'), &
-                                                   pl%get_as1dFloat('h_sl-tw'), &
+                                                   phase_lattice(ph))
-                                                   phase%get_asString('lattice'))
+      if (prm%fccTwinTransNucleation .and. size(N_tw) /= 1) extmsg = trim(extmsg)//' N_tw: nucleation'
-      if (prm%fccTwinTransNucleation .and. size(N_tw) /= 1) extmsg = trim(extmsg)//' interaction_sliptwin'
     endif slipAndTwinActive
 
     slipAndTransActive: if (prm%sum_N_sl * prm%sum_N_tr > 0) then
-      prm%h_sl_tr = lattice_interaction_SlipByTrans(N_sl,N_tr,&
+      prm%h_sl_tr = lattice_interaction_SlipByTrans(N_sl,N_tr,pl%get_as1dFloat('h_sl-tr'), &
-                                                    pl%get_as1dFloat('h_sl-tr'), &
+                                                    phase_lattice(ph))
-                                                    phase%get_asString('lattice'))
+      if (prm%fccTwinTransNucleation .and. size(N_tr) /= 1) extmsg = trim(extmsg)//' N_tr: nucleation'
-      if (prm%fccTwinTransNucleation .and. size(N_tr) /= 1) extmsg = trim(extmsg)//' interaction_sliptrans'
     endif slipAndTransActive
 
 !--------------------------------------------------------------------------------------------------
@@ -491,22 +481,21 @@ module function plastic_dislotwin_homogenizedC(ph,en) result(homogenizedC)
   real(pReal) :: f_unrotated
 
 
-  associate(prm => param(ph),&
+  associate(prm => param(ph), stt => state(ph))
-            stt => state(ph))
 
-  f_unrotated = 1.0_pReal &
+    f_unrotated = 1.0_pReal &
-              - sum(stt%f_tw(1:prm%sum_N_tw,en)) &
+                - sum(stt%f_tw(1:prm%sum_N_tw,en)) &
-              - sum(stt%f_tr(1:prm%sum_N_tr,en))
+                - sum(stt%f_tr(1:prm%sum_N_tr,en))
 
-  homogenizedC = f_unrotated * prm%C66
+    homogenizedC = f_unrotated * prm%C66
-  do i=1,prm%sum_N_tw
+    do i=1,prm%sum_N_tw
-    homogenizedC = homogenizedC &
+      homogenizedC = homogenizedC &
-                 + stt%f_tw(i,en)*prm%C66_tw(1:6,1:6,i)
+                   + stt%f_tw(i,en)*prm%C66_tw(1:6,1:6,i)
-  enddo
+    enddo
-  do i=1,prm%sum_N_tr
+    do i=1,prm%sum_N_tr
-    homogenizedC = homogenizedC &
+      homogenizedC = homogenizedC &
-                 + stt%f_tr(i,en)*prm%C66_tr(1:6,1:6,i)
+                   + stt%f_tr(i,en)*prm%C66_tr(1:6,1:6,i)
-  enddo
+    enddo
 
   end associate
 
@@ -531,7 +520,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
      ddot_gamma_dtau, &
      tau
   real(pReal), dimension(param(ph)%sum_N_sl) :: &
-    dot_gamma_sl,ddot_gamma_dtau_slip
+    dot_gamma_sl,ddot_gamma_dtau_sl
   real(pReal), dimension(param(ph)%sum_N_tw) :: &
     dot_gamma_tw,ddot_gamma_dtau_tw
   real(pReal), dimension(param(ph)%sum_N_tr) :: &
@@ -568,15 +557,15 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
   Lp = 0.0_pReal
   dLp_dMp = 0.0_pReal
 
-  call kinetics_slip(Mp,T,ph,en,dot_gamma_sl,ddot_gamma_dtau_slip)
+  call kinetics_sl(Mp,T,ph,en,dot_gamma_sl,ddot_gamma_dtau_sl)
   slipContribution: do i = 1, prm%sum_N_sl
     Lp = Lp + dot_gamma_sl(i)*prm%P_sl(1:3,1:3,i)
     forall (k=1:3,l=1:3,m=1:3,n=1:3) &
       dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
-                       + ddot_gamma_dtau_slip(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i)
+                       + ddot_gamma_dtau_sl(i) * prm%P_sl(k,l,i) * prm%P_sl(m,n,i)
   enddo slipContribution
 
-  call kinetics_twin(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw,ddot_gamma_dtau_tw)
+  call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw,ddot_gamma_dtau_tw)
   twinContibution: do i = 1, prm%sum_N_tw
     Lp = Lp + dot_gamma_tw(i)*prm%P_tw(1:3,1:3,i)
     forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@@ -584,7 +573,7 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
                        + ddot_gamma_dtau_tw(i)* prm%P_tw(k,l,i)*prm%P_tw(m,n,i)
   enddo twinContibution
 
-  call kinetics_trans(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr,ddot_gamma_dtau_tr)
+  call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr,ddot_gamma_dtau_tr)
   transContibution: do i = 1, prm%sum_N_tr
     Lp = Lp + dot_gamma_tr(i)*prm%P_tr(1:3,1:3,i)
     forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@@ -609,8 +598,8 @@ module subroutine dislotwin_LpAndItsTangent(Lp,dLp_dMp,Mp,T,ph,en)
         StressRatio_p = (abs(tau)/prm%xi_sb)**prm%p_sb
         dot_gamma_sb = sign(prm%v_sb*exp(-BoltzmannRatio*(1-StressRatio_p)**prm%q_sb), tau)
         ddot_gamma_dtau = abs(dot_gamma_sb)*BoltzmannRatio* prm%p_sb*prm%q_sb/ prm%xi_sb &
-                   * (abs(tau)/prm%xi_sb)**(prm%p_sb-1.0_pReal) &
+                        * (abs(tau)/prm%xi_sb)**(prm%p_sb-1.0_pReal) &
-                   * (1.0_pReal-StressRatio_p)**(prm%q_sb-1.0_pReal)
+                        * (1.0_pReal-StressRatio_p)**(prm%q_sb-1.0_pReal)
 
         Lp = Lp + dot_gamma_sb * P_sb
         forall (k=1:3,l=1:3,m=1:3,n=1:3) &
@@ -664,7 +653,7 @@ module subroutine dislotwin_dotState(Mp,T,ph,en)
               - sum(stt%f_tw(1:prm%sum_N_tw,en)) &
               - sum(stt%f_tr(1:prm%sum_N_tr,en))
 
-  call kinetics_slip(Mp,T,ph,en,dot_gamma_sl)
+  call kinetics_sl(Mp,T,ph,en,dot_gamma_sl)
   dot%gamma_sl(:,en) = abs(dot_gamma_sl)
 
   rho_dip_distance_min = prm%D_a*prm%b_sl
@@ -709,10 +698,10 @@ module subroutine dislotwin_dotState(Mp,T,ph,en)
                     - 2.0_pReal*rho_dip_distance_min/prm%b_sl * stt%rho_dip(:,en)*abs(dot_gamma_sl) &
                     - dot_rho_dip_climb
 
-  call kinetics_twin(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw)
+  call kinetics_tw(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tw)
   dot%f_tw(:,en) = f_unrotated*dot_gamma_tw/prm%gamma_char
 
-  call kinetics_trans(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr)
+  call kinetics_tr(Mp,T,dot_gamma_sl,ph,en,dot_gamma_tr)
   dot%f_tr(:,en) = f_unrotated*dot_gamma_tr
 
   end associate
@@ -745,54 +734,52 @@ module subroutine dislotwin_dependentState(T,ph,en)
     x0
 
 
-  associate(prm => param(ph),&
+  associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
-            stt => state(ph),&
-            dst => dependentState(ph))
 
-  sumf_tw  = sum(stt%f_tw(1:prm%sum_N_tw,en))
+    sumf_tw = sum(stt%f_tw(1:prm%sum_N_tw,en))
-  sumf_tr = sum(stt%f_tr(1:prm%sum_N_tr,en))
+    sumf_tr = sum(stt%f_tr(1:prm%sum_N_tr,en))
 
-  Gamma = prm%Gamma_sf(1) + prm%Gamma_sf(2) * T
+    Gamma = prm%Gamma_sf(1) + prm%Gamma_sf(2) * T
 
-  !* rescaled volume fraction for topology
+    !* rescaled volume fraction for topology
-  f_over_t_tw = stt%f_tw(1:prm%sum_N_tw,en)/prm%t_tw                                                ! this is per system ...
+    f_over_t_tw = stt%f_tw(1:prm%sum_N_tw,en)/prm%t_tw                                              ! this is per system ...
-  f_over_t_tr = sumf_tr/prm%t_tr                                                                    ! but this not
+    f_over_t_tr = sumf_tr/prm%t_tr                                                                  ! but this not
                                                                                                     ! ToDo ...Physically correct, but naming could be adjusted
 
-  inv_lambda_sl = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,en)+stt%rho_dip(:,en)))/prm%i_sl
+    inv_lambda_sl = sqrt(matmul(prm%forestProjection,stt%rho_mob(:,en)+stt%rho_dip(:,en)))/prm%i_sl
-  if (prm%sum_N_tw > 0 .and. prm%sum_N_sl > 0) &
+    if (prm%sum_N_tw > 0 .and. prm%sum_N_sl > 0) &
-    inv_lambda_sl = inv_lambda_sl + matmul(prm%h_sl_tw,f_over_t_tw)/(1.0_pReal-sumf_tw)
+      inv_lambda_sl = inv_lambda_sl + matmul(prm%h_sl_tw,f_over_t_tw)/(1.0_pReal-sumf_tw)
-  if (prm%sum_N_tr > 0 .and. prm%sum_N_sl > 0) &
+    if (prm%sum_N_tr > 0 .and. prm%sum_N_sl > 0) &
-    inv_lambda_sl = inv_lambda_sl + matmul(prm%h_sl_tr,f_over_t_tr)/(1.0_pReal-sumf_tr)
+      inv_lambda_sl = inv_lambda_sl + matmul(prm%h_sl_tr,f_over_t_tr)/(1.0_pReal-sumf_tr)
-  dst%Lambda_sl(:,en) = prm%D / (1.0_pReal+prm%D*inv_lambda_sl)
+    dst%Lambda_sl(:,en) = prm%D / (1.0_pReal+prm%D*inv_lambda_sl)
 
-  inv_lambda_tw_tw = matmul(prm%h_tw_tw,f_over_t_tw)/(1.0_pReal-sumf_tw)
+    inv_lambda_tw_tw = matmul(prm%h_tw_tw,f_over_t_tw)/(1.0_pReal-sumf_tw)
-  dst%Lambda_tw(:,en) = prm%i_tw*prm%D/(1.0_pReal+prm%D*inv_lambda_tw_tw)
+    dst%Lambda_tw(:,en) = prm%i_tw*prm%D/(1.0_pReal+prm%D*inv_lambda_tw_tw)
 
-  inv_lambda_tr_tr = matmul(prm%h_tr_tr,f_over_t_tr)/(1.0_pReal-sumf_tr)
+    inv_lambda_tr_tr = matmul(prm%h_tr_tr,f_over_t_tr)/(1.0_pReal-sumf_tr)
-  dst%Lambda_tr(:,en) = prm%i_tr*prm%D/(1.0_pReal+prm%D*inv_lambda_tr_tr)
+    dst%Lambda_tr(:,en) = prm%i_tr*prm%D/(1.0_pReal+prm%D*inv_lambda_tr_tr)
 
-  !* threshold stress for dislocation motion
+    !* threshold stress for dislocation motion
-  dst%tau_pass(:,en) = prm%mu*prm%b_sl* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en)))
+    dst%tau_pass(:,en) = prm%mu*prm%b_sl* sqrt(matmul(prm%h_sl_sl,stt%rho_mob(:,en)+stt%rho_dip(:,en)))
 
-  !* threshold stress for growing twin/martensite
+    !* threshold stress for growing twin/martensite
-  if(prm%sum_N_tw == prm%sum_N_sl) &
+    if(prm%sum_N_tw == prm%sum_N_sl) &
-    dst%tau_hat_tw(:,en) = Gamma/(3.0_pReal*prm%b_tw) &
+      dst%tau_hat_tw(:,en) = Gamma/(3.0_pReal*prm%b_tw) &
-                         + 3.0_pReal*prm%b_tw*prm%mu/(prm%L_tw*prm%b_sl) ! slip Burgers here correct?
+                           + 3.0_pReal*prm%b_tw*prm%mu/(prm%L_tw*prm%b_sl) ! slip Burgers here correct?
-  if(prm%sum_N_tr == prm%sum_N_sl) &
+    if(prm%sum_N_tr == prm%sum_N_sl) &
-    dst%tau_hat_tr(:,en) = Gamma/(3.0_pReal*prm%b_tr) &
+      dst%tau_hat_tr(:,en) = Gamma/(3.0_pReal*prm%b_tr) &
-                         + 3.0_pReal*prm%b_tr*prm%mu/(prm%L_tr*prm%b_sl) & ! slip Burgers here correct?
+                           + 3.0_pReal*prm%b_tr*prm%mu/(prm%L_tr*prm%b_sl) & ! slip Burgers here correct?
-                         + prm%h*prm%delta_G/ (3.0_pReal*prm%b_tr)
+                           + prm%h*prm%delta_G/ (3.0_pReal*prm%b_tr)
 
-  dst%V_tw(:,en) = (PI/4.0_pReal)*prm%t_tw*dst%Lambda_tw(:,en)**2.0_pReal
+    dst%V_tw(:,en) = (PI/4.0_pReal)*prm%t_tw*dst%Lambda_tw(:,en)**2.0_pReal
-  dst%V_tr(:,en) = (PI/4.0_pReal)*prm%t_tr*dst%Lambda_tr(:,en)**2.0_pReal
+    dst%V_tr(:,en) = (PI/4.0_pReal)*prm%t_tr*dst%Lambda_tr(:,en)**2.0_pReal
 
 
-  x0 = prm%mu*prm%b_tw**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu)        ! ToDo: In the paper, this is the Burgers vector for slip and is the same for twin and trans
+    x0 = prm%mu*prm%b_tw**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu)      ! ToDo: In the paper, this is the Burgers vector for slip and is the same for twin and trans
-  dst%tau_r_tw(:,en) = prm%mu*prm%b_tw/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(pi/3.0_pReal)/x0)
+    dst%tau_r_tw(:,en) = prm%mu*prm%b_tw/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tw)+cos(pi/3.0_pReal)/x0)
 
-  x0 = prm%mu*prm%b_tr**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu)        ! ToDo: In the paper, this is the Burgers vector for slip
+    x0 = prm%mu*prm%b_tr**2.0_pReal/(Gamma*8.0_pReal*PI)*(2.0_pReal+prm%nu)/(1.0_pReal-prm%nu)      ! ToDo: In the paper, this is the Burgers vector for slip
-  dst%tau_r_tr(:,en) = prm%mu*prm%b_tr/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(pi/3.0_pReal)/x0)
+    dst%tau_r_tr(:,en) = prm%mu*prm%b_tr/(2.0_pReal*PI)*(1.0_pReal/(x0+prm%x_c_tr)+cos(pi/3.0_pReal)/x0)
 
   end associate
 
@@ -857,8 +844,8 @@ end subroutine plastic_dislotwin_results
 ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
 ! have the optional arguments at the end
 !--------------------------------------------------------------------------------------------------
-pure subroutine kinetics_slip(Mp,T,ph,en, &
+pure subroutine kinetics_sl(Mp,T,ph,en, &
-                              dot_gamma_sl,ddot_gamma_dtau_slip,tau_slip)
+                            dot_gamma_sl,ddot_gamma_dtau_sl,tau_sl)
 
   real(pReal), dimension(3,3),  intent(in) :: &
     Mp                                                                                              !< Mandel stress
@@ -871,8 +858,8 @@ pure subroutine kinetics_slip(Mp,T,ph,en, &
   real(pReal), dimension(param(ph)%sum_N_sl), intent(out) :: &
     dot_gamma_sl
   real(pReal), dimension(param(ph)%sum_N_sl), optional, intent(out) :: &
-    ddot_gamma_dtau_slip, &
+    ddot_gamma_dtau_sl, &
-    tau_slip
+    tau_sl
   real(pReal), dimension(param(ph)%sum_N_sl) :: &
     ddot_gamma_dtau
 
@@ -891,9 +878,7 @@ pure subroutine kinetics_slip(Mp,T,ph,en, &
 
   associate(prm => param(ph), stt => state(ph), dst => dependentState(ph))
 
-  do i = 1, prm%sum_N_sl
+  tau = [(math_tensordot(Mp,prm%P_sl(1:3,1:3,i)),i = 1, prm%sum_N_sl)]
-    tau(i) = math_tensordot(Mp,prm%P_sl(1:3,1:3,i))
-  enddo
 
   tau_eff = abs(tau)-dst%tau_pass(:,en)
 
@@ -921,10 +906,10 @@ pure subroutine kinetics_slip(Mp,T,ph,en, &
 
   end associate
 
-  if(present(ddot_gamma_dtau_slip)) ddot_gamma_dtau_slip = ddot_gamma_dtau
+  if(present(ddot_gamma_dtau_sl)) ddot_gamma_dtau_sl = ddot_gamma_dtau
-  if(present(tau_slip))             tau_slip             = tau
+  if(present(tau_sl))             tau_sl             = tau
 
-end subroutine kinetics_slip
+end subroutine kinetics_sl
 
 
 !--------------------------------------------------------------------------------------------------
@@ -934,8 +919,8 @@ end subroutine kinetics_slip
 ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
 ! have the optional arguments at the end.
 !--------------------------------------------------------------------------------------------------
-pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,ph,en,&
+pure subroutine kinetics_tw(Mp,T,dot_gamma_sl,ph,en,&
-                              dot_gamma_tw,ddot_gamma_dtau_tw)
+                            dot_gamma_tw,ddot_gamma_dtau_tw)
 
   real(pReal), dimension(3,3),  intent(in) :: &
     Mp                                                                                              !< Mandel stress
@@ -993,7 +978,7 @@ pure subroutine kinetics_twin(Mp,T,dot_gamma_sl,ph,en,&
 
   if(present(ddot_gamma_dtau_tw)) ddot_gamma_dtau_tw = ddot_gamma_dtau
 
-end subroutine kinetics_twin
+end subroutine kinetics_tw
 
 
 !--------------------------------------------------------------------------------------------------
@@ -1003,8 +988,8 @@ end subroutine kinetics_twin
 ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
 ! have the optional arguments at the end.
 !--------------------------------------------------------------------------------------------------
-pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,ph,en,&
+pure subroutine kinetics_tr(Mp,T,dot_gamma_sl,ph,en,&
-                              dot_gamma_tr,ddot_gamma_dtau_tr)
+                            dot_gamma_tr,ddot_gamma_dtau_tr)
 
   real(pReal), dimension(3,3),  intent(in) :: &
     Mp                                                                                              !< Mandel stress
@@ -1061,6 +1046,6 @@ pure subroutine kinetics_trans(Mp,T,dot_gamma_sl,ph,en,&
 
   if(present(ddot_gamma_dtau_tr)) ddot_gamma_dtau_tr = ddot_gamma_dtau
 
-end subroutine kinetics_trans
+end subroutine kinetics_tr
 
 end submodule dislotwin

src/phase_mechanical_plastic_isotropic.f90

@@ -170,27 +170,28 @@ module subroutine isotropic_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
   integer :: &
     k, l, m, n
 
+
   associate(prm => param(ph), stt => state(ph))
 
-  Mp_dev = math_deviatoric33(Mp)
+    Mp_dev = math_deviatoric33(Mp)
-  squarenorm_Mp_dev = math_tensordot(Mp_dev,Mp_dev)
+    squarenorm_Mp_dev = math_tensordot(Mp_dev,Mp_dev)
-  norm_Mp_dev = sqrt(squarenorm_Mp_dev)
+    norm_Mp_dev = sqrt(squarenorm_Mp_dev)
 
-  if (norm_Mp_dev > 0.0_pReal) then
+    if (norm_Mp_dev > 0.0_pReal) then
-    dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp_dev/(prm%M*stt%xi(en))) **prm%n
+      dot_gamma = prm%dot_gamma_0 * (sqrt(1.5_pReal) * norm_Mp_dev/(prm%M*stt%xi(en))) **prm%n
 
-    Lp = dot_gamma * Mp_dev/norm_Mp_dev
+      Lp = dot_gamma * Mp_dev/norm_Mp_dev
-    forall (k=1:3,l=1:3,m=1:3,n=1:3) &
+      forall (k=1:3,l=1:3,m=1:3,n=1:3) &
-      dLp_dMp(k,l,m,n) = (prm%n-1.0_pReal) * Mp_dev(k,l)*Mp_dev(m,n) / squarenorm_Mp_dev
+        dLp_dMp(k,l,m,n) = (prm%n-1.0_pReal) * Mp_dev(k,l)*Mp_dev(m,n) / squarenorm_Mp_dev
-    forall (k=1:3,l=1:3) &
+      forall (k=1:3,l=1:3) &
-      dLp_dMp(k,l,k,l) = dLp_dMp(k,l,k,l) + 1.0_pReal
+        dLp_dMp(k,l,k,l) = dLp_dMp(k,l,k,l) + 1.0_pReal
-    forall (k=1:3,m=1:3) &
+      forall (k=1:3,m=1:3) &
-      dLp_dMp(k,k,m,m) = dLp_dMp(k,k,m,m) - 1.0_pReal/3.0_pReal
+        dLp_dMp(k,k,m,m) = dLp_dMp(k,k,m,m) - 1.0_pReal/3.0_pReal
-    dLp_dMp = dot_gamma * dLp_dMp / norm_Mp_dev
+      dLp_dMp = dot_gamma * dLp_dMp / norm_Mp_dev
-  else
+    else
-    Lp = 0.0_pReal
+      Lp = 0.0_pReal
-    dLp_dMp = 0.0_pReal
+      dLp_dMp = 0.0_pReal
-  end if
+    end if
 
   end associate
 
@@ -218,19 +219,20 @@ module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,ph,en)
   integer :: &
     k, l, m, n
 
+
   associate(prm => param(ph), stt => state(ph))
 
-  tr=math_trace33(math_spherical33(Mi))
+    tr=math_trace33(math_spherical33(Mi))
 
-  if (prm%dilatation .and. abs(tr) > 0.0_pReal) then                                                ! no stress or J2 plasticity --> Li and its derivative are zero
+    if (prm%dilatation .and. abs(tr) > 0.0_pReal) then                                              ! no stress or J2 plasticity --> Li and its derivative are zero
-    Li = math_I3 &
+      Li = math_I3 &
-       * prm%dot_gamma_0 * (3.0_pReal*prm%M*stt%xi(en))**(-prm%n) &
+         * prm%dot_gamma_0 * (3.0_pReal*prm%M*stt%xi(en))**(-prm%n) &
-       * tr * abs(tr)**(prm%n-1.0_pReal)
+         * tr * abs(tr)**(prm%n-1.0_pReal)
-    forall (k=1:3,l=1:3,m=1:3,n=1:3) dLi_dMi(k,l,m,n) = prm%n / tr * Li(k,l) * math_I3(m,n)
+      forall (k=1:3,l=1:3,m=1:3,n=1:3) dLi_dMi(k,l,m,n) = prm%n / tr * Li(k,l) * math_I3(m,n)
-  else
+    else
-    Li      = 0.0_pReal
+      Li      = 0.0_pReal
-    dLi_dMi = 0.0_pReal
+      dLi_dMi = 0.0_pReal
-  endif
+    endif
 
   end associate
 

src/phase_mechanical_plastic_kinehardening.f90

@@ -19,11 +19,11 @@ submodule(phase:plastic) kinehardening
       xi_inf_f, &
       xi_inf_b
     real(pReal),              allocatable, dimension(:,:) :: &
-      interaction_slipslip                                                                          !< slip resistance from slip activity
+      h_sl_sl                                                                                       !< slip resistance from slip activity
     real(pReal),              allocatable, dimension(:,:,:) :: &
       P, &
-      nonSchmid_pos, &
+      P_nS_pos, &
-      nonSchmid_neg
+      P_nS_neg
     integer :: &
       sum_N_sl
     logical :: &
@@ -33,13 +33,14 @@ submodule(phase:plastic) kinehardening
   end type tParameters
 
   type :: tKinehardeningState
-    real(pReal), pointer, dimension(:,:) :: &                                                       !< vectors along NipcMyInstance
+    real(pReal), pointer, dimension(:,:) :: &
-      crss, &                                                                                       !< critical resolved stress
+      xi, &                                                                                         !< resistance against plastic slip
-      crss_back, &                                                                                  !< critical resolved back stress
+      chi, &                                                                                        !< back stress
-      sense, &                                                                                      !< sense of acting shear stress (-1 or +1)
+      chi_0, &                                                                                      !< back stress at last switch of stress sense
-      chi0, &                                                                                       !< backstress at last switch of stress sense
+      gamma, &                                                                                      !< accumulated (absolute) shear
-      gamma0, &                                                                                     !< accumulated shear at last switch of stress sense
+      gamma_0, &                                                                                    !< accumulated shear at last switch of stress sense
-      accshear                                                                                      !< accumulated (absolute) shear
+      sgn_gamma                                                                                     !< sense of acting shear stress (-1 or +1)
+
   end type tKinehardeningState
 
 !--------------------------------------------------------------------------------------------------
@@ -112,21 +113,19 @@ module function plastic_kinehardening_init() result(myPlasticity)
     N_sl         = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray)
     prm%sum_N_sl = sum(abs(N_sl))
     slipActive: if (prm%sum_N_sl > 0) then
-      prm%P = lattice_SchmidMatrix_slip(N_sl,phase%get_asString('lattice'),&
+      prm%P = lattice_SchmidMatrix_slip(N_sl,phase_lattice(ph),phase_cOverA(ph))
-                                        phase%get_asFloat('c/a',defaultVal=0.0_pReal))
 
-      if(trim(phase%get_asString('lattice')) == 'cI') then
+      if (phase_lattice(ph) == 'cI') then
         a = pl%get_as1dFloat('a_nonSchmid',defaultVal = emptyRealArray)
         if(size(a) > 0) prm%nonSchmidActive = .true.
-        prm%nonSchmid_pos  = lattice_nonSchmidMatrix(N_sl,a,+1)
+        prm%P_nS_pos = lattice_nonSchmidMatrix(N_sl,a,+1)
-        prm%nonSchmid_neg  = lattice_nonSchmidMatrix(N_sl,a,-1)
+        prm%P_nS_neg = lattice_nonSchmidMatrix(N_sl,a,-1)
       else
-        prm%nonSchmid_pos  = prm%P
+        prm%P_nS_pos = prm%P
-        prm%nonSchmid_neg  = prm%P
+        prm%P_nS_neg = prm%P
       endif
-      prm%interaction_SlipSlip = lattice_interaction_SlipBySlip(N_sl, &
+      prm%h_sl_sl = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dFloat('h_sl-sl'), &
-                                                                pl%get_as1dFloat('h_sl-sl'), &
+                                                   phase_lattice(ph))
-                                                                phase%get_asString('lattice'))
 
       xi_0          = pl%get_as1dFloat('xi_0',       requiredSize=size(N_sl))
       prm%xi_inf_f  = pl%get_as1dFloat('xi_inf_f',   requiredSize=size(N_sl))
@@ -156,18 +155,17 @@ module function plastic_kinehardening_init() result(myPlasticity)
       if (any(prm%xi_inf_f     <= 0.0_pReal))  extmsg = trim(extmsg)//' xi_inf_f'
       if (any(prm%xi_inf_b     <= 0.0_pReal))  extmsg = trim(extmsg)//' xi_inf_b'
 
-      !ToDo: Any sensible checks for theta?
     else slipActive
       xi_0 = emptyRealArray
       allocate(prm%xi_inf_f,prm%xi_inf_b,prm%h_0_f,prm%h_inf_f,prm%h_0_b,prm%h_inf_b,source=emptyRealArray)
-      allocate(prm%interaction_SlipSlip(0,0))
+      allocate(prm%h_sl_sl(0,0))
     endif slipActive
 
 !--------------------------------------------------------------------------------------------------
 ! allocate state arrays
     Nmembers = count(material_phaseID == ph)
-    sizeDotState   = size(['crss     ','crss_back', 'accshear ']) * prm%sum_N_sl !ToDo: adjust names like in material.yaml
+    sizeDotState   = size(['xi   ','chi  ', 'gamma']) * prm%sum_N_sl
-    sizeDeltaState = size(['sense ',   'chi0  ',    'gamma0'   ]) * prm%sum_N_sl !ToDo: adjust names like in material.yaml
+    sizeDeltaState = size(['sgn_gamma',   'chi_0    ',    'gamma_0  ']) * prm%sum_N_sl
     sizeState = sizeDotState + sizeDeltaState
 
     call phase_allocateState(plasticState(ph),Nmembers,sizeState,sizeDotState,sizeDeltaState)
@@ -176,40 +174,40 @@ module function plastic_kinehardening_init() result(myPlasticity)
 ! state aliases and initialization
     startIndex = 1
     endIndex   = prm%sum_N_sl
-    stt%crss => plasticState(ph)%state   (startIndex:endIndex,:)
+    stt%xi => plasticState(ph)%state   (startIndex:endIndex,:)
-    stt%crss = spread(xi_0, 2, Nmembers)
+    stt%xi = spread(xi_0, 2, Nmembers)
-    dot%crss => plasticState(ph)%dotState(startIndex:endIndex,:)
+    dot%xi => plasticState(ph)%dotState(startIndex:endIndex,:)
     plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
     if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
 
     startIndex = endIndex + 1
     endIndex   = endIndex + prm%sum_N_sl
-    stt%crss_back => plasticState(ph)%state   (startIndex:endIndex,:)
+    stt%chi => plasticState(ph)%state   (startIndex:endIndex,:)
-    dot%crss_back => plasticState(ph)%dotState(startIndex:endIndex,:)
+    dot%chi => plasticState(ph)%dotState(startIndex:endIndex,:)
     plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
 
     startIndex = endIndex + 1
     endIndex   = endIndex + prm%sum_N_sl
-    stt%accshear => plasticState(ph)%state   (startIndex:endIndex,:)
+    stt%gamma => plasticState(ph)%state   (startIndex:endIndex,:)
-    dot%accshear => plasticState(ph)%dotState(startIndex:endIndex,:)
+    dot%gamma => plasticState(ph)%dotState(startIndex:endIndex,:)
     plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
     if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
 
     o = plasticState(ph)%offsetDeltaState
     startIndex = endIndex + 1
     endIndex   = endIndex + prm%sum_N_sl
-    stt%sense => plasticState(ph)%state     (startIndex  :endIndex  ,:)
+    stt%sgn_gamma => plasticState(ph)%state     (startIndex  :endIndex  ,:)
-    dlt%sense => plasticState(ph)%deltaState(startIndex-o:endIndex-o,:)
+    dlt%sgn_gamma => plasticState(ph)%deltaState(startIndex-o:endIndex-o,:)
 
     startIndex = endIndex + 1
     endIndex   = endIndex +  prm%sum_N_sl
-    stt%chi0 => plasticState(ph)%state     (startIndex  :endIndex  ,:)
+    stt%chi_0 => plasticState(ph)%state     (startIndex  :endIndex  ,:)
-    dlt%chi0 => plasticState(ph)%deltaState(startIndex-o:endIndex-o,:)
+    dlt%chi_0 => plasticState(ph)%deltaState(startIndex-o:endIndex-o,:)
 
     startIndex = endIndex + 1
     endIndex   = endIndex +  prm%sum_N_sl
-    stt%gamma0 => plasticState(ph)%state     (startIndex  :endIndex  ,:)
+    stt%gamma_0 => plasticState(ph)%state     (startIndex  :endIndex  ,:)
-    dlt%gamma0 => plasticState(ph)%deltaState(startIndex-o:endIndex-o,:)
+    dlt%gamma_0 => plasticState(ph)%deltaState(startIndex-o:endIndex-o,:)
 
     end associate
 
@@ -242,22 +240,22 @@ pure module subroutine kinehardening_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
   integer :: &
     i,k,l,m,n
   real(pReal), dimension(param(ph)%sum_N_sl) :: &
-    gdot_pos,gdot_neg, &
+    dot_gamma_pos,dot_gamma_neg, &
-    dgdot_dtau_pos,dgdot_dtau_neg
+    ddot_gamma_dtau_pos,ddot_gamma_dtau_neg
 
   Lp = 0.0_pReal
   dLp_dMp = 0.0_pReal
 
   associate(prm => param(ph))
 
-  call kinetics(Mp,ph,en,gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg)
+  call kinetics(Mp,ph,en,dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg)
 
   do i = 1, prm%sum_N_sl
-    Lp = Lp + (gdot_pos(i)+gdot_neg(i))*prm%P(1:3,1:3,i)
+    Lp = Lp + (dot_gamma_pos(i)+dot_gamma_neg(i))*prm%P(1:3,1:3,i)
     forall (k=1:3,l=1:3,m=1:3,n=1:3) &
       dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
-                       + dgdot_dtau_pos(i) * prm%P(k,l,i) * prm%nonSchmid_pos(m,n,i) &
+                       + ddot_gamma_dtau_pos(i) * prm%P(k,l,i) * prm%P_nS_pos(m,n,i) &
-                       + dgdot_dtau_neg(i) * prm%P(k,l,i) * prm%nonSchmid_neg(m,n,i)
+                       + ddot_gamma_dtau_neg(i) * prm%P(k,l,i) * prm%P_nS_neg(m,n,i)
   enddo
 
   end associate
@@ -279,29 +277,28 @@ module subroutine plastic_kinehardening_dotState(Mp,ph,en)
   real(pReal) :: &
     sumGamma
   real(pReal), dimension(param(ph)%sum_N_sl) :: &
-    gdot_pos,gdot_neg
+    dot_gamma_pos,dot_gamma_neg
 
 
-  associate(prm => param(ph), stt => state(ph),&
+  associate(prm => param(ph), stt => state(ph),dot => dotState(ph))
-            dot => dotState(ph))
 
-  call kinetics(Mp,ph,en,gdot_pos,gdot_neg)
+    call kinetics(Mp,ph,en,dot_gamma_pos,dot_gamma_neg)
-  dot%accshear(:,en) = abs(gdot_pos+gdot_neg)
+    dot%gamma(:,en) = abs(dot_gamma_pos+dot_gamma_neg)
-  sumGamma = sum(stt%accshear(:,en))
+    sumGamma = sum(stt%gamma(:,en))
 
 
-  dot%crss(:,en) = matmul(prm%interaction_SlipSlip,dot%accshear(:,en)) &
+    dot%xi(:,en) = matmul(prm%h_sl_sl,dot%gamma(:,en)) &
-                 * (  prm%h_inf_f &
+                   * (  prm%h_inf_f &
-                     + (prm%h_0_f - prm%h_inf_f + prm%h_0_f*prm%h_inf_f*sumGamma/prm%xi_inf_f) &
+                       + (prm%h_0_f - prm%h_inf_f + prm%h_0_f*prm%h_inf_f*sumGamma/prm%xi_inf_f) &
-                     * exp(-sumGamma*prm%h_0_f/prm%xi_inf_f) &
+                       * exp(-sumGamma*prm%h_0_f/prm%xi_inf_f) &
-                   )
+                     )
 
-  dot%crss_back(:,en) = stt%sense(:,en)*dot%accshear(:,en) * &
+    dot%chi(:,en) = stt%sgn_gamma(:,en)*dot%gamma(:,en) * &
-           ( prm%h_inf_b + &
+             ( prm%h_inf_b + &
-             (prm%h_0_b - prm%h_inf_b &
+               (prm%h_0_b - prm%h_inf_b &
-               + prm%h_0_b*prm%h_inf_b/(prm%xi_inf_b+stt%chi0(:,en))*(stt%accshear(:,en)-stt%gamma0(:,en))&
+                 + prm%h_0_b*prm%h_inf_b/(prm%xi_inf_b+stt%chi_0(:,en))*(stt%gamma(:,en)-stt%gamma_0(:,en))&
-             ) *exp(-(stt%accshear(:,en)-stt%gamma0(:,en)) *prm%h_0_b/(prm%xi_inf_b+stt%chi0(:,en))) &
+               ) *exp(-(stt%gamma(:,en)-stt%gamma_0(:,en)) *prm%h_0_b/(prm%xi_inf_b+stt%chi_0(:,en))) &
-           )
+             )
 
   end associate
 
@@ -320,28 +317,26 @@ module subroutine plastic_kinehardening_deltaState(Mp,ph,en)
     en
 
   real(pReal), dimension(param(ph)%sum_N_sl) :: &
-    gdot_pos,gdot_neg, &
+    dot_gamma_pos,dot_gamma_neg, &
-    sense
+    sgn_gamma
+
 
   associate(prm => param(ph), stt => state(ph), dlt => deltaState(ph))
 
-  call kinetics(Mp,ph,en,gdot_pos,gdot_neg)
+    call kinetics(Mp,ph,en,dot_gamma_pos,dot_gamma_neg)
-  sense = merge(state(ph)%sense(:,en), &                                                            ! keep existing...
+    sgn_gamma = merge(state(ph)%sgn_gamma(:,en), &
-                sign(1.0_pReal,gdot_pos+gdot_neg), &                                                ! ...or have a defined
+                      sign(1.0_pReal,dot_gamma_pos+dot_gamma_neg), &
-                dEq0(gdot_pos+gdot_neg,1e-10_pReal))                                                ! current sense of shear direction
+                      dEq0(dot_gamma_pos+dot_gamma_neg,1e-10_pReal))
 
-
+    where(dNeq(sgn_gamma,stt%sgn_gamma(:,en),0.1_pReal)) ! ToDo sgn_gamma*stt%sgn_gamma(:,en)<0
-!--------------------------------------------------------------------------------------------------
+      dlt%sgn_gamma (:,en) = sgn_gamma - stt%sgn_gamma(:,en)
-! switch in sense of shear?
+      dlt%chi_0  (:,en) = abs(stt%chi(:,en)) - stt%chi_0(:,en)
-  where(dNeq(sense,stt%sense(:,en),0.1_pReal))
+      dlt%gamma_0(:,en) = stt%gamma(:,en) - stt%gamma_0(:,en)
-    dlt%sense (:,en) = sense - stt%sense(:,en)                                                      ! switch sense
+    else where
-    dlt%chi0  (:,en) = abs(stt%crss_back(:,en)) - stt%chi0(:,en)                                    ! remember current backstress magnitude
+      dlt%sgn_gamma (:,en) = 0.0_pReal
-    dlt%gamma0(:,en) = stt%accshear(:,en) - stt%gamma0(:,en)                                        ! remember current accumulated shear
+      dlt%chi_0  (:,en) = 0.0_pReal
-  else where
+      dlt%gamma_0(:,en) = 0.0_pReal
-    dlt%sense (:,en) = 0.0_pReal
+    end where
-    dlt%chi0  (:,en) = 0.0_pReal
-    dlt%gamma0(:,en) = 0.0_pReal
-  end where
 
   end associate
 
@@ -362,22 +357,22 @@ module subroutine plastic_kinehardening_results(ph,group)
   outputsLoop: do o = 1,size(prm%output)
     select case(trim(prm%output(o)))
      case ('xi')
-       if(prm%sum_N_sl>0) call results_writeDataset(stt%crss,group,trim(prm%output(o)), &
+       if(prm%sum_N_sl>0) call results_writeDataset(stt%xi,group,trim(prm%output(o)), &
                                                     'resistance against plastic slip','Pa')
-     case ('tau_b')
+     case ('tau_b') !ToDo: chi
-       if(prm%sum_N_sl>0) call results_writeDataset(stt%crss_back,group,trim(prm%output(o)), &
+       if(prm%sum_N_sl>0) call results_writeDataset(stt%chi,group,trim(prm%output(o)), &
-                                                    'back stress against plastic slip','Pa')
+                                                    'back stress','Pa')
      case ('sgn(gamma)')
-       if(prm%sum_N_sl>0) call results_writeDataset(stt%sense,group,trim(prm%output(o)), & ! ToDo: could be int
+       if(prm%sum_N_sl>0) call results_writeDataset(stt%sgn_gamma,group,trim(prm%output(o)), & ! ToDo: could be int
                                                     'sense of shear','1')
      case ('chi_0')
-       if(prm%sum_N_sl>0) call results_writeDataset(stt%chi0,group,trim(prm%output(o)), &
+       if(prm%sum_N_sl>0) call results_writeDataset(stt%chi_0,group,trim(prm%output(o)), &
-                                                    'tbd','Pa')
+                                                    'back stress at last switch of stress sense','Pa')
      case ('gamma_0')
-       if(prm%sum_N_sl>0) call results_writeDataset(stt%gamma0,group,trim(prm%output(o)), &
+       if(prm%sum_N_sl>0) call results_writeDataset(stt%gamma_0,group,trim(prm%output(o)), &
-                                                    'tbd','1')
+                                                    'plastic shear at last switch of stress sense','1')
      case ('gamma')
-       if(prm%sum_N_sl>0) call results_writeDataset(stt%accshear,group,trim(prm%output(o)), &
+       if(prm%sum_N_sl>0) call results_writeDataset(stt%gamma,group,trim(prm%output(o)), &
                                                     'plastic shear','1')
     end select
   enddo outputsLoop
@@ -394,62 +389,64 @@ end subroutine plastic_kinehardening_results
 ! have the optional arguments at the end.
 !--------------------------------------------------------------------------------------------------
 pure subroutine kinetics(Mp,ph,en, &
-                         gdot_pos,gdot_neg,dgdot_dtau_pos,dgdot_dtau_neg)
+                         dot_gamma_pos,dot_gamma_neg,ddot_gamma_dtau_pos,ddot_gamma_dtau_neg)
 
   real(pReal), dimension(3,3),  intent(in) :: &
     Mp                                                                                              !< Mandel stress
-  integer,                intent(in) :: &
+  integer,                      intent(in) :: &
     ph, &
     en
 
   real(pReal),                  intent(out), dimension(param(ph)%sum_N_sl) :: &
-    gdot_pos, &
+    dot_gamma_pos, &
-    gdot_neg
+    dot_gamma_neg
-  real(pReal),                  intent(out), optional, dimension(param(ph)%sum_N_sl) :: &
+  real(pReal),                  intent(out), dimension(param(ph)%sum_N_sl), optional :: &
-    dgdot_dtau_pos, &
+    ddot_gamma_dtau_pos, &
-    dgdot_dtau_neg
+    ddot_gamma_dtau_neg
 
   real(pReal), dimension(param(ph)%sum_N_sl) :: &
     tau_pos, &
     tau_neg
   integer :: i
 
+
   associate(prm => param(ph), stt => state(ph))
 
-  do i = 1, prm%sum_N_sl
+    do i = 1, prm%sum_N_sl
-    tau_pos(i) =       math_tensordot(Mp,prm%nonSchmid_pos(1:3,1:3,i)) - stt%crss_back(i,en)
+      tau_pos(i) =       math_tensordot(Mp,prm%P_nS_pos(1:3,1:3,i)) - stt%chi(i,en)
-    tau_neg(i) = merge(math_tensordot(Mp,prm%nonSchmid_neg(1:3,1:3,i)) - stt%crss_back(i,en), &
+      tau_neg(i) = merge(math_tensordot(Mp,prm%P_nS_neg(1:3,1:3,i)) - stt%chi(i,en), &
-                       0.0_pReal, prm%nonSchmidActive)
+                         0.0_pReal, prm%nonSchmidActive)
-  enddo
+    enddo
 
-  where(dNeq0(tau_pos))
+    where(dNeq0(tau_pos))
-    gdot_pos = prm%dot_gamma_0 * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) &                  ! 1/2 if non-Schmid active
+      dot_gamma_pos = prm%dot_gamma_0 * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) &           ! 1/2 if non-Schmid active
-             * sign(abs(tau_pos/stt%crss(:,en))**prm%n,  tau_pos)
+                    * sign(abs(tau_pos/stt%xi(:,en))**prm%n,  tau_pos)
-  else where
-    gdot_pos = 0.0_pReal
-  end where
-
-  where(dNeq0(tau_neg))
-    gdot_neg = prm%dot_gamma_0 * 0.5_pReal &                                                        ! only used if non-Schmid active, always 1/2
-             * sign(abs(tau_neg/stt%crss(:,en))**prm%n,  tau_neg)
-  else where
-    gdot_neg = 0.0_pReal
-  end where
-
-  if (present(dgdot_dtau_pos)) then
-    where(dNeq0(gdot_pos))
-      dgdot_dtau_pos = gdot_pos*prm%n/tau_pos
     else where
-      dgdot_dtau_pos = 0.0_pReal
+      dot_gamma_pos = 0.0_pReal
     end where
-  endif
+
-  if (present(dgdot_dtau_neg)) then
+    where(dNeq0(tau_neg))
-    where(dNeq0(gdot_neg))
+      dot_gamma_neg = prm%dot_gamma_0 * 0.5_pReal &                                                 ! only used if non-Schmid active, always 1/2
-      dgdot_dtau_neg = gdot_neg*prm%n/tau_neg
+                    * sign(abs(tau_neg/stt%xi(:,en))**prm%n,  tau_neg)
     else where
-      dgdot_dtau_neg = 0.0_pReal
+      dot_gamma_neg = 0.0_pReal
     end where
-  endif
+
+    if (present(ddot_gamma_dtau_pos)) then
+      where(dNeq0(dot_gamma_pos))
+        ddot_gamma_dtau_pos = dot_gamma_pos*prm%n/tau_pos
+      else where
+        ddot_gamma_dtau_pos = 0.0_pReal
+      end where
+    endif
+    if (present(ddot_gamma_dtau_neg)) then
+      where(dNeq0(dot_gamma_neg))
+        ddot_gamma_dtau_neg = dot_gamma_neg*prm%n/tau_neg
+      else where
+        ddot_gamma_dtau_neg = 0.0_pReal
+      end where
+    endif
+
   end associate
 
 end subroutine kinetics

src/phase_mechanical_plastic_nonlocal.f90

@@ -108,9 +108,9 @@ submodule(phase:plastic) nonlocal
       forestProjection_Edge, &                                                                      !< matrix of forest projections of edge dislocations
       forestProjection_Screw                                                                        !< matrix of forest projections of screw dislocations
     real(pReal), dimension(:,:,:), allocatable :: &
-      Schmid, &                                                                                     !< Schmid contribution
+      P_sl, &                                                                                       !< Schmid contribution
-      nonSchmid_pos, &
+      P_nS_pos, &
-      nonSchmid_neg                                                                                 !< combined projection of Schmid and non-Schmid contributions to the resolved shear stress (only for screws)
+      P_nS_neg                                                                                      !< combined projection of Schmid and non-Schmid contributions to the resolved shear stress (only for screws)
     integer :: &
       sum_N_sl = 0
     integer,     dimension(:),     allocatable :: &
@@ -240,41 +240,35 @@ 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 = elastic_mu(ph)
-    prm%mu = lattice_mu(ph)
+    prm%nu = elastic_nu(ph)
-    prm%nu = lattice_nu(ph)
 
     ini%N_sl     = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray)
     prm%sum_N_sl = sum(abs(ini%N_sl))
     slipActive: if (prm%sum_N_sl > 0) then
-      prm%Schmid = lattice_SchmidMatrix_slip(ini%N_sl,phase%get_asString('lattice'),&
+      prm%P_sl = lattice_SchmidMatrix_slip(ini%N_sl,phase_lattice(ph), phase_cOverA(ph))
-                                             phase%get_asFloat('c/a',defaultVal=0.0_pReal))
 
-      if(trim(phase%get_asString('lattice')) == 'cI') then
+      if (phase_lattice(ph) == 'cI') then
         a = pl%get_as1dFloat('a_nonSchmid',defaultVal = emptyRealArray)
         if(size(a) > 0) prm%nonSchmidActive = .true.
-        prm%nonSchmid_pos = lattice_nonSchmidMatrix(ini%N_sl,a,+1)
+        prm%P_nS_pos = lattice_nonSchmidMatrix(ini%N_sl,a,+1)
-        prm%nonSchmid_neg = lattice_nonSchmidMatrix(ini%N_sl,a,-1)
+        prm%P_nS_neg = lattice_nonSchmidMatrix(ini%N_sl,a,-1)
       else
-        prm%nonSchmid_pos = prm%Schmid
+        prm%P_nS_pos = prm%P_sl
-        prm%nonSchmid_neg = prm%Schmid
+        prm%P_nS_neg = prm%P_sl
       endif
 
-      prm%h_sl_sl = lattice_interaction_SlipBySlip(ini%N_sl, &
+      prm%h_sl_sl = lattice_interaction_SlipBySlip(ini%N_sl,pl%get_as1dFloat('h_sl-sl'), &
-                                                   pl%get_as1dFloat('h_sl-sl'), &
+                                                   phase_lattice(ph))
-                                                   phase%get_asString('lattice'))
 
-      prm%forestProjection_edge  = lattice_forestProjection_edge (ini%N_sl,phase%get_asString('lattice'),&
+      prm%forestProjection_edge  = lattice_forestProjection_edge (ini%N_sl,phase_lattice(ph),&
-                                                                  phase%get_asFloat('c/a',defaultVal=0.0_pReal))
+                                                                  phase_cOverA(ph))
-      prm%forestProjection_screw = lattice_forestProjection_screw(ini%N_sl,phase%get_asString('lattice'),&
+      prm%forestProjection_screw = lattice_forestProjection_screw(ini%N_sl,phase_lattice(ph),&
-                                                                  phase%get_asFloat('c/a',defaultVal=0.0_pReal))
+                                                                  phase_cOverA(ph))
 
-      prm%slip_direction  = lattice_slip_direction (ini%N_sl,phase%get_asString('lattice'),&
+      prm%slip_direction  = lattice_slip_direction (ini%N_sl,phase_lattice(ph),phase_cOverA(ph))
-                                                    phase%get_asFloat('c/a',defaultVal=0.0_pReal))
+      prm%slip_transverse = lattice_slip_transverse(ini%N_sl,phase_lattice(ph),phase_cOverA(ph))
-      prm%slip_transverse = lattice_slip_transverse(ini%N_sl,phase%get_asString('lattice'),&
+      prm%slip_normal     = lattice_slip_normal    (ini%N_sl,phase_lattice(ph),phase_cOverA(ph))
-                                                    phase%get_asFloat('c/a',defaultVal=0.0_pReal))
-      prm%slip_normal     = lattice_slip_normal    (ini%N_sl,phase%get_asString('lattice'),&
-                                                    phase%get_asFloat('c/a',defaultVal=0.0_pReal))
 
       ! collinear systems (only for octahedral slip systems in fcc)
       allocate(prm%colinearSystem(prm%sum_N_sl), source = -1)
@@ -761,11 +755,7 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
     Mp
                                                                                                     !< derivative of Lp with respect to Mp
   integer :: &
-    ns, &                                                                                           !< short notation for the total number of active slip systems
+    i, j, k, l, &
-    i, &
-    j, &
-    k, &
-    l, &
     t, &                                                                                            !< dislocation type
     s                                                                                               !< index of my current slip system
   real(pReal), dimension(param(ph)%sum_N_sl,8) :: &
@@ -779,26 +769,24 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
     dv_dtauNS                                                                                       !< velocity derivative with respect to the shear stress
   real(pReal), dimension(param(ph)%sum_N_sl) :: &
     tau, &                                                                                          !< resolved shear stress including backstress terms
-    gdotTotal                                                                                       !< shear rate
+    dot_gamma                                                                                       !< shear rate
 
-  associate(prm => param(ph),dst=>microstructure(ph),&
+  associate(prm => param(ph),dst=>microstructure(ph),stt=>state(ph))
-  stt=>state(ph))
-  ns = prm%sum_N_sl
 
   !*** shortcut to state variables
   rho = getRho(ph,en)
   rhoSgl = rho(:,sgl)
 
-  do s = 1,ns
+  do s = 1,prm%sum_N_sl
-    tau(s) = math_tensordot(Mp, prm%Schmid(1:3,1:3,s))
+    tau(s) = math_tensordot(Mp, prm%P_sl(1:3,1:3,s))
     tauNS(s,1) = tau(s)
     tauNS(s,2) = tau(s)
     if (tau(s) > 0.0_pReal) then
-      tauNS(s,3) = math_tensordot(Mp, +prm%nonSchmid_pos(1:3,1:3,s))
+      tauNS(s,3) = math_tensordot(Mp, +prm%P_nS_pos(1:3,1:3,s))
-      tauNS(s,4) = math_tensordot(Mp, -prm%nonSchmid_neg(1:3,1:3,s))
+      tauNS(s,4) = math_tensordot(Mp, -prm%P_nS_neg(1:3,1:3,s))
     else
-      tauNS(s,3) = math_tensordot(Mp, +prm%nonSchmid_neg(1:3,1:3,s))
+      tauNS(s,3) = math_tensordot(Mp, +prm%P_nS_neg(1:3,1:3,s))
-      tauNS(s,4) = math_tensordot(Mp, -prm%nonSchmid_pos(1:3,1:3,s))
+      tauNS(s,4) = math_tensordot(Mp, -prm%P_nS_pos(1:3,1:3,s))
     endif
   enddo
   tauNS = tauNS + spread(dst%tau_back(:,en),2,4)
@@ -826,22 +814,22 @@ module subroutine nonlocal_LpAndItsTangent(Lp,dLp_dMp, &
   stt%v(:,en) = pack(v,.true.)
 
   !*** Bauschinger effect
-  forall (s = 1:ns, t = 5:8, rhoSgl(s,t) * v(s,t-4) < 0.0_pReal) &
+  forall (s = 1:prm%sum_N_sl, t = 5:8, rhoSgl(s,t) * v(s,t-4) < 0.0_pReal) &
     rhoSgl(s,t-4) = rhoSgl(s,t-4) + abs(rhoSgl(s,t))
 
-  gdotTotal = sum(rhoSgl(:,1:4) * v, 2) * prm%b_sl
+  dot_gamma = sum(rhoSgl(:,1:4) * v, 2) * prm%b_sl
 
   Lp = 0.0_pReal
   dLp_dMp = 0.0_pReal
-  do s = 1,ns
+  do s = 1,prm%sum_N_sl
-    Lp = Lp + gdotTotal(s) * prm%Schmid(1:3,1:3,s)
+    Lp = Lp + dot_gamma(s) * prm%P_sl(1:3,1:3,s)
     forall (i=1:3,j=1:3,k=1:3,l=1:3) &
       dLp_dMp(i,j,k,l) = dLp_dMp(i,j,k,l) &
-          + prm%Schmid(i,j,s) * prm%Schmid(k,l,s) &
+          + prm%P_sl(i,j,s) * prm%P_sl(k,l,s) &
           * sum(rhoSgl(s,1:4) * dv_dtau(s,1:4)) * prm%b_sl(s) &
-          + prm%Schmid(i,j,s) &
+          + prm%P_sl(i,j,s) &
-          * (+ prm%nonSchmid_pos(k,l,s) * rhoSgl(s,3) * dv_dtauNS(s,3) &
+          * (+ prm%P_nS_pos(k,l,s) * rhoSgl(s,3) * dv_dtauNS(s,3) &
-             - prm%nonSchmid_neg(k,l,s) * rhoSgl(s,4) * dv_dtauNS(s,4))  * prm%b_sl(s)
+             - prm%P_nS_neg(k,l,s) * rhoSgl(s,4) * dv_dtauNS(s,4))  * prm%b_sl(s)
   enddo
 
   end associate
@@ -861,7 +849,6 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en)
     en
 
   integer :: &
-    ns, &                                                                                           ! short notation for the total number of active slip systems
     c, &                                                                                            ! character of dislocation
     t, &                                                                                            ! type of dislocation
     s                                                                                               ! index of my current slip system
@@ -881,11 +868,10 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en)
     deltaDUpper                                                                                     ! change in maximum stable dipole distance for edges and screws
 
   associate(prm => param(ph),dst => microstructure(ph),del => deltaState(ph))
-  ns = prm%sum_N_sl
 
   !*** shortcut to state variables
-  forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en)
+  forall (s = 1:prm%sum_N_sl, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en)
-  forall (s = 1:ns, c = 1:2) dUpperOld(s,c) = plasticState(ph)%state(iD(s,c,ph),en)
+  forall (s = 1:prm%sum_N_sl, c = 1:2) dUpperOld(s,c) = plasticState(ph)%state(iD(s,c,ph),en)
 
   rho =  getRho(ph,en)
   rhoDip = rho(:,dip)
@@ -908,7 +894,7 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en)
 
   !*** calculate limits for stable dipole height
   do s = 1,prm%sum_N_sl
-    tau(s) = math_tensordot(Mp, prm%Schmid(1:3,1:3,s)) +dst%tau_back(s,en)
+    tau(s) = math_tensordot(Mp, prm%P_sl(1:3,1:3,s)) +dst%tau_back(s,en)
     if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal
   enddo
 
@@ -926,16 +912,16 @@ module subroutine plastic_nonlocal_deltaState(Mp,ph,en)
 
   !*** dissociation by stress increase
   deltaRhoDipole2SingleStress = 0.0_pReal
-  forall (c=1:2, s=1:ns, deltaDUpper(s,c) < 0.0_pReal .and. &
+  forall (c=1:2, s=1:prm%sum_N_sl, deltaDUpper(s,c) < 0.0_pReal .and. &
                                           dNeq0(dUpperOld(s,c) - prm%minDipoleHeight(s,c))) &
     deltaRhoDipole2SingleStress(s,8+c) = rhoDip(s,c) * deltaDUpper(s,c) &
                                        / (dUpperOld(s,c) - prm%minDipoleHeight(s,c))
 
   forall (t=1:4) deltaRhoDipole2SingleStress(:,t) = -0.5_pReal * deltaRhoDipole2SingleStress(:,(t-1)/2+9)
-  forall (s = 1:ns, c = 1:2) plasticState(ph)%state(iD(s,c,ph),en) = dUpper(s,c)
+  forall (s = 1:prm%sum_N_sl, c = 1:2) plasticState(ph)%state(iD(s,c,ph),en) = dUpper(s,c)
 
   plasticState(ph)%deltaState(:,en) = 0.0_pReal
-  del%rho(:,en) = reshape(deltaRhoRemobilization + deltaRhoDipole2SingleStress, [10*ns])
+  del%rho(:,en) = reshape(deltaRhoRemobilization + deltaRhoDipole2SingleStress, [10*prm%sum_N_sl])
 
   end associate
 
@@ -946,7 +932,7 @@ end subroutine plastic_nonlocal_deltaState
 !> @brief calculates the rate of change of microstructure
 !---------------------------------------------------------------------------------------------------
 module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
-                                            ph,en,ip,el)
+                                    ph,en,ip,el)
 
   real(pReal), dimension(3,3), intent(in) :: &
     Mp                                                                                              !< MandelStress
@@ -960,7 +946,6 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
     el                                                                                              !< current element number
 
   integer ::  &
-    ns, &                                                                                           !< short notation for the total number of active slip systems
     c, &                                                                                            !< character of dislocation
     t, &                                                                                            !< type of dislocation
     s                                                                                               !< index of my current slip system
@@ -978,7 +963,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
   real(pReal), dimension(param(ph)%sum_N_sl,4) :: &
     v, &                                                                                            !< current dislocation glide velocity
     v0, &
-    gdot                                                                                            !< shear rates
+    dot_gamma                                                                                       !< shear rates
   real(pReal), dimension(param(ph)%sum_N_sl) :: &
     tau, &                                                                                          !< current resolved shear stress
     v_climb                                                                                         !< climb velocity of edge dipoles
@@ -994,14 +979,10 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
     return
   endif
 
-  associate(prm => param(ph), &
+  associate(prm => param(ph), dst => microstructure(ph), dot => dotState(ph), stt => state(ph))
-            dst => microstructure(ph), &
-            dot => dotState(ph), &
-            stt => state(ph))
-  ns = prm%sum_N_sl
 
   tau = 0.0_pReal
-  gdot = 0.0_pReal
+  dot_gamma = 0.0_pReal
 
   rho  = getRho(ph,en)
   rhoSgl = rho(:,sgl)
@@ -1009,14 +990,14 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
   rho0 = getRho0(ph,en)
   my_rhoSgl0 = rho0(:,sgl)
 
-  forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en)
+  forall (s = 1:prm%sum_N_sl, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en)
-  gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4)
+  dot_gamma = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4)
 
 
 
   ! limits for stable dipole height
-  do s = 1,ns
+  do s = 1,prm%sum_N_sl
-    tau(s) = math_tensordot(Mp, prm%Schmid(1:3,1:3,s)) + dst%tau_back(s,en)
+    tau(s) = math_tensordot(Mp, prm%P_sl(1:3,1:3,s)) + dst%tau_back(s,en)
     if (abs(tau(s)) < 1.0e-15_pReal) tau(s) = 1.0e-15_pReal
   enddo
 
@@ -1035,22 +1016,22 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
   ! dislocation multiplication
   rhoDotMultiplication = 0.0_pReal
   isBCC: if (phase_lattice(ph) == 'cI') then
-    forall (s = 1:ns, sum(abs(v(s,1:4))) > 0.0_pReal)
+    forall (s = 1:prm%sum_N_sl, sum(abs(v(s,1:4))) > 0.0_pReal)
-      rhoDotMultiplication(s,1:2) = sum(abs(gdot(s,3:4))) / prm%b_sl(s) &                           ! assuming double-cross-slip of screws to be decisive for multiplication
+      rhoDotMultiplication(s,1:2) = sum(abs(dot_gamma(s,3:4))) / prm%b_sl(s) &                      ! assuming double-cross-slip of screws to be decisive for multiplication
                                   * sqrt(stt%rho_forest(s,en)) / prm%i_sl(s) ! &                    ! mean free path
                                   ! * 2.0_pReal * sum(abs(v(s,3:4))) / sum(abs(v(s,1:4)))           ! ratio of screw to overall velocity determines edge generation
-      rhoDotMultiplication(s,3:4) = sum(abs(gdot(s,3:4))) /prm%b_sl(s) &                            ! assuming double-cross-slip of screws to be decisive for multiplication
+      rhoDotMultiplication(s,3:4) = sum(abs(dot_gamma(s,3:4))) /prm%b_sl(s) &                       ! assuming double-cross-slip of screws to be decisive for multiplication
                                   * sqrt(stt%rho_forest(s,en)) / prm%i_sl(s) ! &                    ! mean free path
                                   ! * 2.0_pReal * sum(abs(v(s,1:2))) / sum(abs(v(s,1:4)))           ! ratio of edge to overall velocity determines screw generation
     endforall
 
   else isBCC
     rhoDotMultiplication(:,1:4) = spread( &
-          (sum(abs(gdot(:,1:2)),2) * prm%f_ed_mult + sum(abs(gdot(:,3:4)),2)) &
+          (sum(abs(dot_gamma(:,1:2)),2) * prm%f_ed_mult + sum(abs(dot_gamma(:,3:4)),2)) &
         * sqrt(stt%rho_forest(:,en)) / prm%i_sl / prm%b_sl, 2, 4)                                   ! eq. 3.26
   endif isBCC
 
-  forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),en)
+  forall (s = 1:prm%sum_N_sl, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),en)
 
 
   !****************************************************************************
@@ -1059,20 +1040,20 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
   ! formation by glide
   do c = 1,2
     rhoDotSingle2DipoleGlide(:,2*c-1) = -2.0_pReal * dUpper(:,c) / prm%b_sl &
-                                                   * (      rhoSgl(:,2*c-1)  * abs(gdot(:,2*c)) &   ! negative mobile --> positive mobile
+                                                   * (      rhoSgl(:,2*c-1)  * abs(dot_gamma(:,2*c)) &   ! negative mobile --> positive mobile
-                                                      +     rhoSgl(:,2*c)    * abs(gdot(:,2*c-1)) & ! positive mobile --> negative mobile
+                                                      +     rhoSgl(:,2*c)    * abs(dot_gamma(:,2*c-1)) & ! positive mobile --> negative mobile
-                                                      + abs(rhoSgl(:,2*c+4)) * abs(gdot(:,2*c-1)))  ! positive mobile --> negative immobile
+                                                      + abs(rhoSgl(:,2*c+4)) * abs(dot_gamma(:,2*c-1)))  ! positive mobile --> negative immobile
 
     rhoDotSingle2DipoleGlide(:,2*c) = -2.0_pReal * dUpper(:,c) / prm%b_sl &
-                                                 * (      rhoSgl(:,2*c-1)  * abs(gdot(:,2*c)) &     ! negative mobile --> positive mobile
+                                                 * (      rhoSgl(:,2*c-1)  * abs(dot_gamma(:,2*c)) &     ! negative mobile --> positive mobile
-                                                    +     rhoSgl(:,2*c)    * abs(gdot(:,2*c-1)) &   ! positive mobile --> negative mobile
+                                                    +     rhoSgl(:,2*c)    * abs(dot_gamma(:,2*c-1)) &   ! positive mobile --> negative mobile
-                                                    + abs(rhoSgl(:,2*c+3)) * abs(gdot(:,2*c)))      ! negative mobile --> positive immobile
+                                                    + abs(rhoSgl(:,2*c+3)) * abs(dot_gamma(:,2*c)))      ! negative mobile --> positive immobile
 
     rhoDotSingle2DipoleGlide(:,2*c+3) = -2.0_pReal * dUpper(:,c) / prm%b_sl &
-                                                   * rhoSgl(:,2*c+3) * abs(gdot(:,2*c))             ! negative mobile --> positive immobile
+                                                   * rhoSgl(:,2*c+3) * abs(dot_gamma(:,2*c))             ! negative mobile --> positive immobile
 
     rhoDotSingle2DipoleGlide(:,2*c+4) = -2.0_pReal * dUpper(:,c) / prm%b_sl &
-                                                   * rhoSgl(:,2*c+4) * abs(gdot(:,2*c-1))           ! positive mobile --> negative immobile
+                                                   * rhoSgl(:,2*c+4) * abs(dot_gamma(:,2*c-1))           ! positive mobile --> negative immobile
 
     rhoDotSingle2DipoleGlide(:,c+8) = abs(rhoDotSingle2DipoleGlide(:,2*c+3)) &
                                     + abs(rhoDotSingle2DipoleGlide(:,2*c+4)) &
@@ -1085,13 +1066,13 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
   rhoDotAthermalAnnihilation = 0.0_pReal
   forall (c=1:2) &
     rhoDotAthermalAnnihilation(:,c+8) = -2.0_pReal * dLower(:,c) / prm%b_sl &
-       * (  2.0_pReal * (rhoSgl(:,2*c-1) * abs(gdot(:,2*c)) + rhoSgl(:,2*c) * abs(gdot(:,2*c-1))) & ! was single hitting single
+       * (  2.0_pReal * (rhoSgl(:,2*c-1) * abs(dot_gamma(:,2*c)) + rhoSgl(:,2*c) * abs(dot_gamma(:,2*c-1))) & ! was single hitting single
-          + 2.0_pReal * (abs(rhoSgl(:,2*c+3)) * abs(gdot(:,2*c)) + abs(rhoSgl(:,2*c+4)) * abs(gdot(:,2*c-1))) & ! was single hitting immobile single or was immobile single hit by single
+          + 2.0_pReal * (abs(rhoSgl(:,2*c+3)) * abs(dot_gamma(:,2*c)) + abs(rhoSgl(:,2*c+4)) * abs(dot_gamma(:,2*c-1))) & ! was single hitting immobile single or was immobile single hit by single
-          + rhoDip(:,c) * (abs(gdot(:,2*c-1)) + abs(gdot(:,2*c))))                                  ! single knocks dipole constituent
+          + rhoDip(:,c) * (abs(dot_gamma(:,2*c-1)) + abs(dot_gamma(:,2*c))))                                  ! single knocks dipole constituent
 
   ! annihilated screw dipoles leave edge jogs behind on the colinear system
   if (phase_lattice(ph) == 'cF') &
-    forall (s = 1:ns, prm%colinearSystem(s) > 0) &
+    forall (s = 1:prm%sum_N_sl, prm%colinearSystem(s) > 0) &
       rhoDotAthermalAnnihilation(prm%colinearSystem(s),1:2) = - rhoDotAthermalAnnihilation(s,10) &
         * 0.25_pReal * sqrt(stt%rho_forest(s,en)) * (dUpper(s,2) + dLower(s,2)) * prm%f_ed
 
@@ -1101,7 +1082,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
   D_SD = prm%D_0 * exp(-prm%Q_cl / (kB * Temperature))                                              ! eq. 3.53
   v_climb = D_SD * prm%mu * prm%V_at &
           / (PI * (1.0_pReal-prm%nu) * (dUpper(:,1) + dLower(:,1)) * kB * Temperature)              ! eq. 3.54
-  forall (s = 1:ns, dUpper(s,1) > dLower(s,1)) &
+  forall (s = 1:prm%sum_N_sl, dUpper(s,1) > dLower(s,1)) &
     rhoDotThermalAnnihilation(s,9) = max(- 4.0_pReal * rhoDip(s,1) * v_climb(s) / (dUpper(s,1) - dLower(s,1)), &
                                          - rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) &
                                                                   - rhoDotSingle2DipoleGlide(s,9))  ! make sure that we do not annihilate more dipoles than we have
@@ -1124,7 +1105,7 @@ module subroutine nonlocal_dotState(Mp, Temperature,timestep, &
     plasticState(ph)%dotState = IEEE_value(1.0_pReal,IEEE_quiet_NaN)
   else
     dot%rho(:,en) = pack(rhoDot,.true.)
-    dot%gamma(:,en) = sum(gdot,2)
+    dot%gamma(:,en) = sum(dot_gamma,2)
   endif
 
   end associate
@@ -1174,7 +1155,7 @@ function rhoDotFlux(timestep,ph,en,ip,el)
     v, &                                                                                            !< current dislocation glide velocity
     v0, &
     neighbor_v0, &                                                                                  !< dislocation glide velocity of enighboring ip
-    gdot                                                                                            !< shear rates
+    dot_gamma                                                                                            !< shear rates
   real(pReal), dimension(3,param(ph)%sum_N_sl,4) :: &
     m                                                                                               !< direction of dislocation motion
   real(pReal), dimension(3,3) :: &
@@ -1200,7 +1181,7 @@ function rhoDotFlux(timestep,ph,en,ip,el)
             stt => state(ph))
   ns = prm%sum_N_sl
 
-  gdot = 0.0_pReal
+  dot_gamma = 0.0_pReal
 
   rho  = getRho(ph,en)
   rhoSgl = rho(:,sgl)
@@ -1208,7 +1189,7 @@ function rhoDotFlux(timestep,ph,en,ip,el)
   my_rhoSgl0 = rho0(:,sgl)
 
   forall (s = 1:ns, t = 1:4) v(s,t) = plasticState(ph)%state(iV(s,t,ph),en)                   !ToDo: MD: I think we should use state0 here
-  gdot = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4)
+  dot_gamma = rhoSgl(:,1:4) * v * spread(prm%b_sl,2,4)
 
   forall (s = 1:ns, t = 1:4) v0(s,t) = plasticState(ph)%state0(iV(s,t,ph),en)
 
@@ -1218,14 +1199,14 @@ function rhoDotFlux(timestep,ph,en,ip,el)
   if (plasticState(ph)%nonlocal) then
 
     !*** check CFL (Courant-Friedrichs-Lewy) condition for flux
-    if (any( abs(gdot) > 0.0_pReal &                                                                ! any active slip system ...
+    if (any( abs(dot_gamma) > 0.0_pReal &                                                                ! any active slip system ...
             .and. prm%C_CFL * abs(v0) * timestep &
                 > IPvolume(ip,el) / maxval(IParea(:,ip,el)))) then                                  ! ...with velocity above critical value (we use the reference volume and area for simplicity here)
 #ifdef DEBUG
     if (debugConstitutive%extensive) then
       print'(a,i5,a,i2)', '<< CONST >> CFL condition not fullfilled at el ',el,' ip ',ip
       print'(a,e10.3,a,e10.3)', '<< CONST >> velocity is at  ', &
-        maxval(abs(v0), abs(gdot) > 0.0_pReal &
+        maxval(abs(v0), abs(dot_gamma) > 0.0_pReal &
                        .and.  prm%C_CFL * abs(v0) * timestep &
                              > IPvolume(ip,el) / maxval(IParea(:,ip,el))), &
         ' at a timestep of ',timestep

src/phase_mechanical_plastic_phenopowerlaw.f90

@@ -33,8 +33,8 @@ submodule(phase:plastic) phenopowerlaw
     real(pReal),               allocatable, dimension(:,:,:) :: &
       P_sl, &
       P_tw, &
-      nonSchmid_pos, &
+      P_nS_pos, &
-      nonSchmid_neg
+      P_nS_neg
     integer :: &
       sum_N_sl, &                                                                                   !< total number of active slip system
       sum_N_tw                                                                                      !< total number of active twin systems
@@ -46,10 +46,10 @@ submodule(phase:plastic) phenopowerlaw
 
   type :: tPhenopowerlawState
     real(pReal), pointer, dimension(:,:) :: &
-      xi_slip, &
+      xi_sl, &
-      xi_twin, &
+      xi_tw, &
-      gamma_slip, &
+      gamma_sl, &
-      gamma_twin
+      gamma_tw
   end type tPhenopowerlawState
 
 !--------------------------------------------------------------------------------------------------
@@ -115,21 +115,19 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
     N_sl         = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray)
     prm%sum_N_sl = sum(abs(N_sl))
     slipActive: if (prm%sum_N_sl > 0) then
-      prm%P_sl = lattice_SchmidMatrix_slip(N_sl,phase%get_asString('lattice'),&
+      prm%P_sl = lattice_SchmidMatrix_slip(N_sl,phase_lattice(ph),phase_cOverA(ph))
-                                           phase%get_asFloat('c/a',defaultVal=0.0_pReal))
 
-      if(phase%get_asString('lattice') == 'cI') then
+      if (phase_lattice(ph) == 'cI') then
         a = pl%get_as1dFloat('a_nonSchmid',defaultVal=emptyRealArray)
         if(size(a) > 0) prm%nonSchmidActive = .true.
-        prm%nonSchmid_pos  = lattice_nonSchmidMatrix(N_sl,a,+1)
+        prm%P_nS_pos = lattice_nonSchmidMatrix(N_sl,a,+1)
-        prm%nonSchmid_neg  = lattice_nonSchmidMatrix(N_sl,a,-1)
+        prm%P_nS_neg = lattice_nonSchmidMatrix(N_sl,a,-1)
       else
-        prm%nonSchmid_pos  = prm%P_sl
+        prm%P_nS_pos = prm%P_sl
-        prm%nonSchmid_neg  = prm%P_sl
+        prm%P_nS_neg = prm%P_sl
       endif
-      prm%h_sl_sl   = lattice_interaction_SlipBySlip(N_sl, &
+      prm%h_sl_sl   = lattice_interaction_SlipBySlip(N_sl,pl%get_as1dFloat('h_sl-sl'), &
-                                                     pl%get_as1dFloat('h_sl-sl'), &
+                                                     phase_lattice(ph))
-                                                     phase%get_asString('lattice'))
 
       xi_0_sl             = pl%get_as1dFloat('xi_0_sl',   requiredSize=size(N_sl))
       prm%xi_inf_sl       = pl%get_as1dFloat('xi_inf_sl', requiredSize=size(N_sl))
@@ -164,13 +162,11 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
     N_tw         = pl%get_as1dInt('N_tw', defaultVal=emptyIntArray)
     prm%sum_N_tw = sum(abs(N_tw))
     twinActive: if (prm%sum_N_tw > 0) then
-      prm%P_tw            = lattice_SchmidMatrix_twin(N_tw,phase%get_asString('lattice'),&
+      prm%P_tw            = lattice_SchmidMatrix_twin(N_tw,phase_lattice(ph),phase_cOverA(ph))
-                                                      phase%get_asFloat('c/a',defaultVal=0.0_pReal))
+      prm%h_tw_tw         = lattice_interaction_TwinByTwin(N_tw,pl%get_as1dFloat('h_tw-tw'), &
-      prm%h_tw_tw         = lattice_interaction_TwinByTwin(N_tw,&
+                                                           phase_lattice(ph))
-                                                           pl%get_as1dFloat('h_tw-tw'), &
+      prm%gamma_char      = lattice_characteristicShear_twin(N_tw,phase_lattice(ph),&
-                                                           phase%get_asString('lattice'))
+                                                             phase_cOverA(ph))
-      prm%gamma_char      = lattice_characteristicShear_twin(N_tw,phase%get_asString('lattice'),&
-                                                             phase%get_asFloat('c/a',defaultVal=0.0_pReal))
 
       xi_0_tw             = pl%get_as1dFloat('xi_0_tw',requiredSize=size(N_tw))
 
@@ -200,12 +196,10 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
 ! slip-twin related parameters
     slipAndTwinActive: if (prm%sum_N_sl > 0 .and. prm%sum_N_tw > 0) then
       prm%h_0_tw_sl  = pl%get_asFloat('h_0_tw-sl')
-      prm%h_sl_tw    = lattice_interaction_SlipByTwin(N_sl,N_tw,&
+      prm%h_sl_tw    = lattice_interaction_SlipByTwin(N_sl,N_tw,pl%get_as1dFloat('h_sl-tw'), &
-                                                      pl%get_as1dFloat('h_sl-tw'), &
+                                                      phase_lattice(ph))
-                                                      phase%get_asString('lattice'))
+      prm%h_tw_sl    = lattice_interaction_TwinBySlip(N_tw,N_sl,pl%get_as1dFloat('h_tw-sl'), &
-      prm%h_tw_sl    = lattice_interaction_TwinBySlip(N_tw,N_sl,&
+                                                      phase_lattice(ph))
-                                                      pl%get_as1dFloat('h_tw-sl'), &
-                                                      phase%get_asString('lattice'))
     else slipAndTwinActive
       allocate(prm%h_sl_tw(prm%sum_N_sl,prm%sum_N_tw))                                              ! at least one dimension is 0
       allocate(prm%h_tw_sl(prm%sum_N_tw,prm%sum_N_sl))                                              ! at least one dimension is 0
@@ -235,30 +229,30 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
 ! state aliases and initialization
     startIndex = 1
     endIndex   = prm%sum_N_sl
-    stt%xi_slip => plasticState(ph)%state   (startIndex:endIndex,:)
+    stt%xi_sl => plasticState(ph)%state   (startIndex:endIndex,:)
-    stt%xi_slip =  spread(xi_0_sl, 2, Nmembers)
+    stt%xi_sl =  spread(xi_0_sl, 2, Nmembers)
-    dot%xi_slip => plasticState(ph)%dotState(startIndex:endIndex,:)
+    dot%xi_sl => plasticState(ph)%dotState(startIndex:endIndex,:)
     plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
     if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
 
     startIndex = endIndex + 1
     endIndex   = endIndex + prm%sum_N_tw
-    stt%xi_twin => plasticState(ph)%state   (startIndex:endIndex,:)
+    stt%xi_tw => plasticState(ph)%state   (startIndex:endIndex,:)
-    stt%xi_twin =  spread(xi_0_tw, 2, Nmembers)
+    stt%xi_tw =  spread(xi_0_tw, 2, Nmembers)
-    dot%xi_twin => plasticState(ph)%dotState(startIndex:endIndex,:)
+    dot%xi_tw => plasticState(ph)%dotState(startIndex:endIndex,:)
     plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
 
     startIndex = endIndex + 1
     endIndex   = endIndex + prm%sum_N_sl
-    stt%gamma_slip => plasticState(ph)%state   (startIndex:endIndex,:)
+    stt%gamma_sl => plasticState(ph)%state   (startIndex:endIndex,:)
-    dot%gamma_slip => plasticState(ph)%dotState(startIndex:endIndex,:)
+    dot%gamma_sl => plasticState(ph)%dotState(startIndex:endIndex,:)
     plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
     if(any(plasticState(ph)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_gamma'
 
     startIndex = endIndex + 1
     endIndex   = endIndex + prm%sum_N_tw
-    stt%gamma_twin => plasticState(ph)%state   (startIndex:endIndex,:)
+    stt%gamma_tw => plasticState(ph)%state   (startIndex:endIndex,:)
-    dot%gamma_twin => plasticState(ph)%dotState(startIndex:endIndex,:)
+    dot%gamma_tw => plasticState(ph)%dotState(startIndex:endIndex,:)
     plasticState(ph)%atol(startIndex:endIndex) = pl%get_asFloat('atol_gamma',defaultVal=1.0e-6_pReal)
 
     end associate
@@ -293,31 +287,31 @@ pure module subroutine phenopowerlaw_LpAndItsTangent(Lp,dLp_dMp,Mp,ph,en)
   integer :: &
     i,k,l,m,n
   real(pReal), dimension(param(ph)%sum_N_sl) :: &
-    gdot_slip_pos,gdot_slip_neg, &
+    dot_gamma_sl_pos,dot_gamma_sl_neg, &
-    dgdot_dtauslip_pos,dgdot_dtauslip_neg
+    ddot_gamma_dtau_sl_pos,ddot_gamma_dtau_sl_neg
   real(pReal), dimension(param(ph)%sum_N_tw) :: &
-    gdot_twin,dgdot_dtautwin
+    dot_gamma_tw,ddot_gamma_dtau_tw
 
   Lp = 0.0_pReal
   dLp_dMp = 0.0_pReal
 
   associate(prm => param(ph))
 
-  call kinetics_slip(Mp,ph,en,gdot_slip_pos,gdot_slip_neg,dgdot_dtauslip_pos,dgdot_dtauslip_neg)
+  call kinetics_sl(Mp,ph,en,dot_gamma_sl_pos,dot_gamma_sl_neg,ddot_gamma_dtau_sl_pos,ddot_gamma_dtau_sl_neg)
   slipSystems: do i = 1, prm%sum_N_sl
-    Lp = Lp + (gdot_slip_pos(i)+gdot_slip_neg(i))*prm%P_sl(1:3,1:3,i)
+    Lp = Lp + (dot_gamma_sl_pos(i)+dot_gamma_sl_neg(i))*prm%P_sl(1:3,1:3,i)
     forall (k=1:3,l=1:3,m=1:3,n=1:3) &
       dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
-                       + dgdot_dtauslip_pos(i) * prm%P_sl(k,l,i) * prm%nonSchmid_pos(m,n,i) &
+                       + ddot_gamma_dtau_sl_pos(i) * prm%P_sl(k,l,i) * prm%P_nS_pos(m,n,i) &
-                       + dgdot_dtauslip_neg(i) * prm%P_sl(k,l,i) * prm%nonSchmid_neg(m,n,i)
+                       + ddot_gamma_dtau_sl_neg(i) * prm%P_sl(k,l,i) * prm%P_nS_neg(m,n,i)
   enddo slipSystems
 
-  call kinetics_twin(Mp,ph,en,gdot_twin,dgdot_dtautwin)
+  call kinetics_tw(Mp,ph,en,dot_gamma_tw,ddot_gamma_dtau_tw)
   twinSystems: do i = 1, prm%sum_N_tw
-    Lp = Lp + gdot_twin(i)*prm%P_tw(1:3,1:3,i)
+    Lp = Lp + dot_gamma_tw(i)*prm%P_tw(1:3,1:3,i)
     forall (k=1:3,l=1:3,m=1:3,n=1:3) &
       dLp_dMp(k,l,m,n) = dLp_dMp(k,l,m,n) &
-                       + dgdot_dtautwin(i)*prm%P_tw(k,l,i)*prm%P_tw(m,n,i)
+                       + ddot_gamma_dtau_tw(i)*prm%P_tw(k,l,i)*prm%P_tw(m,n,i)
   enddo twinSystems
 
   end associate
@@ -337,46 +331,32 @@ module subroutine phenopowerlaw_dotState(Mp,ph,en)
     en
 
   real(pReal) :: &
-    c_SlipSlip,c_TwinSlip,c_TwinTwin, &
+    xi_sl_sat_offset,&
-    xi_slip_sat_offset,&
+    sumF
-    sumGamma,sumF
   real(pReal), dimension(param(ph)%sum_N_sl) :: &
-    left_SlipSlip,right_SlipSlip, &
+    dot_gamma_sl_pos,dot_gamma_sl_neg, &
-    gdot_slip_pos,gdot_slip_neg
+    right_SlipSlip
 
-  associate(prm => param(ph), stt => state(ph), &
-  dot => dotState(ph))
 
-  sumGamma = sum(stt%gamma_slip(:,en))
+  associate(prm => param(ph), stt => state(ph), dot => dotState(ph))
-  sumF     = sum(stt%gamma_twin(:,en)/prm%gamma_char)
 
-!--------------------------------------------------------------------------------------------------
+    call kinetics_sl(Mp,ph,en,dot_gamma_sl_pos,dot_gamma_sl_neg)
-! system-independent (nonlinear) prefactors to M_Xx (X influenced by x) matrices
+    dot%gamma_sl(:,en) = abs(dot_gamma_sl_pos+dot_gamma_sl_neg)
-  c_SlipSlip = prm%h_0_sl_sl * (1.0_pReal + prm%c_1*sumF** prm%c_2)
+    call kinetics_tw(Mp,ph,en,dot%gamma_tw(:,en))
-  c_TwinSlip = prm%h_0_tw_sl * sumGamma**prm%c_3
-  c_TwinTwin = prm%h_0_tw_tw * sumF**prm%c_4
 
-!--------------------------------------------------------------------------------------------------
+    sumF = sum(stt%gamma_tw(:,en)/prm%gamma_char)
-!  calculate left and right vectors
+    xi_sl_sat_offset = prm%f_sat_sl_tw*sqrt(sumF)
-  left_SlipSlip  = 1.0_pReal + prm%h_int
+    right_SlipSlip = sign(abs(1.0_pReal-stt%xi_sl(:,en) / (prm%xi_inf_sl+xi_sl_sat_offset))**prm%a_sl, &
-  xi_slip_sat_offset = prm%f_sat_sl_tw*sqrt(sumF)
+                          1.0_pReal-stt%xi_sl(:,en) / (prm%xi_inf_sl+xi_sl_sat_offset))
-  right_SlipSlip = sign(abs(1.0_pReal-stt%xi_slip(:,en) / (prm%xi_inf_sl+xi_slip_sat_offset)) **prm%a_sl, &
-                        1.0_pReal-stt%xi_slip(:,en) / (prm%xi_inf_sl+xi_slip_sat_offset))
 
-!--------------------------------------------------------------------------------------------------
+    dot%xi_sl(:,en) = prm%h_0_sl_sl * (1.0_pReal + prm%c_1*sumF** prm%c_2) * (1.0_pReal + prm%h_int) &
-! shear rates
+                      * matmul(prm%h_sl_sl,dot%gamma_sl(:,en)*right_SlipSlip) &
-  call kinetics_slip(Mp,ph,en,gdot_slip_pos,gdot_slip_neg)
+                    + matmul(prm%h_sl_tw,dot%gamma_tw(:,en))
-  dot%gamma_slip(:,en) = abs(gdot_slip_pos+gdot_slip_neg)
-  call kinetics_twin(Mp,ph,en,dot%gamma_twin(:,en))
 
-!--------------------------------------------------------------------------------------------------
+    dot%xi_tw(:,en) = prm%h_0_tw_sl * sum(stt%gamma_sl(:,en))**prm%c_3 &
-! hardening
+                      * matmul(prm%h_tw_sl,dot%gamma_sl(:,en)) &
-  dot%xi_slip(:,en) = c_SlipSlip * left_SlipSlip * &
+                    + prm%h_0_tw_tw * sumF**prm%c_4 * matmul(prm%h_tw_tw,dot%gamma_tw(:,en))
-                      matmul(prm%h_sl_sl,dot%gamma_slip(:,en)*right_SlipSlip) &
-                    + matmul(prm%h_sl_tw,dot%gamma_twin(:,en))
 
-  dot%xi_twin(:,en) = c_TwinSlip * matmul(prm%h_tw_sl,dot%gamma_slip(:,en)) &
-                    + c_TwinTwin * matmul(prm%h_tw_tw,dot%gamma_twin(:,en))
   end associate
 
 end subroutine phenopowerlaw_dotState
@@ -397,17 +377,17 @@ module subroutine plastic_phenopowerlaw_results(ph,group)
     select case(trim(prm%output(o)))
 
       case('xi_sl')
-        if(prm%sum_N_sl>0) call results_writeDataset(stt%xi_slip,group,trim(prm%output(o)), &
+        if(prm%sum_N_sl>0) call results_writeDataset(stt%xi_sl,group,trim(prm%output(o)), &
                                                      'resistance against plastic slip','Pa')
       case('gamma_sl')
-        if(prm%sum_N_sl>0) call results_writeDataset(stt%gamma_slip,group,trim(prm%output(o)), &
+        if(prm%sum_N_sl>0) call results_writeDataset(stt%gamma_sl,group,trim(prm%output(o)), &
                                                      'plastic shear','1')
 
       case('xi_tw')
-        if(prm%sum_N_tw>0) call results_writeDataset(stt%xi_twin,group,trim(prm%output(o)), &
+        if(prm%sum_N_tw>0) call results_writeDataset(stt%xi_tw,group,trim(prm%output(o)), &
                                                      'resistance against twinning','Pa')
       case('gamma_tw')
-        if(prm%sum_N_tw>0) call results_writeDataset(stt%gamma_twin,group,trim(prm%output(o)), &
+        if(prm%sum_N_tw>0) call results_writeDataset(stt%gamma_tw,group,trim(prm%output(o)), &
                                                      'twinning shear','1')
 
     end select
@@ -424,8 +404,8 @@ end subroutine plastic_phenopowerlaw_results
 ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
 ! have the optional arguments at the end.
 !--------------------------------------------------------------------------------------------------
-pure subroutine kinetics_slip(Mp,ph,en, &
+pure subroutine kinetics_sl(Mp,ph,en, &
-                              gdot_slip_pos,gdot_slip_neg,dgdot_dtau_slip_pos,dgdot_dtau_slip_neg)
+                            dot_gamma_sl_pos,dot_gamma_sl_neg,ddot_gamma_dtau_sl_pos,ddot_gamma_dtau_sl_neg)
 
   real(pReal), dimension(3,3),  intent(in) :: &
     Mp                                                                                              !< Mandel stress
@@ -434,56 +414,57 @@ pure subroutine kinetics_slip(Mp,ph,en, &
     en
 
   real(pReal),                  intent(out), dimension(param(ph)%sum_N_sl) :: &
-    gdot_slip_pos, &
+    dot_gamma_sl_pos, &
-    gdot_slip_neg
+    dot_gamma_sl_neg
   real(pReal),                  intent(out), optional, dimension(param(ph)%sum_N_sl) :: &
-    dgdot_dtau_slip_pos, &
+    ddot_gamma_dtau_sl_pos, &
-    dgdot_dtau_slip_neg
+    ddot_gamma_dtau_sl_neg
 
   real(pReal), dimension(param(ph)%sum_N_sl) :: &
-    tau_slip_pos, &
+    tau_sl_pos, &
-    tau_slip_neg
+    tau_sl_neg
   integer :: i
 
   associate(prm => param(ph), stt => state(ph))
 
-  do i = 1, prm%sum_N_sl
+    do i = 1, prm%sum_N_sl
-    tau_slip_pos(i) =       math_tensordot(Mp,prm%nonSchmid_pos(1:3,1:3,i))
+      tau_sl_pos(i) =       math_tensordot(Mp,prm%P_nS_pos(1:3,1:3,i))
-    tau_slip_neg(i) = merge(math_tensordot(Mp,prm%nonSchmid_neg(1:3,1:3,i)), &
+      tau_sl_neg(i) = merge(math_tensordot(Mp,prm%P_nS_neg(1:3,1:3,i)), &
                             0.0_pReal, prm%nonSchmidActive)
-  enddo
+    enddo
 
-  where(dNeq0(tau_slip_pos))
+    where(dNeq0(tau_sl_pos))
-    gdot_slip_pos = prm%dot_gamma_0_sl * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) &          ! 1/2 if non-Schmid active
+      dot_gamma_sl_pos = prm%dot_gamma_0_sl * merge(0.5_pReal,1.0_pReal, prm%nonSchmidActive) &     ! 1/2 if non-Schmid active
-                  * sign(abs(tau_slip_pos/stt%xi_slip(:,en))**prm%n_sl,  tau_slip_pos)
+                       * sign(abs(tau_sl_pos/stt%xi_sl(:,en))**prm%n_sl,  tau_sl_pos)
-  else where
-    gdot_slip_pos = 0.0_pReal
-  end where
-
-  where(dNeq0(tau_slip_neg))
-    gdot_slip_neg = prm%dot_gamma_0_sl * 0.5_pReal &                                                ! only used if non-Schmid active, always 1/2
-                  * sign(abs(tau_slip_neg/stt%xi_slip(:,en))**prm%n_sl,  tau_slip_neg)
-  else where
-    gdot_slip_neg = 0.0_pReal
-  end where
-
-  if (present(dgdot_dtau_slip_pos)) then
-    where(dNeq0(gdot_slip_pos))
-      dgdot_dtau_slip_pos = gdot_slip_pos*prm%n_sl/tau_slip_pos
     else where
-      dgdot_dtau_slip_pos = 0.0_pReal
+      dot_gamma_sl_pos = 0.0_pReal
     end where
-  endif
+
-  if (present(dgdot_dtau_slip_neg)) then
+    where(dNeq0(tau_sl_neg))
-    where(dNeq0(gdot_slip_neg))
+      dot_gamma_sl_neg = prm%dot_gamma_0_sl * 0.5_pReal &                                           ! only used if non-Schmid active, always 1/2
-      dgdot_dtau_slip_neg = gdot_slip_neg*prm%n_sl/tau_slip_neg
+                       * sign(abs(tau_sl_neg/stt%xi_sl(:,en))**prm%n_sl,  tau_sl_neg)
     else where
-      dgdot_dtau_slip_neg = 0.0_pReal
+      dot_gamma_sl_neg = 0.0_pReal
     end where
-  endif
+
+    if (present(ddot_gamma_dtau_sl_pos)) then
+      where(dNeq0(dot_gamma_sl_pos))
+        ddot_gamma_dtau_sl_pos = dot_gamma_sl_pos*prm%n_sl/tau_sl_pos
+      else where
+        ddot_gamma_dtau_sl_pos = 0.0_pReal
+      end where
+    endif
+    if (present(ddot_gamma_dtau_sl_neg)) then
+      where(dNeq0(dot_gamma_sl_neg))
+        ddot_gamma_dtau_sl_neg = dot_gamma_sl_neg*prm%n_sl/tau_sl_neg
+      else where
+        ddot_gamma_dtau_sl_neg = 0.0_pReal
+      end where
+    endif
+
   end associate
 
-end subroutine kinetics_slip
+end subroutine kinetics_sl
 
 
 !--------------------------------------------------------------------------------------------------
@@ -493,8 +474,8 @@ end subroutine kinetics_slip
 ! NOTE: Against the common convention, the result (i.e. intent(out)) variables are the last to
 ! have the optional arguments at the end.
 !--------------------------------------------------------------------------------------------------
-pure subroutine kinetics_twin(Mp,ph,en,&
+pure subroutine kinetics_tw(Mp,ph,en,&
-                              gdot_twin,dgdot_dtau_twin)
+                            dot_gamma_tw,ddot_gamma_dtau_tw)
 
   real(pReal), dimension(3,3),  intent(in) :: &
     Mp                                                                                              !< Mandel stress
@@ -503,37 +484,36 @@ pure subroutine kinetics_twin(Mp,ph,en,&
     en
 
   real(pReal), dimension(param(ph)%sum_N_tw), intent(out) :: &
-    gdot_twin
+    dot_gamma_tw
   real(pReal), dimension(param(ph)%sum_N_tw), intent(out), optional :: &
-    dgdot_dtau_twin
+    ddot_gamma_dtau_tw
 
   real(pReal), dimension(param(ph)%sum_N_tw) :: &
-    tau_twin
+    tau_tw
   integer :: i
 
+
   associate(prm => param(ph), stt => state(ph))
 
-  do i = 1, prm%sum_N_tw
+    tau_tw = [(math_tensordot(Mp,prm%P_tw(1:3,1:3,i)),i=1,prm%sum_N_tw)]
-    tau_twin(i)  = math_tensordot(Mp,prm%P_tw(1:3,1:3,i))
-  enddo
 
-  where(tau_twin > 0.0_pReal)
+    where(tau_tw > 0.0_pReal)
-    gdot_twin = (1.0_pReal-sum(stt%gamma_twin(:,en)/prm%gamma_char)) &                              ! only twin in untwinned volume fraction
+      dot_gamma_tw = (1.0_pReal-sum(stt%gamma_tw(:,en)/prm%gamma_char)) &                           ! only twin in untwinned volume fraction
-              * prm%dot_gamma_0_tw*(abs(tau_twin)/stt%xi_twin(:,en))**prm%n_tw
+                   * prm%dot_gamma_0_tw*(abs(tau_tw)/stt%xi_tw(:,en))**prm%n_tw
-  else where
-    gdot_twin = 0.0_pReal
-  end where
-
-  if (present(dgdot_dtau_twin)) then
-    where(dNeq0(gdot_twin))
-      dgdot_dtau_twin = gdot_twin*prm%n_tw/tau_twin
     else where
-      dgdot_dtau_twin = 0.0_pReal
+      dot_gamma_tw = 0.0_pReal
     end where
-  endif
+
+    if (present(ddot_gamma_dtau_tw)) then
+      where(dNeq0(dot_gamma_tw))
+        ddot_gamma_dtau_tw = dot_gamma_tw*prm%n_tw/tau_tw
+      else where
+        ddot_gamma_dtau_tw = 0.0_pReal
+      end where
+    endif
 
   end associate
 
-end subroutine kinetics_twin
+end subroutine kinetics_tw
 
 end submodule phenopowerlaw

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
@@ -216,7 +216,7 @@ module function phase_K_T(co,ce) result(K)
 end function phase_K_T
 
 
-module function thermal_stress(Delta_t,ph,en) result(converged_)           ! ?? why is this called "stress" when it seems closer to "updateState" ??
+module function phase_thermal_constitutive(Delta_t,ph,en) result(converged_)
 
   real(pReal), intent(in) :: Delta_t
   integer, intent(in) :: ph, en
@@ -225,7 +225,7 @@ module function thermal_stress(Delta_t,ph,en) result(converged_)           ! ??
 
   converged_ = .not. integrateThermalState(Delta_t,ph,en)
 
-end function thermal_stress
+end function phase_thermal_constitutive
 
 
 !--------------------------------------------------------------------------------------------------

src/prec.f90

@@ -38,16 +38,14 @@ module prec
       sizeDotState     = 0, &                                                                       !< size of dot state, i.e. state(1:sizeDot) follows time evolution by dotState rates
       offsetDeltaState = 0, &                                                                       !< index offset of delta state
       sizeDeltaState   = 0                                                                          !< size of delta state, i.e. state(offset+1:offset+sizeDelta) follows time evolution by deltaState increments
-    ! http://stackoverflow.com/questions/3948210
+    real(pReal), allocatable, dimension(:) :: &
-    real(pReal), pointer,     dimension(:), contiguous :: &
       atol
-    real(pReal), pointer,     dimension(:,:), contiguous :: &                                       ! a pointer is needed here because we might point to state/doState. However, they will never point to something, but are rather allocated and, hence, contiguous
+    ! http://stackoverflow.com/questions/3948210
+    real(pReal), pointer,     dimension(:,:), contiguous :: &                                       !< is basically an allocatable+target, but in a type needs to be pointer
       state0, &
       state, &                                                                                      !< state
       dotState, &                                                                                   !< rate of state change
       deltaState                                                                                    !< increment of state change
-    real(pReal), allocatable, dimension(:,:) :: &
-      subState0
   end type
 
   type, extends(tState) :: tPlasticState
@@ -61,12 +59,9 @@ module prec
   real(pReal), private, parameter :: PREAL_EPSILON = epsilon(0.0_pReal)                             !< minimum positive number such that 1.0 + EPSILON /= 1.0.
   real(pReal), private, parameter :: PREAL_MIN     = tiny(0.0_pReal)                                !< smallest normalized floating point number
 
-  integer,                   dimension(0), parameter :: &
+  integer,                   dimension(0), parameter :: emptyIntArray    = [integer::]
-    emptyIntArray    = [integer::]
+  real(pReal),               dimension(0), parameter :: emptyRealArray   = [real(pReal)::]
-  real(pReal),               dimension(0), parameter :: &
+  character(len=pStringLen), dimension(0), parameter :: emptyStringArray = [character(len=pStringLen)::]
-    emptyRealArray   = [real(pReal)::]
-  character(len=pStringLen), dimension(0), parameter :: &
-    emptyStringArray = [character(len=pStringLen)::]
 
   private :: &
     selfTest