Merge branch 'Fortran-cleaning' into 'development'

Fortran cleaning

See merge request damask/DAMASK!390

examples/mesh/numerics.config

@@ -1,3 +0,0 @@
-residualStiffness 0.001
-charLength        0.02
-petsc_options -mech_snes_type newtonls -mech_ksp_type fgmres -mech_pc_type ml -mech_ksp_monitor

src/CPFEM.f90

@@ -177,11 +177,6 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
 
   if (iand(mode, CPFEM_AGERESULTS) /= 0_pInt) call CPFEM_forward
 
-    !chosenThermal1: select case (thermal_type(material_homogenizationAt(elCP)))
-    !  case (THERMAL_conduction_ID) chosenThermal1
-    !    temperature(material_homogenizationAt(elCP))%p(material_homogenizationMemberAt(ip,elCP)) = &
-    !      temperature_inp
-    !end select chosenThermal1
     homogenization_F0(1:3,1:3,ce) = ffn
     homogenization_F(1:3,1:3,ce) = ffn1
 

src/IO.f90

@@ -20,6 +20,9 @@ module IO
   character, parameter, public :: &
     IO_EOL = new_line('DAMASK'), &                                                                  !< end of line character
     IO_COMMENT = '#'
+  character, parameter :: &
+    CR = achar(13), &
+    LF = IO_EOL
   character(len=*), parameter :: &
     IO_DIVIDER = '───────────────────'//&
                  '───────────────────'//&
@@ -112,8 +115,8 @@ end function IO_readlines
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief Read whole file.
+!> @brief Read ASCII file.
-!> @details ensures that the string ends with a new line (expected UNIX behavior)
+!> @details Proper Unix style (LF line endings and LF at EOF) is ensured.
 !--------------------------------------------------------------------------------------------------
 function IO_read(fileName) result(fileContent)
 
@@ -124,7 +127,6 @@ function IO_read(fileName) result(fileContent)
     fileLength, &
     fileUnit, &
     myStat
-  character, parameter :: CR = achar(13)
 
 
   inquire(file = fileName, size=fileLength)
@@ -141,15 +143,7 @@ function IO_read(fileName) result(fileContent)
   if(myStat /= 0) call IO_error(102,ext_msg=trim(fileName))
   close(fileUnit)
 
-  foundCRLF: if (scan(fileContent(:index(fileContent,IO_EOL)),CR) /= 0) then
+  if (scan(fileContent(:index(fileContent,LF)),CR//LF) /= 0) fileContent = CRLF2LF(fileContent)
-    CRLF2LF: block
-      integer :: c
-      do c=1, len(fileContent)
-        if (fileContent(c:c) == CR) fileContent(c:c) = ' '
-      enddo
-    end block CRLF2LF
-  endif foundCRLF
-
   if(fileContent(fileLength:fileLength) /= IO_EOL)           fileContent = fileContent//IO_EOL      ! ensure EOL@EOF
 
 end function IO_read
@@ -633,10 +627,36 @@ subroutine IO_warning(warning_ID,el,ip,g,ext_msg)
 end subroutine IO_warning
 
 
+!--------------------------------------------------------------------------------------------------
+!> @brief Convert Windows (CRLF) to Unix (LF) line endings.
+!--------------------------------------------------------------------------------------------------
+pure function CRLF2LF(string)
+
+  character(len=*), intent(in)  :: string
+  character(len=:), allocatable :: CRLF2LF
+
+  integer :: c,n
+
+
+  allocate(character(len=len_trim(string))::CRLF2LF)
+  if (len(CRLF2LF) == 0) return
+
+  n = 0
+  do c=1, len_trim(string)
+    CRLF2LF(c-n:c-n) = string(c:c)
+    if (c == len_trim(string)) exit
+    if (string(c:c+1) == CR//LF) n = n + 1
+  enddo
+
+  CRLF2LF = CRLF2LF(:c-n)
+
+end function
+
+
 !--------------------------------------------------------------------------------------------------
 !> @brief Check correctness of some IO functions.
 !--------------------------------------------------------------------------------------------------
-subroutine selfTest
+subroutine selfTest()
 
   integer, dimension(:), allocatable :: chunkPos
   character(len=:),      allocatable :: str
@@ -671,6 +691,15 @@ subroutine selfTest
   chunkPos = IO_stringPos(str)
   if(3112019 /= IO_intValue(str,chunkPos,2))        error stop 'IO_intValue'
 
+  if (CRLF2LF('') /= '')                            error stop 'CRLF2LF/0'
+  if (CRLF2LF(LF)     /= LF)                        error stop 'CRLF2LF/1a'
+  if (CRLF2LF(CR//LF) /= LF)                        error stop 'CRLF2LF/1b'
+  if (CRLF2LF(' '//LF)     /= ' '//LF)              error stop 'CRLF2LF/2a'
+  if (CRLF2LF(' '//CR//LF) /= ' '//LF)              error stop 'CRLF2LF/2b'
+  if (CRLF2LF('A'//CR//LF//'B') /= 'A'//LF//'B')    error stop 'CRLF2LF/3'
+  if (CRLF2LF('A'//CR//LF//'B'//CR//LF) /= &
+              'A'//LF//'B'//LF)                     error stop 'CRLF2LF/4'
+
   if(.not. IO_isBlank('  '))                        error stop 'IO_isBlank/1'
   if(.not. IO_isBlank('  #isBlank'))                error stop 'IO_isBlank/2'
   if(      IO_isBlank('  i#s'))                     error stop 'IO_isBlank/3'

src/config.f90

@@ -92,7 +92,7 @@ end subroutine parse_debug
 
 !--------------------------------------------------------------------------------------------------
 !> @brief Deallocate config_material.
-!ToDo: deallocation of numerics debug (optional)
+!ToDo: deallocation of numerics and debug (optional)
 !--------------------------------------------------------------------------------------------------
 subroutine config_deallocate
 

src/discretization.f90

@@ -15,7 +15,7 @@ module discretization
     discretization_Nelems
 
   integer,     public, protected, dimension(:),   allocatable :: &
-    discretization_materialAt   
+    discretization_materialAt                                                                       !ToDo: discretization_ID_material
 
   real(pReal), public, protected, dimension(:,:), allocatable :: &
     discretization_IPcoords0, &

src/grid/DAMASK_grid.f90

@@ -284,7 +284,7 @@ program DAMASK_grid
       if (loadCases(l)%f_restart < 1)  errorID = 839
 
       if (dEq(loadCases(l)%r,1.0_pReal,1.e-9_pReal)) then
-        print'(a)', '  r: 1 (constant step widths)'
+        print'(a)', '  r: 1 (constant step width)'
       else
         print'(a,f0.3)', '  r: ', loadCases(l)%r
       endif

src/homogenization.f90

@@ -236,20 +236,20 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
     NiterationMPstate, &
     ip, &                                                                                            !< integration point number
     el, &                                                                                            !< element number
-    myNgrains, co, ce, ho, en, ph
+    co, ce, ho, en, ph
   logical :: &
     converged
   logical, dimension(2) :: &
     doneAndHappy
 
   !$OMP PARALLEL
-  !$OMP DO PRIVATE(ce,en,ho,myNgrains,NiterationMPstate,converged,doneAndHappy)
+  !$OMP DO PRIVATE(ce,en,ho,NiterationMPstate,converged,doneAndHappy)
   do el = FEsolving_execElem(1),FEsolving_execElem(2)
-    ho = material_homogenizationAt(el)
+
-    myNgrains = homogenization_Nconstituents(ho)
     do ip = FEsolving_execIP(1),FEsolving_execIP(2)
       ce = (el-1)*discretization_nIPs + ip
       en = material_homogenizationEntry(ce)
+      ho = material_homogenizationID(ce)
 
       call phase_restore(ce,.false.) ! wrong name (is more a forward function)
 
@@ -266,7 +266,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
 
         call mechanical_partition(homogenization_F(1:3,1:3,ce),ce)
         converged = .true.
-        do co = 1, myNgrains
+        do co = 1, homogenization_Nconstituents(ho)
           converged = converged .and. crystallite_stress(dt,co,ip,el)
         enddo
 
@@ -290,12 +290,12 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
     !$OMP DO PRIVATE(ho,ph,ce)
     do el = FEsolving_execElem(1),FEsolving_execElem(2)
       if (terminallyIll) continue
-      ho = material_homogenizationAt(el)
       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)
-          ph = material_phaseAt(co,el)
+          ph = material_phaseID(co,ce)
           if (.not. thermal_stress(dt,ph,material_phaseMemberAt(co,ip,el))) then
             if (.not. terminallyIll) &                                                              ! so first signals terminally ill...
               print*, ' Integration point ', ip,' at element ', el, ' terminally ill'
@@ -308,9 +308,9 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
 
     !$OMP DO PRIVATE(ho,ce)
     elementLooping3: do el = FEsolving_execElem(1),FEsolving_execElem(2)
-      ho = material_homogenizationAt(el)
       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

src/homogenization_damage.f90

@@ -39,7 +39,7 @@ module subroutine damage_init()
     configHomogenization, &
     configHomogenizationDamage, &
     num_generic
-  integer :: ho,Nmaterialpoints
+  integer :: ho,Nmembers
 
 
   print'(/,a)', ' <<<+-  homogenization:damage init  -+>>>'
@@ -50,7 +50,8 @@ module subroutine damage_init()
   allocate(current(configHomogenizations%length))
 
   do ho = 1, configHomogenizations%length
-    allocate(current(ho)%phi(count(material_homogenizationID==ho)), source=1.0_pReal)
+    Nmembers = count(material_homogenizationID == ho)
+    allocate(current(ho)%phi(Nmembers), source=1.0_pReal)
     configHomogenization => configHomogenizations%get(ho)
     associate(prm => param(ho))
       if (configHomogenization%contains('damage')) then
@@ -60,10 +61,9 @@ module subroutine damage_init()
 #else
         prm%output = configHomogenizationDamage%get_as1dString('output',defaultVal=emptyStringArray)
 #endif
-        Nmaterialpoints = count(material_homogenizationAt == ho)
         damageState_h(ho)%sizeState = 1
-        allocate(damageState_h(ho)%state0(1,Nmaterialpoints), source=1.0_pReal)
+        allocate(damageState_h(ho)%state0(1,Nmembers), source=1.0_pReal)
-        allocate(damageState_h(ho)%state (1,Nmaterialpoints), source=1.0_pReal)
+        allocate(damageState_h(ho)%state (1,Nmembers), source=1.0_pReal)
       else
         prm%output = emptyStringArray
       endif

src/homogenization_mechanical_RGC.f90

@@ -78,7 +78,7 @@ module subroutine RGC_init(num_homogMech)
 
   integer :: &
     ho, &
-    Nmaterialpoints, &
+    Nmembers, &
     sizeState, nIntFaceTot
 
   class (tNode), pointer :: &
@@ -161,28 +161,28 @@ module subroutine RGC_init(num_homogMech)
     prm%D_alpha  = homogMech%get_as1dFloat('D_alpha', requiredSize=3)
     prm%a_g      = homogMech%get_as1dFloat('a_g',     requiredSize=3)
 
-    Nmaterialpoints = count(material_homogenizationAt == ho)
+    Nmembers = count(material_homogenizationID == ho)
     nIntFaceTot = 3*(  (prm%N_constituents(1)-1)*prm%N_constituents(2)*prm%N_constituents(3) &
                       + prm%N_constituents(1)*(prm%N_constituents(2)-1)*prm%N_constituents(3) &
                       + prm%N_constituents(1)*prm%N_constituents(2)*(prm%N_constituents(3)-1))
     sizeState = nIntFaceTot
 
     homogState(ho)%sizeState = sizeState
-    allocate(homogState(ho)%state0   (sizeState,Nmaterialpoints), source=0.0_pReal)
+    allocate(homogState(ho)%state0   (sizeState,Nmembers), source=0.0_pReal)
-    allocate(homogState(ho)%state    (sizeState,Nmaterialpoints), source=0.0_pReal)
+    allocate(homogState(ho)%state    (sizeState,Nmembers), source=0.0_pReal)
 
     stt%relaxationVector   => homogState(ho)%state(1:nIntFaceTot,:)
     st0%relaxationVector   => homogState(ho)%state0(1:nIntFaceTot,:)
 
-    allocate(dst%volumeDiscrepancy(   Nmaterialpoints), source=0.0_pReal)
+    allocate(dst%volumeDiscrepancy(   Nmembers), source=0.0_pReal)
-    allocate(dst%relaxationRate_avg(  Nmaterialpoints), source=0.0_pReal)
+    allocate(dst%relaxationRate_avg(  Nmembers), source=0.0_pReal)
-    allocate(dst%relaxationRate_max(  Nmaterialpoints), source=0.0_pReal)
+    allocate(dst%relaxationRate_max(  Nmembers), source=0.0_pReal)
-    allocate(dst%mismatch(          3,Nmaterialpoints), source=0.0_pReal)
+    allocate(dst%mismatch(          3,Nmembers), source=0.0_pReal)
 
 !--------------------------------------------------------------------------------------------------
 ! assigning cluster orientations
-    dependentState(ho)%orientation = spread(eu2om(prm%a_g*inRad),3,Nmaterialpoints)
+    dependentState(ho)%orientation = spread(eu2om(prm%a_g*inRad),3,Nmembers)
-    !dst%orientation = spread(eu2om(prm%a_g*inRad),3,Nmaterialpoints) ifort version 18.0.1 crashes (for whatever reason)
+    !dst%orientation = spread(eu2om(prm%a_g*inRad),3,Nmembers) ifort version 18.0.1 crashes (for whatever reason)
 
     end associate
 

src/homogenization_mechanical_isostrain.f90

@@ -15,7 +15,7 @@ module subroutine isostrain_init
 
   integer :: &
     ho, &
-    Nmaterialpoints
+    Nmembers
 
   print'(/,a)', ' <<<+-  homogenization:mechanical:isostrain init  -+>>>'
 
@@ -25,10 +25,10 @@ module subroutine isostrain_init
   do ho = 1, size(homogenization_type)
     if (homogenization_type(ho) /= HOMOGENIZATION_ISOSTRAIN_ID) cycle
 
-    Nmaterialpoints = count(material_homogenizationAt == ho)
+    Nmembers = count(material_homogenizationID == ho)
     homogState(ho)%sizeState = 0
-    allocate(homogState(ho)%state0(0,Nmaterialpoints))
+    allocate(homogState(ho)%state0(0,Nmembers))
-    allocate(homogState(ho)%state (0,Nmaterialpoints))
+    allocate(homogState(ho)%state (0,Nmembers))
 
   enddo
 

src/homogenization_mechanical_pass.f90

@@ -15,7 +15,7 @@ module subroutine pass_init
 
   integer :: &
     ho, &
-    Nmaterialpoints
+    Nmembers
 
   print'(/,a)', ' <<<+-  homogenization:mechanical:pass init  -+>>>'
 
@@ -28,10 +28,10 @@ module subroutine pass_init
     if(homogenization_Nconstituents(ho) /= 1) &
       call IO_error(211,ext_msg='N_constituents (pass)')
 
-    Nmaterialpoints = count(material_homogenizationAt == ho)
+    Nmembers = count(material_homogenizationID == ho)
     homogState(ho)%sizeState = 0
-    allocate(homogState(ho)%state0(0,Nmaterialpoints))
+    allocate(homogState(ho)%state0(0,Nmembers))
-    allocate(homogState(ho)%state (0,Nmaterialpoints))
+    allocate(homogState(ho)%state (0,Nmembers))
 
   enddo
 

src/lattice.f90

@@ -781,7 +781,7 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,structure) resul
        2, 4, 7, 5, 5, 7, 4, 2, 3, 1, 6, 6,  13,14,17, 9, 9,14,17,13,18, 8,11,18,  28,25,28,28,25,28,28,28,28,25,28,28,25,28,28,28,28,28,28,25,28,28,28,25, &
        7, 5, 2, 4, 7, 5, 2, 4, 6, 6, 1, 3,  17, 9,13,14,17,13, 9,14,11,18,18, 8,  28,28,28,25,28,28,25,28,28,28,28,25,28,28,25,28,28,25,28,28,28,25,28,28, &
        4, 2, 5, 7, 4, 2, 5, 7, 6, 6, 3, 1,  14,13, 9,17,14, 9,13,17, 8,18,18,11,  25,28,28,28,28,25,28,28,25,28,28,28,28,25,28,28,28,28,25,28,28,28,25,28, &
-      !                                                                                                                                                                                             
+
       19,19,10, 8, 9,12,16,15, 9,12,16,15,   1,20,24,24,23,22,21, 2,23,22, 2,21,  28,28,26,28,28,28,28,26,28,28,26,28,28,28,28,26,26,28,28,28,26,28,28,28, &
       19,19, 8,10,16,15, 9,12,16,15, 9,12,  20, 1,24,24,22,23, 2,21,22,23,21, 2,  28,28,28,26,28,28,26,28,28,28,28,26,28,28,26,28,28,26,28,28,28,26,28,28, &
       10, 8,19,19,12, 9,15,16,15,16,12, 9,  24,24, 1,20,21, 2,23,22, 2,21,23,22,  26,28,28,28,28,26,28,28,26,28,28,28,28,26,28,28,28,28,26,28,28,28,26,28, &
@@ -794,7 +794,7 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,structure) resul
        9,12,15,16,16,15,12, 9,10, 8,19,19,  21,23,22, 2, 2,23,22,21,24, 1,20,24,  28,26,28,28,26,28,28,28,28,26,28,28,26,28,28,28,28,28,28,26,28,28,28,26, &
       16,15,12, 9, 9,12,15,16, 8,10,19,19,   2,22,23,21,21,22,23, 2,24,20, 1,24,  28,28,26,28,28,28,28,26,28,28,26,28,28,28,28,26,26,28,28,28,26,28,28,28, &
       15,16, 9,12,15,16, 9,12,19,19, 8,10,  22, 2,21,23,22,21, 2,23,20,24,24, 1,  28,28,28,26,28,28,26,28,28,28,28,26,28,28,26,28,28,26,28,28,28,26,28,28, &
-      !                                                                                                                                                                                             
+                                                                                                                                  
       28,25,28,28,28,25,28,28,28,28,28,25,  28,28,26,28,28,26,28,28,26,28,28,28,   1,28,28,28,28,27,28,28,27,28,28,28,28,27,28,28,28,28,27,28,28,28,27,28, &
       25,28,28,28,28,28,28,25,28,25,28,28,  28,28,28,26,26,28,28,28,28,26,28,28,  28, 1,28,28,27,28,28,28,28,27,28,28,27,28,28,28,28,28,28,27,28,28,28,27, &
       28,28,28,25,25,28,28,28,25,28,28,28,  26,28,28,28,28,28,28,26,28,28,26,28,  28,28, 1,28,28,28,28,27,28,28,27,28,28,28,28,27,27,28,28,28,27,28,28,28, &
@@ -824,12 +824,12 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,structure) resul
                                                                                                     !< 2: collinear interaction --> alpha 1
                                                                                                     !< 3: coplanar interaction  --> alpha 2
                                                                                                     !< 4-7: other coefficients
-                                                                                                    !< 8: {110}-{112}, collinear and perpendicular planes --> alpha 6
+                                                                                                    !< 8: {110}-{112} collinear and perpendicular planes --> alpha 6
-                                                                                                    !< 9: {110}-{112}, just collinear                     --> alpha 7
+                                                                                                    !< 9: {110}-{112} collinear                          --> alpha 7
                                                                                                     !< 10-24: other coefficients
-                                                                                                    !< 25: {110}-{123}, just collinear
+                                                                                                    !< 25: {110}-{123} collinear
-                                                                                                    !< 26: {112}-{123}, just collinear
+                                                                                                    !< 26: {112}-{123} collinear
-                                                                                                    !< 27: {123}-{123}, just collinear
+                                                                                                    !< 27: {123}-{123} collinear
                                                                                                     !< 28: other interaction
 
   integer, dimension(HEX_NSLIP,HEX_NSLIP), parameter :: &
@@ -837,22 +837,22 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,structure) resul
        1, 2, 2,   3, 3, 3,   7, 7, 7,  13,13,13,13,13,13,  21,21,21,21,21,21,21,21,21,21,21,21,  31,31,31,31,31,31, & ! -----> acting (forest)
        2, 1, 2,   3, 3, 3,   7, 7, 7,  13,13,13,13,13,13,  21,21,21,21,21,21,21,21,21,21,21,21,  31,31,31,31,31,31, & ! |
        2, 2, 1,   3, 3, 3,   7, 7, 7,  13,13,13,13,13,13,  21,21,21,21,21,21,21,21,21,21,21,21,  31,31,31,31,31,31, & ! |
-     !                                                                                                                ! v
+                                                                                                                      ! v
        6, 6, 6,   4, 5, 5,   8, 8, 8,  14,14,14,14,14,14,  22,22,22,22,22,22,22,22,22,22,22,22,  32,32,32,32,32,32, & ! reacting (primary)
        6, 6, 6,   5, 4, 5,   8, 8, 8,  14,14,14,14,14,14,  22,22,22,22,22,22,22,22,22,22,22,22,  32,32,32,32,32,32, &
        6, 6, 6,   5, 5, 4,   8, 8, 8,  14,14,14,14,14,14,  22,22,22,22,22,22,22,22,22,22,22,22,  32,32,32,32,32,32, &
-     !
+
       12,12,12,  11,11,11,   9,10,10,  15,15,15,15,15,15,  23,23,23,23,23,23,23,23,23,23,23,23,  33,33,33,33,33,33, &
       12,12,12,  11,11,11,  10, 9,10,  15,15,15,15,15,15,  23,23,23,23,23,23,23,23,23,23,23,23,  33,33,33,33,33,33, &
       12,12,12,  11,11,11,  10,10, 9,  15,15,15,15,15,15,  23,23,23,23,23,23,23,23,23,23,23,23,  33,33,33,33,33,33, &
-     !
+
       20,20,20,  19,19,19,  18,18,18,  16,17,17,17,17,17,  24,24,24,24,24,24,24,24,24,24,24,24,  34,34,34,34,34,34, &
       20,20,20,  19,19,19,  18,18,18,  17,16,17,17,17,17,  24,24,24,24,24,24,24,24,24,24,24,24,  34,34,34,34,34,34, &
       20,20,20,  19,19,19,  18,18,18,  17,17,16,17,17,17,  24,24,24,24,24,24,24,24,24,24,24,24,  34,34,34,34,34,34, &
       20,20,20,  19,19,19,  18,18,18,  17,17,17,16,17,17,  24,24,24,24,24,24,24,24,24,24,24,24,  34,34,34,34,34,34, &
       20,20,20,  19,19,19,  18,18,18,  17,17,17,17,16,17,  24,24,24,24,24,24,24,24,24,24,24,24,  34,34,34,34,34,34, &
       20,20,20,  19,19,19,  18,18,18,  17,17,17,17,17,16,  24,24,24,24,24,24,24,24,24,24,24,24,  34,34,34,34,34,34, &
-     !
+
       30,30,30,  29,29,29,  28,28,28,  27,27,27,27,27,27,  25,26,26,26,26,26,26,26,26,26,26,26,  35,35,35,35,35,35, &
       30,30,30,  29,29,29,  28,28,28,  27,27,27,27,27,27,  26,25,26,26,26,26,26,26,26,26,26,26,  35,35,35,35,35,35, &
       30,30,30,  29,29,29,  28,28,28,  27,27,27,27,27,27,  26,26,25,26,26,26,26,26,26,26,26,26,  35,35,35,35,35,35, &
@@ -865,7 +865,7 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,structure) resul
       30,30,30,  29,29,29,  28,28,28,  27,27,27,27,27,27,  26,26,26,26,26,26,26,26,26,25,26,26,  35,35,35,35,35,35, &
       30,30,30,  29,29,29,  28,28,28,  27,27,27,27,27,27,  26,26,26,26,26,26,26,26,26,26,25,26,  35,35,35,35,35,35, &
       30,30,30,  29,29,29,  28,28,28,  27,27,27,27,27,27,  26,26,26,26,26,26,26,26,26,26,26,25,  35,35,35,35,35,35, &
-     !
+
       42,42,42,  41,41,41,  40,40,40,  39,39,39,39,39,39,  38,38,38,38,38,38,38,38,38,38,38,38,  36,37,37,37,37,37, &
       42,42,42,  41,41,41,  40,40,40,  39,39,39,39,39,39,  38,38,38,38,38,38,38,38,38,38,38,38,  37,36,37,37,37,37, &
       42,42,42,  41,41,41,  40,40,40,  39,39,39,39,39,39,  38,38,38,38,38,38,38,38,38,38,38,38,  37,37,36,37,37,37, &
@@ -881,34 +881,34 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,structure) resul
                                                                                                                                                                                                                                         ! |
         6,  6,   4,  5,   8,  8,  14, 14, 14, 14,  22, 22,  32, 32, 32, 32,  44, 44,  58, 58,  74, 74, 74, 74,  92, 92, 92, 92, 92, 92, 92, 92, 112, 112, 112, 112, 134,134,134,134,134,134,134,134, 158,158,158,158,158,158,158,158, & ! v
         6,  6,   5,  4,   8,  8,  14, 14, 14, 14,  22, 22,  32, 32, 32, 32,  44, 44,  58, 58,  74, 74, 74, 74,  92, 92, 92, 92, 92, 92, 92, 92, 112, 112, 112, 112, 134,134,134,134,134,134,134,134, 158,158,158,158,158,158,158,158, & ! reacting
-     !
+
        12, 12,  11, 11,   9, 10,  15, 15, 15, 15,  23, 23,  33, 33, 33, 33,  45, 45,  59, 59,  75, 75, 75, 75,  93, 93, 93, 93, 93, 93, 93, 93, 113, 113, 113, 113, 135,135,135,135,135,135,135,135, 159,159,159,159,159,159,159,159, &
        12, 12,  11, 11,  10,  9,  15, 15, 15, 15,  23, 23,  33, 33, 33, 33,  45, 45,  59, 59,  75, 75, 75, 75,  93, 93, 93, 93, 93, 93, 93, 93, 113, 113, 113, 113, 135,135,135,135,135,135,135,135, 159,159,159,159,159,159,159,159, &
-     !
+
        20, 20,  19, 19,  18, 18,  16, 17, 17, 17,  24, 24,  34, 34, 34, 34,  46, 46,  60, 60,  76, 76, 76, 76,  94, 94, 94, 94, 94, 94, 94, 94, 114, 114, 114, 114, 136,136,136,136,136,136,136,136, 160,160,160,160,160,160,160,160, &
        20, 20,  19, 19,  18, 18,  17, 16, 17, 17,  24, 24,  34, 34, 34, 34,  46, 46,  60, 60,  76, 76, 76, 76,  94, 94, 94, 94, 94, 94, 94, 94, 114, 114, 114, 114, 136,136,136,136,136,136,136,136, 160,160,160,160,160,160,160,160, &
        20, 20,  19, 19,  18, 18,  17, 17, 16, 17,  24, 24,  34, 34, 34, 34,  46, 46,  60, 60,  76, 76, 76, 76,  94, 94, 94, 94, 94, 94, 94, 94, 114, 114, 114, 114, 136,136,136,136,136,136,136,136, 160,160,160,160,160,160,160,160, &
        20, 20,  19, 19,  18, 18,  17, 17, 17, 16,  24, 24,  34, 34, 34, 34,  46, 46,  60, 60,  76, 76, 76, 76,  94, 94, 94, 94, 94, 94, 94, 94, 114, 114, 114, 114, 136,136,136,136,136,136,136,136, 160,160,160,160,160,160,160,160, &
-     !
+
        30, 30,  29, 29,  28, 28,  27, 27, 27, 27,  25, 26,  35, 35, 35, 35,  47, 47,  61, 61,  77, 77, 77, 77,  95, 95, 95, 95, 95, 95, 95, 95, 115, 115, 115, 115, 137,137,137,137,137,137,137,137, 161,161,161,161,161,161,161,161, &
        30, 30,  29, 29,  28, 28,  27, 27, 27, 27,  26, 25,  35, 35, 35, 35,  47, 47,  61, 61,  77, 77, 77, 77,  95, 95, 95, 95, 95, 95, 95, 95, 115, 115, 115, 115, 137,137,137,137,137,137,137,137, 161,161,161,161,161,161,161,161, &
-     !
+
        42, 42,  41, 41,  40, 40,  39, 39, 39, 39,  38, 38,  36, 37, 37, 37,  48, 48,  62, 62,  78, 78, 78, 78,  96, 96, 96, 96, 96, 96, 96, 96, 116, 116, 116, 116, 138,138,138,138,138,138,138,138, 162,162,162,162,162,162,162,162, &
        42, 42,  41, 41,  40, 40,  39, 39, 39, 39,  38, 38,  37, 36, 37, 37,  48, 48,  62, 62,  78, 78, 78, 78,  96, 96, 96, 96, 96, 96, 96, 96, 116, 116, 116, 116, 138,138,138,138,138,138,138,138, 162,162,162,162,162,162,162,162, &
        42, 42,  41, 41,  40, 40,  39, 39, 39, 39,  38, 38,  37, 37, 36, 37,  48, 48,  62, 62,  78, 78, 78, 78,  96, 96, 96, 96, 96, 96, 96, 96, 116, 116, 116, 116, 138,138,138,138,138,138,138,138, 162,162,162,162,162,162,162,162, &
        42, 42,  41, 41,  40, 40,  39, 39, 39, 39,  38, 38,  37, 37, 37, 36,  48, 48,  62, 62,  78, 78, 78, 78,  96, 96, 96, 96, 96, 96, 96, 96, 116, 116, 116, 116, 138,138,138,138,138,138,138,138, 162,162,162,162,162,162,162,162, &
-     !
+
        56, 56,  55, 55,  54, 54,  53, 53, 53, 53,  52, 52,  51, 51, 51, 51,  49, 50,  63, 63,  79, 79, 79, 79,  97, 97, 97, 97, 97, 97, 97, 97, 117, 117, 117, 117, 139,139,139,139,139,139,139,139, 163,163,163,163,163,163,163,163, &
        56, 56,  55, 55,  54, 54,  53, 53, 53, 53,  52, 52,  51, 51, 51, 51,  50, 49,  63, 63,  79, 79, 79, 79,  97, 97, 97, 97, 97, 97, 97, 97, 117, 117, 117, 117, 139,139,139,139,139,139,139,139, 163,163,163,163,163,163,163,163, &
-     !
+
        72, 72,  71, 71,  70, 70,  69, 69, 69, 69,  68, 68,  67, 67, 67, 67,  66, 66,  64, 65,  80, 80, 80, 80,  98, 98, 98, 98, 98, 98, 98, 98, 118, 118, 118, 118, 140,140,140,140,140,140,140,140, 164,164,164,164,164,164,164,164, &
        72, 72,  71, 71,  70, 70,  69, 69, 69, 69,  68, 68,  67, 67, 67, 67,  66, 66,  65, 64,  80, 80, 80, 80,  98, 98, 98, 98, 98, 98, 98, 98, 118, 118, 118, 118, 140,140,140,140,140,140,140,140, 164,164,164,164,164,164,164,164, &
-     !
+
        90, 90,  89, 89,  88, 88,  87, 87, 87, 87,  86, 86,  85, 85, 85, 85,  84, 84,  83, 83,  81, 82, 82, 82,  99, 99, 99, 99, 99, 99, 99, 99, 119, 119, 119, 119, 141,141,141,141,141,141,141,141, 165,165,165,165,165,165,165,165, &
        90, 90,  89, 89,  88, 88,  87, 87, 87, 87,  86, 86,  85, 85, 85, 85,  84, 84,  83, 83,  82, 81, 82, 82,  99, 99, 99, 99, 99, 99, 99, 99, 119, 119, 119, 119, 141,141,141,141,141,141,141,141, 165,165,165,165,165,165,165,165, &
        90, 90,  89, 89,  88, 88,  87, 87, 87, 87,  86, 86,  85, 85, 85, 85,  84, 84,  83, 83,  82, 82, 81, 82,  99, 99, 99, 99, 99, 99, 99, 99, 119, 119, 119, 119, 141,141,141,141,141,141,141,141, 165,165,165,165,165,165,165,165, &
        90, 90,  89, 89,  88, 88,  87, 87, 87, 87,  86, 86,  85, 85, 85, 85,  84, 84,  83, 83,  82, 82, 82, 81,  99, 99, 99, 99, 99, 99, 99, 99, 119, 119, 119, 119, 141,141,141,141,141,141,141,141, 165,165,165,165,165,165,165,165, &
-     !
+
       110,110, 109,109, 108,108, 107,107,107,107, 106,106, 105,105,105,105, 104,104, 103,103, 102,102,102,102, 100,101,101,101,101,101,101,101, 120, 120, 120, 120, 142,142,142,142,142,142,142,142, 166,166,166,166,166,166,166,166, &
       110,110, 109,109, 108,108, 107,107,107,107, 106,106, 105,105,105,105, 104,104, 103,103, 102,102,102,102, 101,100,101,101,101,101,101,101, 120, 120, 120, 120, 142,142,142,142,142,142,142,142, 166,166,166,166,166,166,166,166, &
       110,110, 109,109, 108,108, 107,107,107,107, 106,106, 105,105,105,105, 104,104, 103,103, 102,102,102,102, 101,101,100,101,101,101,101,101, 120, 120, 120, 120, 142,142,142,142,142,142,142,142, 166,166,166,166,166,166,166,166, &
@@ -917,12 +917,12 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,structure) resul
       110,110, 109,109, 108,108, 107,107,107,107, 106,106, 105,105,105,105, 104,104, 103,103, 102,102,102,102, 101,101,101,101,101,100,101,101, 120, 120, 120, 120, 142,142,142,142,142,142,142,142, 166,166,166,166,166,166,166,166, &
       110,110, 109,109, 108,108, 107,107,107,107, 106,106, 105,105,105,105, 104,104, 103,103, 102,102,102,102, 101,101,101,101,101,101,100,101, 120, 120, 120, 120, 142,142,142,142,142,142,142,142, 166,166,166,166,166,166,166,166, &
       110,110, 109,109, 108,108, 107,107,107,107, 106,106, 105,105,105,105, 104,104, 103,103, 102,102,102,102, 101,101,101,101,101,101,101,100, 120, 120, 120, 120, 142,142,142,142,142,142,142,142, 166,166,166,166,166,166,166,166, &
-     !
+
       132,132, 131,131, 130,130, 129,129,129,129, 128,128, 127,127,127,127, 126,126, 125,125, 124,124,124,124, 123,123,123,123,123,123,123,123, 121, 122, 122, 122, 143,143,143,143,143,143,143,143, 167,167,167,167,167,167,167,167, &
       132,132, 131,131, 130,130, 129,129,129,129, 128,128, 127,127,127,127, 126,126, 125,125, 124,124,124,124, 123,123,123,123,123,123,123,123, 121, 121, 122, 122, 143,143,143,143,143,143,143,143, 167,167,167,167,167,167,167,167, &
       132,132, 131,131, 130,130, 129,129,129,129, 128,128, 127,127,127,127, 126,126, 125,125, 124,124,124,124, 123,123,123,123,123,123,123,123, 121, 122, 121, 122, 143,143,143,143,143,143,143,143, 167,167,167,167,167,167,167,167, &
       132,132, 131,131, 130,130, 129,129,129,129, 128,128, 127,127,127,127, 126,126, 125,125, 124,124,124,124, 123,123,123,123,123,123,123,123, 121, 122, 122, 121, 143,143,143,143,143,143,143,143, 167,167,167,167,167,167,167,167, &
-     !
+
       156,156, 155,155, 154,154, 153,153,153,153, 152,152, 151,151,151,151, 150,150, 149,149, 148,148,148,148, 147,147,147,147,147,147,147,147, 146, 146, 146, 146, 144,145,145,145,145,145,145,145, 168,168,168,168,168,168,168,168, &
       156,156, 155,155, 154,154, 153,153,153,153, 152,152, 151,151,151,151, 150,150, 149,149, 148,148,148,148, 147,147,147,147,147,147,147,147, 146, 146, 146, 146, 145,144,145,145,145,145,145,145, 168,168,168,168,168,168,168,168, &
       156,156, 155,155, 154,154, 153,153,153,153, 152,152, 151,151,151,151, 150,150, 149,149, 148,148,148,148, 147,147,147,147,147,147,147,147, 146, 146, 146, 146, 145,145,144,145,145,145,145,145, 168,168,168,168,168,168,168,168, &
@@ -931,7 +931,7 @@ function lattice_interaction_SlipBySlip(Nslip,interactionValues,structure) resul
       156,156, 155,155, 154,154, 153,153,153,153, 152,152, 151,151,151,151, 150,150, 149,149, 148,148,148,148, 147,147,147,147,147,147,147,147, 146, 146, 146, 146, 145,145,145,145,145,144,145,145, 168,168,168,168,168,168,168,168, &
       156,156, 155,155, 154,154, 153,153,153,153, 152,152, 151,151,151,151, 150,150, 149,149, 148,148,148,148, 147,147,147,147,147,147,147,147, 146, 146, 146, 146, 145,145,145,145,145,145,144,145, 168,168,168,168,168,168,168,168, &
       156,156, 155,155, 154,154, 153,153,153,153, 152,152, 151,151,151,151, 150,150, 149,149, 148,148,148,148, 147,147,147,147,147,147,147,147, 146, 146, 146, 146, 145,145,145,145,145,145,145,144, 168,168,168,168,168,168,168,168, &
-     !
+
       182,182, 181,181, 180,180, 179,179,179,179, 178,178, 177,177,177,177, 176,176, 175,175, 174,174,174,174, 173,173,173,173,173,173,173,173, 172, 172, 172, 172, 171,171,171,171,171,171,171,171, 169,170,170,170,170,170,170,170, &
       182,182, 181,181, 180,180, 179,179,179,179, 178,178, 177,177,177,177, 176,176, 175,175, 174,174,174,174, 173,173,173,173,173,173,173,173, 172, 172, 172, 172, 171,171,171,171,171,171,171,171, 170,169,170,170,170,170,170,170, &
       182,182, 181,181, 180,180, 179,179,179,179, 178,178, 177,177,177,177, 176,176, 175,175, 174,174,174,174, 173,173,173,173,173,173,173,173, 172, 172, 172, 172, 171,171,171,171,171,171,171,171, 170,170,169,170,170,170,170,170, &
@@ -1021,21 +1021,21 @@ function lattice_interaction_TwinByTwin(Ntwin,interactionValues,structure) resul
        2, 2, 2, 1, 2, 2,   3, 3, 3, 3, 3, 3,   7, 7, 7, 7, 7, 7,  13,13,13,13,13,13, & ! v
        2, 2, 2, 2, 1, 2,   3, 3, 3, 3, 3, 3,   7, 7, 7, 7, 7, 7,  13,13,13,13,13,13, & ! reacting
        2, 2, 2, 2, 2, 1,   3, 3, 3, 3, 3, 3,   7, 7, 7, 7, 7, 7,  13,13,13,13,13,13, &
-     !
+
        6, 6, 6, 6, 6, 6,   4, 5, 5, 5, 5, 5,   8, 8, 8, 8, 8, 8,  14,14,14,14,14,14, &
        6, 6, 6, 6, 6, 6,   5, 4, 5, 5, 5, 5,   8, 8, 8, 8, 8, 8,  14,14,14,14,14,14, &
        6, 6, 6, 6, 6, 6,   5, 5, 4, 5, 5, 5,   8, 8, 8, 8, 8, 8,  14,14,14,14,14,14, &
        6, 6, 6, 6, 6, 6,   5, 5, 5, 4, 5, 5,   8, 8, 8, 8, 8, 8,  14,14,14,14,14,14, &
        6, 6, 6, 6, 6, 6,   5, 5, 5, 5, 4, 5,   8, 8, 8, 8, 8, 8,  14,14,14,14,14,14, &
        6, 6, 6, 6, 6, 6,   5, 5, 5, 5, 5, 4,   8, 8, 8, 8, 8, 8,  14,14,14,14,14,14, &
-     !
+
       12,12,12,12,12,12,  11,11,11,11,11,11,   9,10,10,10,10,10,  15,15,15,15,15,15, &
       12,12,12,12,12,12,  11,11,11,11,11,11,  10, 9,10,10,10,10,  15,15,15,15,15,15, &
       12,12,12,12,12,12,  11,11,11,11,11,11,  10,10, 9,10,10,10,  15,15,15,15,15,15, &
       12,12,12,12,12,12,  11,11,11,11,11,11,  10,10,10, 9,10,10,  15,15,15,15,15,15, &
       12,12,12,12,12,12,  11,11,11,11,11,11,  10,10,10,10, 9,10,  15,15,15,15,15,15, &
       12,12,12,12,12,12,  11,11,11,11,11,11,  10,10,10,10,10, 9,  15,15,15,15,15,15, &
-     !
+
       20,20,20,20,20,20,  19,19,19,19,19,19,  18,18,18,18,18,18,  16,17,17,17,17,17, &
       20,20,20,20,20,20,  19,19,19,19,19,19,  18,18,18,18,18,18,  17,16,17,17,17,17, &
       20,20,20,20,20,20,  19,19,19,19,19,19,  18,18,18,18,18,18,  17,17,16,17,17,17, &
@@ -1160,7 +1160,7 @@ function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,structure)
       3,3,3,2,2,3,3,3,3,2,3,3, &
       3,2,3,3,3,3,2,3,3,3,3,2, &
       3,3,2,3,3,2,3,3,2,3,3,3, &
-      !
+
       1,3,3,3,3,3,3,2,3,3,2,3, &
       3,1,3,3,3,3,2,3,3,3,3,2, &
       3,3,1,3,3,2,3,3,2,3,3,3, &
@@ -1173,7 +1173,7 @@ function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,structure)
       3,3,3,2,2,3,3,3,3,1,3,3, &
       2,3,3,3,3,3,3,2,3,3,1,3, &
       3,2,3,3,3,3,2,3,3,3,3,1, &
-      !
+
       4,4,4,4,4,4,4,4,4,4,4,4, &
       4,4,4,4,4,4,4,4,4,4,4,4, &
       4,4,4,4,4,4,4,4,4,4,4,4, &
@@ -1213,18 +1213,18 @@ function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,structure)
        5, 5, 5, 5, 5, 5,   6, 6, 6, 6, 6, 6,   7, 7, 7, 7, 7, 7,   8, 8, 8, 8, 8, 8, & ! slip (reacting)
        5, 5, 5, 5, 5, 5,   6, 6, 6, 6, 6, 6,   7, 7, 7, 7, 7, 7,   8, 8, 8, 8, 8, 8, &
        5, 5, 5, 5, 5, 5,   6, 6, 6, 6, 6, 6,   7, 7, 7, 7, 7, 7,   8, 8, 8, 8, 8, 8, &
-     !
+
        9, 9, 9, 9, 9, 9,  10,10,10,10,10,10,  11,11,11,11,11,11,  12,12,12,12,12,12, &
        9, 9, 9, 9, 9, 9,  10,10,10,10,10,10,  11,11,11,11,11,11,  12,12,12,12,12,12, &
        9, 9, 9, 9, 9, 9,  10,10,10,10,10,10,  11,11,11,11,11,11,  12,12,12,12,12,12, &
-     !
+
       13,13,13,13,13,13,  14,14,14,14,14,14,  15,15,15,15,15,15,  16,16,16,16,16,16, &
       13,13,13,13,13,13,  14,14,14,14,14,14,  15,15,15,15,15,15,  16,16,16,16,16,16, &
       13,13,13,13,13,13,  14,14,14,14,14,14,  15,15,15,15,15,15,  16,16,16,16,16,16, &
       13,13,13,13,13,13,  14,14,14,14,14,14,  15,15,15,15,15,15,  16,16,16,16,16,16, &
       13,13,13,13,13,13,  14,14,14,14,14,14,  15,15,15,15,15,15,  16,16,16,16,16,16, &
       13,13,13,13,13,13,  14,14,14,14,14,14,  15,15,15,15,15,15,  16,16,16,16,16,16, &
-     !
+
       17,17,17,17,17,17,  18,18,18,18,18,18,  19,19,19,19,19,19,  20,20,20,20,20,20, &
       17,17,17,17,17,17,  18,18,18,18,18,18,  19,19,19,19,19,19,  20,20,20,20,20,20, &
       17,17,17,17,17,17,  18,18,18,18,18,18,  19,19,19,19,19,19,  20,20,20,20,20,20, &
@@ -1237,14 +1237,13 @@ function lattice_interaction_SlipByTwin(Nslip,Ntwin,interactionValues,structure)
       17,17,17,17,17,17,  18,18,18,18,18,18,  19,19,19,19,19,19,  20,20,20,20,20,20, &
       17,17,17,17,17,17,  18,18,18,18,18,18,  19,19,19,19,19,19,  20,20,20,20,20,20, &
       17,17,17,17,17,17,  18,18,18,18,18,18,  19,19,19,19,19,19,  20,20,20,20,20,20, &
-     !
+
       21,21,21,21,21,21,  22,22,22,22,22,22,  23,23,23,23,23,23,  24,24,24,24,24,24, &
       21,21,21,21,21,21,  22,22,22,22,22,22,  23,23,23,23,23,23,  24,24,24,24,24,24, &
       21,21,21,21,21,21,  22,22,22,22,22,22,  23,23,23,23,23,23,  24,24,24,24,24,24, &
       21,21,21,21,21,21,  22,22,22,22,22,22,  23,23,23,23,23,23,  24,24,24,24,24,24, &
       21,21,21,21,21,21,  22,22,22,22,22,22,  23,23,23,23,23,23,  24,24,24,24,24,24, &
       21,21,21,21,21,21,  22,22,22,22,22,22,  23,23,23,23,23,23,  24,24,24,24,24,24  &
-     !
       ],shape(HEX_INTERACTIONSLIPTWIN))                                                             !< Slip-twin interaction types for hex
 
   select case(structure)
@@ -1352,21 +1351,21 @@ function lattice_interaction_TwinBySlip(Ntwin,Nslip,interactionValues,structure)
       1, 1, 1,   5, 5, 5,   9, 9, 9,  13,13,13,13,13,13,  17,17,17,17,17,17,17,17,17,17,17,17,  21,21,21,21,21,21, & ! v
       1, 1, 1,   5, 5, 5,   9, 9, 9,  13,13,13,13,13,13,  17,17,17,17,17,17,17,17,17,17,17,17,  21,21,21,21,21,21, & ! twin (reacting)
       1, 1, 1,   5, 5, 5,   9, 9, 9,  13,13,13,13,13,13,  17,17,17,17,17,17,17,17,17,17,17,17,  21,21,21,21,21,21, &
-     !
+
       2, 2, 2,   6, 6, 6,  10,10,10,  14,14,14,14,14,14,  18,18,18,18,18,18,18,18,18,18,18,18,  22,22,22,22,22,22, &
       2, 2, 2,   6, 6, 6,  10,10,10,  14,14,14,14,14,14,  18,18,18,18,18,18,18,18,18,18,18,18,  22,22,22,22,22,22, &
       2, 2, 2,   6, 6, 6,  10,10,10,  14,14,14,14,14,14,  18,18,18,18,18,18,18,18,18,18,18,18,  22,22,22,22,22,22, &
       2, 2, 2,   6, 6, 6,  10,10,10,  14,14,14,14,14,14,  18,18,18,18,18,18,18,18,18,18,18,18,  22,22,22,22,22,22, &
       2, 2, 2,   6, 6, 6,  10,10,10,  14,14,14,14,14,14,  18,18,18,18,18,18,18,18,18,18,18,18,  22,22,22,22,22,22, &
       2, 2, 2,   6, 6, 6,  10,10,10,  14,14,14,14,14,14,  18,18,18,18,18,18,18,18,18,18,18,18,  22,22,22,22,22,22, &
-     !
+
       3, 3, 3,   7, 7, 7,  11,11,11,  15,15,15,15,15,15,  19,19,19,19,19,19,19,19,19,19,19,19,  23,23,23,23,23,23, &
       3, 3, 3,   7, 7, 7,  11,11,11,  15,15,15,15,15,15,  19,19,19,19,19,19,19,19,19,19,19,19,  23,23,23,23,23,23, &
       3, 3, 3,   7, 7, 7,  11,11,11,  15,15,15,15,15,15,  19,19,19,19,19,19,19,19,19,19,19,19,  23,23,23,23,23,23, &
       3, 3, 3,   7, 7, 7,  11,11,11,  15,15,15,15,15,15,  19,19,19,19,19,19,19,19,19,19,19,19,  23,23,23,23,23,23, &
       3, 3, 3,   7, 7, 7,  11,11,11,  15,15,15,15,15,15,  19,19,19,19,19,19,19,19,19,19,19,19,  23,23,23,23,23,23, &
       3, 3, 3,   7, 7, 7,  11,11,11,  15,15,15,15,15,15,  19,19,19,19,19,19,19,19,19,19,19,19,  23,23,23,23,23,23, &
-     !
+
       4, 4, 4,   8, 8, 8,  12,12,12,  16,16,16,16,16,16,  20,20,20,20,20,20,20,20,20,20,20,20,  24,24,24,24,24,24, &
       4, 4, 4,   8, 8, 8,  12,12,12,  16,16,16,16,16,16,  20,20,20,20,20,20,20,20,20,20,20,20,  24,24,24,24,24,24, &
       4, 4, 4,   8, 8, 8,  12,12,12,  16,16,16,16,16,16,  20,20,20,20,20,20,20,20,20,20,20,20,  24,24,24,24,24,24, &

src/material.f90

@@ -16,30 +16,33 @@ module material
   implicit none
   private
 
+  type :: tRotationContainer
+    type(Rotation), dimension(:),  allocatable :: data
+  end type
+
+  type(tRotationContainer), dimension(:), allocatable :: material_orientation0
+
   integer, dimension(:), allocatable, public, protected :: &
     homogenization_Nconstituents                                                                    !< number of grains in each homogenization
+  integer, public, protected :: &
+    homogenization_maxNconstituents                                                                 !< max number of grains in any homogenization
 
   character(len=:), public, protected, allocatable, dimension(:) :: &
     material_name_phase, &                                                                          !< name of each phase
     material_name_homogenization                                                                    !< name of each homogenization
 
-  integer, public, protected :: &
-    homogenization_maxNconstituents                                                                 !< max number of grains in any USED homogenization
-
   integer, dimension(:),     allocatable, public, protected :: &                                    ! (elem)
-    material_homogenizationAt, &                                                                    !< homogenization ID of each element
+    material_homogenizationID, &                                                                    !< per cell TODO: material_ID_homogenization
-    material_homogenizationID, &                                                                    !< per cell
+    material_homogenizationEntry                                                                    !< per cell TODO: material_entry_homogenization
-    material_homogenizationEntry                                                                    !< per cell
-  integer, dimension(:,:),   allocatable :: &                                                       ! (ip,elem)
-    material_homogenizationMemberAt                                                                 !< position of the element within its homogenization instance
   integer, dimension(:,:),   allocatable, public, protected :: &                                    ! (constituent,elem)
-    material_phaseAt, &                                                                             !< phase ID of each element
+    material_phaseAt, &                                                                             !< phase ID of each element TODO: remove
-    material_phaseID, &                                                                             !< per (constituent,cell)
+    material_phaseID, &                                                                             !< per (constituent,cell) TODO: material_ID_phase
-    material_phaseEntry                                                                             !< per (constituent,cell)
+    material_phaseEntry                                                                             !< per (constituent,cell) TODO: material_entry_phase
   integer, dimension(:,:,:), allocatable, public, protected :: &                                    ! (constituent,IP,elem)
-    material_phaseMemberAt                                                                          !< position of the element within its phase instance
+    material_phaseMemberAt                                                                          !TODO: remove
-
   public :: &
+    tRotationContainer, &
+    material_orientation0, &
     material_init
 
 contains
@@ -55,25 +58,24 @@ subroutine material_init(restart)
   print'(/,a)', ' <<<+-  material init  -+>>>'; flush(IO_STDOUT)
 
 
-  call material_parseMaterial
+  call parse
-  print*, 'Material parsed'
+  print*, 'parsed material.yaml'
 
 
   if (.not. restart) then
     call results_openJobFile
-    call results_mapping_phase(material_phaseAt,material_phaseMemberAt,material_name_phase)
+    call results_mapping_phase(material_phaseID,material_phaseEntry,material_name_phase)
-    call results_mapping_homogenization(material_homogenizationAt,material_homogenizationMemberAt,material_name_homogenization)
+    call results_mapping_homogenization(material_homogenizationID,material_homogenizationEntry,material_name_homogenization)
     call results_closeJobFile
   endif
 
 end subroutine material_init
 
 
-
 !--------------------------------------------------------------------------------------------------
-!> @brief parses the material part in the material configuration file
+!> @brief Parse material.yaml to get the global structure
 !--------------------------------------------------------------------------------------------------
-subroutine material_parseMaterial
+subroutine parse()
 
   class(tNode), pointer :: materials, &                                                             !> list of materials
                            material, &                                                              !> material definition
@@ -90,7 +92,7 @@ subroutine material_parseMaterial
   real(pReal) :: &
     frac
   integer :: &
-    el, ip, co, &
+    el, ip, co, ma, &
     h, ce
 
   materials       => config_material%get('material')
@@ -111,8 +113,6 @@ subroutine material_parseMaterial
   allocate(counterPhase(phases%length),source=0)
   allocate(counterHomogenization(homogenizations%length),source=0)
 
-  allocate(material_homogenizationAt(discretization_Nelems),source=0)
-  allocate(material_homogenizationMemberAt(discretization_nIPs,discretization_Nelems),source=0)
   allocate(material_phaseAt(homogenization_maxNconstituents,discretization_Nelems),source=0)
   allocate(material_phaseMemberAt(homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems),source=0)
 
@@ -126,13 +126,11 @@ subroutine material_parseMaterial
     material     => materials%get(discretization_materialAt(el))
     constituents => material%get('constituents')
 
-    material_homogenizationAt(el) = homogenizations%getIndex(material%get_asString('homogenization'))
     do ip = 1, discretization_nIPs
       ce = (el-1)*discretization_nIPs + ip
-      counterHomogenization(material_homogenizationAt(el)) = counterHomogenization(material_homogenizationAt(el)) + 1
+      material_homogenizationID(ce) = homogenizations%getIndex(material%get_asString('homogenization'))
-      material_homogenizationMemberAt(ip,el)               = counterHomogenization(material_homogenizationAt(el))
+      counterHomogenization(material_homogenizationID(ce)) = counterHomogenization(material_homogenizationID(ce)) + 1
-      material_homogenizationID(ce)    = material_homogenizationAt(el)
+      material_homogenizationEntry(ce) = counterHomogenization(material_homogenizationID(ce))
-      material_homogenizationEntry(ce) = material_homogenizationMemberAt(ip,el)
     enddo
 
     frac = 0.0_pReal
@@ -145,9 +143,8 @@ subroutine material_parseMaterial
         ce = (el-1)*discretization_nIPs + ip
         counterPhase(material_phaseAt(co,el)) = counterPhase(material_phaseAt(co,el)) + 1
         material_phaseMemberAt(co,ip,el)      = counterPhase(material_phaseAt(co,el))
-
+        material_phaseEntry(co,ce) = counterPhase(material_phaseAt(co,el))
         material_phaseID(co,ce) = material_phaseAt(co,el)
-        material_phaseEntry(co,ce) = material_phaseMemberAt(co,ip,el)
       enddo
 
     enddo
@@ -155,7 +152,19 @@ subroutine material_parseMaterial
 
   enddo
 
-end subroutine material_parseMaterial
+  allocate(material_orientation0(materials%length))
+
+  do ma = 1, materials%length
+    material     => materials%get(ma)
+    constituents => material%get('constituents')
+    allocate(material_orientation0(ma)%data(constituents%length))
+    do co = 1, constituents%length
+      constituent => constituents%get(co)
+      call material_orientation0(ma)%data(co)%fromQuaternion(constituent%get_as1dFloat('O',requiredSize=4))
+    enddo
+ enddo
+
+end subroutine parse
 
 
 !--------------------------------------------------------------------------------------------------

src/math.f90

@@ -1140,8 +1140,7 @@ end function math_areaTriangle
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief limits a scalar value to a certain range (either one or two sided)
+!> @brief Limit a scalar value to a certain range (either one or two sided).
-! Will return NaN if left > right
 !--------------------------------------------------------------------------------------------------
 real(pReal) pure elemental function math_clip(a, left, right)
 

src/parallelization.f90

@@ -27,7 +27,7 @@ module parallelization
 contains
 
 !--------------------------------------------------------------------------------------------------
-!> @brief calls subroutines that reads material, numerics and debug configuration files
+!> @brief Initialize shared memory (openMP) and distributed memory (MPI) parallelization.
 !--------------------------------------------------------------------------------------------------
 subroutine parallelization_init
 

src/phase.f90

@@ -19,11 +19,9 @@ module phase
   implicit none
   private
 
-  type(Rotation), dimension(:,:,:), allocatable :: &
+  type(tRotationContainer), dimension(:), allocatable :: &
-    material_orientation0                                                                           !< initial orientation of each grain,IP,element
+    phase_orientation0, &
-
+    phase_orientation
-  type(rotation),            dimension(:,:,:),        allocatable :: &
-    crystallite_orientation                                                                         !< current orientation
 
   type :: tTensorContainer
     real(pReal), dimension(:,:,:), allocatable :: data
@@ -216,12 +214,11 @@ module phase
 
     end function thermal_stress
 
-    module function integrateDamageState(dt,co,ip,el) result(broken)
+    module function integrateDamageState(dt,co,ce) result(broken)
       real(pReal), intent(in) :: dt
       integer, intent(in) :: &
-        el, &                                                                                            !< element index in element loop
+        ce, &
-        ip, &                                                                                            !< integration point index in ip loop
+        co
-        co                                                                                               !< grain index in grain loop
       logical :: broken
     end function integrateDamageState
 
@@ -260,8 +257,7 @@ module phase
         ph, &
         i, &
         e
-      type(rotation), dimension(1,discretization_nIPs,discretization_Nelems), intent(in) :: &
+        type(tRotationContainer), dimension(:), intent(in) :: orientation
-        orientation                                                                                 !< crystal orientation
     end subroutine plastic_nonlocal_updateCompatibility
 
     module subroutine plastic_dependentState(co,ip,el)
@@ -348,7 +344,7 @@ contains
 subroutine phase_init
 
   integer :: &
-    ph
+    ph, ce, co, ma
   class (tNode), pointer :: &
     debug_constitutive, &
     materials, &
@@ -369,6 +365,25 @@ subroutine phase_init
   materials => config_material%get('material')
   phases    => config_material%get('phase')
 
+
+  allocate(phase_orientation0(phases%length))
+  do ph = 1,phases%length
+    allocate(phase_orientation0(ph)%data(count(material_phaseID==ph)))
+  enddo
+
+  do ce = 1, size(material_phaseID,2)
+    ma = discretization_materialAt((ce-1)/discretization_nIPs+1)
+    do co = 1,homogenization_Nconstituents(material_homogenizationID(ce))
+      ph = material_phaseID(co,ce)
+      phase_orientation0(ph)%data(material_phaseEntry(co,ce)) = material_orientation0(ma)%data(co)
+    enddo
+  enddo
+
+  allocate(phase_orientation(phases%length))
+  do ph = 1,phases%length
+    phase_orientation(ph)%data = phase_orientation0(ph)%data
+  enddo
+
   call mechanical_init(materials,phases)
   call damage_init
   call thermal_init(phases)
@@ -474,6 +489,7 @@ subroutine crystallite_init()
 
   integer :: &
     ph, &
+    ce, &
     co, &                                                                                           !< counter in integration point component loop
     ip, &                                                                                           !< counter in integration point loop
     el, &                                                                                           !< counter in element loop
@@ -492,8 +508,6 @@ subroutine crystallite_init()
   iMax = discretization_nIPs
   eMax = discretization_Nelems
 
-  allocate(crystallite_orientation(cMax,iMax,eMax))
-
   num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict)
 
   num%subStepMinCryst        = num_crystallite%get_asFloat ('subStepMin',       defaultVal=1.0e-3_pReal)
@@ -537,10 +551,11 @@ subroutine crystallite_init()
   flush(IO_STDOUT)
 
 
-  !$OMP PARALLEL DO
+  !$OMP PARALLEL DO PRIVATE(ce)
-  do el = 1, size(material_phaseMemberAt,3)
+  do el = 1, eMax
-    do ip = 1, size(material_phaseMemberAt,2)
+    do ip = 1, iMax
-      do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
+      ce = (el-1)*discretization_nIPs + ip
+      do co = 1,homogenization_Nconstituents(material_homogenizationID(ce))
         call crystallite_orientations(co,ip,el)
         call plastic_dependentState(co,ip,el)                                                       ! update dependent state variables to be consistent with basic states
      enddo
@@ -562,13 +577,16 @@ subroutine crystallite_orientations(co,ip,el)
     ip, &                                                                                           !< counter in integration point loop
     el                                                                                              !< counter in element loop
 
+  integer :: ph, en
 
-  call crystallite_orientation(co,ip,el)%fromMatrix(transpose(math_rotationalPart(&
+
-    mechanical_F_e(material_phaseAt(co,el),material_phaseMemberAt(co,ip,el)))))
+  ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
+  en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
+
+  call phase_orientation(ph)%data(en)%fromMatrix(transpose(math_rotationalPart(mechanical_F_e(ph,en))))
 
   if (plasticState(material_phaseAt(1,el))%nonlocal) &
-    call plastic_nonlocal_updateCompatibility(crystallite_orientation, &
+    call plastic_nonlocal_updateCompatibility(phase_orientation,material_phaseAt(1,el),ip,el)
-                                              material_phaseAt(1,el),ip,el)
 
 
 end subroutine crystallite_orientations
@@ -590,7 +608,7 @@ function crystallite_push33ToRef(co,ce, tensor33)
 
   ph = material_phaseID(co,ce)
   en = material_phaseEntry(co,ce)
-  T = matmul(material_orientation0(co,ph,en)%asMatrix(),transpose(math_inv33(phase_F(co,ce)))) ! ToDo: initial orientation correct?
+  T = matmul(phase_orientation0(ph)%data(en)%asMatrix(),transpose(math_inv33(phase_F(co,ce))))      ! ToDo: initial orientation correct?
 
   crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T))
 

src/phase_damage.f90

@@ -193,13 +193,12 @@ 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,ip,el) result(broken)
+module function integrateDamageState(dt,co,ce) result(broken)
 
   real(pReal), intent(in) :: dt
   integer, intent(in) :: &
-    el, &                                                                                            !< element index in element loop
+    ce, &
-    ip, &                                                                                            !< integration point index in ip loop
+    co
-    co                                                                                               !< grain index in grain loop
   logical :: broken
 
   integer :: &
@@ -215,8 +214,8 @@ module function integrateDamageState(dt,co,ip,el) result(broken)
   logical :: &
     converged_
 
-  ph = material_phaseAt(co,el)
+  ph = material_phaseID(co,ce)
-  me = material_phaseMemberAt(co,ip,el)
+  me = material_phaseEntry(co,ce)
 
   if (damageState(ph)%sizeState == 0) then
     broken = .false.

src/phase_damage_anisobrittle.f90

@@ -92,7 +92,7 @@ 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_phaseAt==p) * discretization_nIPs
+        Nmembers = count(material_phaseID==p)
         call phase_allocateState(damageState(p),Nmembers,1,1,0)
         damageState(p)%atol = src%get_asFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal)
         if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol'

src/phase_mechanical.f90

@@ -193,19 +193,13 @@ module subroutine mechanical_init(materials,phases)
     phases
 
   integer :: &
-    el, &
-    ip, &
-    co, &
-    ce, &
     ph, &
     en, &
-    stiffDegradationCtr, &
     Nmembers
   class(tNode), pointer :: &
     num_crystallite, &
     material, &
     constituents, &
-    constituent, &
     phase, &
     mech
 
@@ -229,8 +223,6 @@ module subroutine mechanical_init(materials,phases)
   allocate(phase_mechanical_P(phases%length))
   allocate(phase_mechanical_S0(phases%length))
 
-  allocate(material_orientation0(homogenization_maxNconstituents,phases%length,maxVal(material_phaseEntry)))
-
   do ph = 1, phases%length
     Nmembers = count(material_phaseID == ph)
 
@@ -258,18 +250,10 @@ module subroutine mechanical_init(materials,phases)
 #endif
   enddo
 
-  do el = 1, size(material_phaseMemberAt,3); do ip = 1, size(material_phaseMemberAt,2)
+  do ph = 1, phases%length
-    do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
+    do en = 1, count(material_phaseID == ph)
-      material     => materials%get(discretization_materialAt(el))
-      constituents => material%get('constituents')
-      constituent => constituents%get(co)
 
-      ph = material_phaseAt(co,el)
+      phase_mechanical_Fp0(ph)%data(1:3,1:3,en) = phase_orientation0(ph)%data(en)%asMatrix()        ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005)
-      en = material_phaseMemberAt(co,ip,el)
-
-      call material_orientation0(co,ph,en)%fromQuaternion(constituent%get_as1dFloat('O',requiredSize=4))
-
-      phase_mechanical_Fp0(ph)%data(1:3,1:3,en) = material_orientation0(co,ph,en)%asMatrix()                         ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005)
       phase_mechanical_Fp0(ph)%data(1:3,1:3,en) = phase_mechanical_Fp0(ph)%data(1:3,1:3,en) &
                                                 / math_det33(phase_mechanical_Fp0(ph)%data(1:3,1:3,en))**(1.0_pReal/3.0_pReal)
       phase_mechanical_Fi0(ph)%data(1:3,1:3,en) = math_I3
@@ -277,12 +261,11 @@ module subroutine mechanical_init(materials,phases)
 
       phase_mechanical_Fe(ph)%data(1:3,1:3,en) = math_inv33(matmul(phase_mechanical_Fi0(ph)%data(1:3,1:3,en), &
                                                                    phase_mechanical_Fp0(ph)%data(1:3,1:3,en)))  ! assuming that euler angles are given in internal strain free configuration
-      phase_mechanical_Fp(ph)%data(1:3,1:3,en) = phase_mechanical_Fp0(ph)%data(1:3,1:3,en)
-      phase_mechanical_Fi(ph)%data(1:3,1:3,en) = phase_mechanical_Fi0(ph)%data(1:3,1:3,en)
-      phase_mechanical_F(ph)%data(1:3,1:3,en)  = phase_mechanical_F0(ph)%data(1:3,1:3,en)
-
     enddo
-  enddo; enddo
+    phase_mechanical_Fp(ph)%data = phase_mechanical_Fp0(ph)%data
+    phase_mechanical_Fi(ph)%data = phase_mechanical_Fi0(ph)%data
+    phase_mechanical_F(ph)%data  = phase_mechanical_F0(ph)%data
+  enddo
 
 
 ! initialize elasticity
@@ -433,8 +416,8 @@ function integrateStress(F,subFp0,subFi0,Delta_t,co,ip,el) result(broken)
 
   broken = .true.
 
-  ph = material_phaseAt(co,el)
+  ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
-  en = material_phaseMemberAt(co,ip,el)
+  en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
 
   call plastic_dependentState(co,ip,el)
 
@@ -609,8 +592,8 @@ function integrateStateFPI(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) resul
     dotState
 
 
-  ph = material_phaseAt(co,el)
+  ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
-  en = material_phaseMemberAt(co,ip,el)
+  en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
 
   broken = plastic_dotState(Delta_t, co,ip,el,ph,en)
   if(broken) return
@@ -694,8 +677,8 @@ function integrateStateEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el) res
     sizeDotState
 
 
-  ph = material_phaseAt(co,el)
+  ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
-  en = material_phaseMemberAt(co,ip,el)
+  en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
 
   broken = plastic_dotState(Delta_t, co,ip,el,ph,en)
   if(broken) return
@@ -734,8 +717,8 @@ function integrateStateAdaptiveEuler(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip
   real(pReal), dimension(phase_plasticity_maxSizeDotState) :: residuum_plastic
 
 
-  ph = material_phaseAt(co,el)
+  ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
-  en = material_phaseMemberAt(co,ip,el)
+  en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
 
   broken = plastic_dotState(Delta_t, co,ip,el,ph,en)
   if(broken) return
@@ -852,8 +835,8 @@ function integrateStateRK(F_0,F,subFp0,subFi0,subState0,Delta_t,co,ip,el,A,B,C,D
   real(pReal), dimension(phase_plasticity_maxSizeDotState,size(B)) :: plastic_RKdotState
 
 
-  ph = material_phaseAt(co,el)
+  ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
-  en = material_phaseMemberAt(co,ip,el)
+  en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
 
   broken = plastic_dotState(Delta_t,co,ip,el,ph,en)
   if(broken) return
@@ -956,7 +939,7 @@ subroutine crystallite_results(group,ph)
             case(lattice_BCT_ID)
               structureLabel = 'tI'
           end select
-          selected_rotations = select_rotations(crystallite_orientation,ph)
+          selected_rotations = select_rotations(phase_orientation(ph)%data)
           call results_writeDataset(group//'/mechanical',selected_rotations,output_constituent(ph)%label(ou),&
                                    'crystal orientation as quaternion','q_0 (q_1 q_2 q_3)')
           call results_addAttribute('lattice',structureLabel,group//'/mechanical/'//output_constituent(ph)%label(ou))
@@ -967,27 +950,16 @@ subroutine crystallite_results(group,ph)
   contains
 
 !--------------------------------------------------------------------------------------------------
-!> @brief select rotations for output
+!> @brief Convert orientation for output: ToDo: implement in HDF5/results
 !--------------------------------------------------------------------------------------------------
-  function select_rotations(dataset,ph)
+  function select_rotations(dataset)
 
-    integer, intent(in) :: ph
+    type(rotation), dimension(:), intent(in) :: dataset
-    type(rotation), dimension(:,:,:), intent(in) :: dataset
+    real(pReal), dimension(4,size(dataset,1)) :: select_rotations
-    real(pReal), allocatable, dimension(:,:) :: select_rotations
+    integer :: en
-    integer :: el,ip,co,j
 
-    allocate(select_rotations(4,count(material_phaseAt==ph)*homogenization_maxNconstituents*discretization_nIPs))
+    do en = 1, size(dataset,1)
-
+      select_rotations(:,en) = dataset(en)%asQuaternion()
-    j=0
-    do el = 1, size(material_phaseAt,2)
-      do ip = 1, discretization_nIPs
-        do co = 1, size(material_phaseAt,1)                                                          !ToDo: this needs to be changed for varying Ngrains
-           if (material_phaseAt(co,el) == ph) then
-             j = j + 1
-             select_rotations(1:4,j) = dataset(co,ip,el)%asQuaternion()
-           endif
-        enddo
-      enddo
    enddo
 
  end function select_rotations
@@ -1045,8 +1017,8 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
   real(pReal), dimension(:), allocatable :: subState0
 
 
-  ph = material_phaseAt(co,el)
+  ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
-  en = material_phaseMemberAt(co,ip,el)
+  en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
   sizeDotState = plasticState(ph)%sizeDotState
 
   subLi0 = phase_mechanical_Li0(ph)%data(1:3,1:3,en)
@@ -1109,7 +1081,7 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
       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_ = converged_ .and. .not. integrateDamageState(subStep * dt,co,ip,el)
+      converged_ = converged_ .and. .not. integrateDamageState(subStep * dt,co,(el-1)*discretization_nIPs + ip)
     endif
 
   enddo cutbackLooping

src/phase_mechanical_plastic.f90

@@ -360,10 +360,10 @@ module subroutine plastic_dependentState(co, ip, el)
     en
 
 
-  ph = material_phaseAt(co,el)
+  ph = material_phaseID(co,(el-1)*discretization_nIPs + ip)
-  en = material_phasememberAt(co,ip,el)
+  en = material_phaseEntry(co,(el-1)*discretization_nIPs + ip)
 
-  plasticType: select case (phase_plasticity(material_phaseAt(co,el)))
+  plasticType: select case (phase_plasticity(ph))
 
     case (PLASTICITY_DISLOTWIN_ID) plasticType
       call dislotwin_dependentState(thermal_T(ph,en),ph,en)

src/phase_mechanical_plastic_nonlocal.f90

@@ -506,7 +506,6 @@ module function plastic_nonlocal_init() result(myPlasticity)
     end associate
 
     if (Nmembers > 0) call stateInit(ini,ph,Nmembers)
-    plasticState(ph)%state0 = plasticState(ph)%state
 
 !--------------------------------------------------------------------------------------------------
 !  exit if any parameter is out of range
@@ -623,7 +622,7 @@ module subroutine nonlocal_dependentState(ph, en, ip, el)
 
   ! coefficients are corrected for the line tension effect
   ! (see Kubin,Devincre,Hoc; 2008; Modeling dislocation storage rates and mean free paths in face-centered cubic crystals)
-  if (any(lattice_structure(material_phaseAt(1,el)) == [LATTICE_bcc_ID,LATTICE_fcc_ID])) then
+  if (any(lattice_structure(ph) == [LATTICE_bcc_ID,LATTICE_fcc_ID])) then
     myInteractionMatrix = prm%h_sl_sl &
                         * spread((  1.0_pReal - prm%f_F &
                                    + prm%f_F &
@@ -644,7 +643,7 @@ module subroutine nonlocal_dependentState(ph, en, ip, el)
   !#################################################################################################
 
   rho0 = getRho0(ph,en)
-  if (.not. phase_localPlasticity(material_phaseAt(1,el)) .and. prm%shortRangeStressCorrection) then
+  if (.not. phase_localPlasticity(ph) .and. prm%shortRangeStressCorrection) then
     invFp = math_inv33(phase_mechanical_Fp(ph)%data(1:3,1:3,en))
     invFe = math_inv33(phase_mechanical_Fe(ph)%data(1:3,1:3,en))
 
@@ -1238,7 +1237,7 @@ function rhoDotFlux(timestep,ph,en,ip,el)
   !****************************************************************************
   !*** calculate dislocation fluxes (only for nonlocal plasticity)
   rhoDotFlux = 0.0_pReal
-  if (.not. phase_localPlasticity(material_phaseAt(1,el))) then
+  if (.not. phase_localPlasticity(ph)) then
 
     !*** check CFL (Courant-Friedrichs-Lewy) condition for flux
     if (any( abs(gdot) > 0.0_pReal &                                                                ! any active slip system ...
@@ -1394,7 +1393,7 @@ end function rhoDotFlux
 !--------------------------------------------------------------------------------------------------
 module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
 
-  type(rotation), dimension(1,discretization_nIPs,discretization_Nelems), intent(in) :: &
+  type(tRotationContainer), dimension(:), intent(in) :: &
     orientation                                                                                     ! crystal orientation
   integer, intent(in) :: &
     ph, &
@@ -1450,8 +1449,8 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
     elseif (prm%chi_GB >= 0.0_pReal) then
       !* GRAIN BOUNDARY !
       !* fixed transmissivity for adjacent ips with different texture (only if explicitly given in material.config)
-      if (any(dNeq(material_orientation0(1,ph,en)%asQuaternion(), &
+      if (any(dNeq(phase_orientation0(ph)%data(en)%asQuaternion(), &
-                   material_orientation0(1,neighbor_phase,neighbor_me)%asQuaternion())) .and. &
+                   phase_orientation0(neighbor_phase)%data(neighbor_me)%asQuaternion())) .and. &
           (.not. phase_localPlasticity(neighbor_phase))) &
         forall(s1 = 1:ns) my_compatibility(:,s1,s1,n) = sqrt(prm%chi_GB)
     else
@@ -1464,7 +1463,7 @@ module subroutine plastic_nonlocal_updateCompatibility(orientation,ph,i,e)
       !* the number of compatible slip systems is minimized with the sum of the original compatibility values exceeding one.
       !* Finally the smallest compatibility value is decreased until the sum is exactly equal to one.
       !* All values below the threshold are set to zero.
-      mis = orientation(1,i,e)%misorientation(orientation(1,neighbor_i,neighbor_e))
+      mis = orientation(ph)%data(en)%misorientation(orientation(neighbor_phase)%data(neighbor_me))
       mySlipSystems: do s1 = 1,ns
         neighborSlipSystems: do s2 = 1,ns
           my_compatibility(1,s2,s1,n) =     math_inner(prm%slip_normal(1:3,s1), &
@@ -1579,21 +1578,20 @@ end subroutine plastic_nonlocal_results
 !--------------------------------------------------------------------------------------------------
 !> @brief populates the initial dislocation density
 !--------------------------------------------------------------------------------------------------
-subroutine stateInit(ini,phase,Nmembers)
+subroutine stateInit(ini,phase,Nentries)
 
   type(tInitialParameters) :: &
     ini
   integer,intent(in) :: &
     phase, &
-    Nmembers
+    Nentries
+
   integer :: &
-    i, &
     e, &
     f, &
     from, &
     upto, &
-    s, &
+    s
-    phasemember
   real(pReal), dimension(2) :: &
     noise, &
     rnd
@@ -1602,33 +1600,26 @@ subroutine stateInit(ini,phase,Nmembers)
     totalVolume, &
     densityBinning, &
     minimumIpVolume
-  real(pReal), dimension(Nmembers) :: &
-    volume
 
 
   associate(stt => state(phase))
 
     if (ini%random_rho_u > 0.0_pReal) then ! randomly distribute dislocation segments on random slip system and of random type in the volume
-    do e = 1,discretization_Nelems
+      totalVolume     = sum(geom(phase)%V_0)
-      do i = 1,discretization_nIPs
+      minimumIPVolume = minval(geom(phase)%V_0)
-        if (material_phaseAt(1,e) == phase) volume(material_phasememberAt(1,i,e)) = IPvolume(i,e)
-      enddo
-    enddo
-    totalVolume     = sum(volume)
-    minimumIPVolume = minval(volume)
       densityBinning  = ini%random_rho_u_binning / minimumIpVolume ** (2.0_pReal / 3.0_pReal)
 
       ! fill random material points with dislocation segments until the desired overall density is reached
       meanDensity = 0.0_pReal
       do while(meanDensity < ini%random_rho_u)
         call random_number(rnd)
-      phasemember = nint(rnd(1)*real(Nmembers,pReal) + 0.5_pReal)
+        e = nint(rnd(1)*real(Nentries,pReal) + 0.5_pReal)
         s = nint(rnd(2)*real(sum(ini%N_sl),pReal)*4.0_pReal + 0.5_pReal)
-      meanDensity = meanDensity + densityBinning * volume(phasemember) / totalVolume
+        meanDensity = meanDensity + densityBinning * geom(phase)%V_0(e) / totalVolume
-      stt%rhoSglMobile(s,phasemember) = densityBinning
+        stt%rhoSglMobile(s,e) = densityBinning
       enddo
     else                                ! homogeneous distribution with noise
-    do e = 1, Nmembers
+      do e = 1, Nentries
         do f = 1,size(ini%N_sl,1)
           from = 1 + sum(ini%N_sl(1:f-1))
           upto = sum(ini%N_sl(1:f))

src/prec.f90

@@ -78,7 +78,7 @@ contains
 
 
 !--------------------------------------------------------------------------------------------------
-!> @brief report precision and do self test
+!> @brief Report precision and do self test.
 !--------------------------------------------------------------------------------------------------
 subroutine prec_init
 

src/results.f90

@@ -415,16 +415,15 @@ end subroutine results_writeTensorDataset_int
 !--------------------------------------------------------------------------------------------------
 !> @brief adds the unique mapping from spatial position and constituent ID to results
 !--------------------------------------------------------------------------------------------------
-subroutine results_mapping_phase(phaseAt,memberAtLocal,label)
+subroutine results_mapping_phase(ID,entry,label)
 
-  integer,          dimension(:,:),   intent(in) :: phaseAt                                         !< phase section at (constituent,element)
+  integer,          dimension(:,:), intent(in) :: ID                                                !< phase ID at (co,ce)
-  integer,          dimension(:,:,:), intent(in) :: memberAtLocal                                   !< phase member at (constituent,IP,element)
+  integer,          dimension(:,:), intent(in) :: entry                                             !< phase entry at (co,ce)
   character(len=*), dimension(:),   intent(in) :: label                                             !< label of each phase section
 
-  integer, dimension(size(memberAtLocal,1),size(memberAtLocal,2),size(memberAtLocal,3)) :: &
+  integer, dimension(size(entry,1),size(entry,2)) :: &
-    phaseAtMaterialpoint, &
+    entryGlobal
-    memberAtGlobal
+  integer, dimension(size(label),0:worldsize-1) :: entryOffset                                      !< offset in entry counting per process
-  integer, dimension(size(label),0:worldsize-1) :: memberOffset                                     !< offset in member counting per process
   integer, dimension(0:worldsize-1)             :: writeSize                                        !< amount of data written per process
   integer(HSIZE_T), dimension(2) :: &
     myShape, &                                                                                      !< shape of the dataset (this process)
@@ -443,51 +442,48 @@ subroutine results_mapping_phase(phaseAt,memberAtLocal,label)
     dt_id
 
   integer(SIZE_T) :: type_size_string, type_size_int
-  integer         :: hdferr, ierr, i
+  integer         :: hdferr, ierr, ce, co
+
+
+  writeSize = 0
+  writeSize(worldrank) = size(entry(1,:))                                                           ! total number of entries of this process
 
-!--------------------------------------------------------------------------------------------------
-! prepare MPI communication (transparent for non-MPI runs)
   call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  memberOffset = 0
-  do i=1, size(label)
-    memberOffset(i,worldrank) = count(phaseAt == i)*size(memberAtLocal,2)                           ! number of points/instance of this process
-  enddo
-  writeSize = 0
-  writeSize(worldrank) = size(memberAtLocal(1,:,:))                                                 ! total number of points by this process
 
+#ifndef PETSc
+  entryGlobal = entry -1                                                                            ! 0-based
+#else
 !--------------------------------------------------------------------------------------------------
 ! MPI settings and communication
-#ifdef PETSc
   call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
   call MPI_allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)        ! get output at each process
   if(ierr /= 0) error stop 'MPI error'
 
-  call MPI_allreduce(MPI_IN_PLACE,memberOffset,size(memberOffset),MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)! get offset at each process
+  entryOffset = 0
+  do co = 1, size(ID,1)
+    do ce = 1, size(ID,2)
+      entryOffset(ID(co,ce),worldrank) = entryOffset(ID(co,ce),worldrank) +1
+    enddo
+  enddo
+  call MPI_allreduce(MPI_IN_PLACE,entryOffset,size(entryOffset),MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)! get offset at each process
   if(ierr /= 0) error stop 'MPI error'
+  entryOffset(:,worldrank) = sum(entryOffset(:,0:worldrank-1),2)
+  do co = 1, size(ID,1)
+    do ce = 1, size(ID,2)
+      entryGlobal(co,ce) = entry(co,ce) -1 + entryOffset(ID(co,ce),worldrank)
+    enddo
+  enddo
 #endif
 
-  myShape    = int([size(phaseAt,1),writeSize(worldrank)],  HSIZE_T)
+  myShape = int([size(ID,1),writeSize(worldrank)], HSIZE_T)
   myOffset = int([0,sum(writeSize(0:worldrank-1))], HSIZE_T)
-  totalShape = int([size(phaseAt,1),sum(writeSize)],        HSIZE_T)
+  totalShape = int([size(ID,1),sum(writeSize)], HSIZE_T)
-
 
 !---------------------------------------------------------------------------------------------------
-! expand phaseAt to consider IPs (is not stored per IP)
+! compound type: label(ID) + entry
-  do i = 1, size(phaseAtMaterialpoint,2)
-    phaseAtMaterialpoint(:,i,:) = phaseAt
-  enddo
-
-!---------------------------------------------------------------------------------------------------
-! renumber member from my process to all processes
-  do i = 1, size(label)
-    where(phaseAtMaterialpoint == i) memberAtGlobal = memberAtLocal + sum(memberOffset(i,0:worldrank-1)) -1     ! convert to 0-based
-  enddo
-
-!---------------------------------------------------------------------------------------------------
-! compound type: name of phase section + position/index within results array
   call h5tcopy_f(H5T_NATIVE_CHARACTER, dt_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
   call h5tset_size_f(dt_id, int(len(label(1)),SIZE_T), hdferr)
@@ -540,10 +536,10 @@ subroutine results_mapping_phase(phaseAt,memberAtLocal,label)
   call h5dcreate_f(loc_id, 'phase', dtype_id, filespace_id, dset_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
-  call h5dwrite_f(dset_id, label_id, reshape(label(pack(phaseAtMaterialpoint,.true.)),myShape), &
+  call h5dwrite_f(dset_id, label_id, reshape(label(pack(ID,.true.)),myShape), &
                   myShape, hdferr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5dwrite_f(dset_id, entry_id, reshape(pack(memberAtGlobal,.true.),myShape), &
+  call h5dwrite_f(dset_id, entry_id, reshape(pack(entryGlobal,.true.),myShape), &
                   myShape, hdferr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
   if(hdferr < 0) error stop 'HDF5 error'
 
@@ -572,16 +568,15 @@ end subroutine results_mapping_phase
 !--------------------------------------------------------------------------------------------------
 !> @brief adds the unique mapping from spatial position and constituent ID to results
 !--------------------------------------------------------------------------------------------------
-subroutine results_mapping_homogenization(homogenizationAt,memberAtLocal,label)
+subroutine results_mapping_homogenization(ID,entry,label)
 
-  integer,          dimension(:),   intent(in) :: homogenizationAt                                  !< homogenization section at (element)
+  integer,          dimension(:), intent(in) :: ID                                                  !< homogenization ID at (ce)
-  integer,          dimension(:,:), intent(in) :: memberAtLocal                                     !< homogenization member at (IP,element)
+  integer,          dimension(:), intent(in) :: entry                                               !< homogenization entry at (ce)
   character(len=*), dimension(:), intent(in) :: label                                               !< label of each homogenization section
 
-  integer, dimension(size(memberAtLocal,1),size(memberAtLocal,2)) :: &
+  integer, dimension(size(entry,1)) :: &
-    homogenizationAtMaterialpoint, &
+    entryGlobal
-    memberAtGlobal
+  integer, dimension(size(label),0:worldsize-1) :: entryOffset                                      !< offset in entry counting per process
-  integer, dimension(size(label),0:worldsize-1) :: memberOffset                                     !< offset in member counting per process
   integer, dimension(0:worldsize-1)             :: writeSize                                        !< amount of data written per process
   integer(HSIZE_T), dimension(1) :: &
     myShape, &                                                                                      !< shape of the dataset (this process)
@@ -600,52 +595,44 @@ subroutine results_mapping_homogenization(homogenizationAt,memberAtLocal,label)
     dt_id
 
   integer(SIZE_T) :: type_size_string, type_size_int
-  integer         :: hdferr, ierr, i
+  integer         :: hdferr, ierr, ce
 
 
-!--------------------------------------------------------------------------------------------------
+  writeSize = 0
-! prepare MPI communication (transparent for non-MPI runs)
+  writeSize(worldrank) = size(entry)                                                                ! total number of entries of this process
+
   call h5pcreate_f(H5P_DATASET_XFER_F, plist_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
-  memberOffset = 0
-  do i=1, size(label)
-    memberOffset(i,worldrank) = count(homogenizationAt == i)*size(memberAtLocal,1)                  ! number of points/instance of this process
-  enddo
-  writeSize = 0
-  writeSize(worldrank) = size(memberAtLocal)                                                        ! total number of points by this process
 
+#ifndef PETSc
+  entryGlobal = entry -1                                                                            ! 0-based
+#else
 !--------------------------------------------------------------------------------------------------
 ! MPI settings and communication
-#ifdef PETSc
   call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
   call MPI_allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)        ! get output at each process
   if(ierr /= 0) error stop 'MPI error'
 
-  call MPI_allreduce(MPI_IN_PLACE,memberOffset,size(memberOffset),MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)! get offset at each process
+  entryOffset = 0
+  do ce = 1, size(ID,1)
+    entryOffset(ID(ce),worldrank) = entryOffset(ID(ce),worldrank) +1
+  enddo
+  call MPI_allreduce(MPI_IN_PLACE,entryOffset,size(entryOffset),MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get offset at each process
   if(ierr /= 0) error stop 'MPI error'
+  entryOffset(:,worldrank) = sum(entryOffset(:,0:worldrank-1),2)
+  do ce = 1, size(ID,1)
+    entryGlobal(ce) = entry(ce) -1 + entryOffset(ID(ce),worldrank)
+  enddo
 #endif
 
   myShape = int([writeSize(worldrank)], HSIZE_T)
   myOffset = int([sum(writeSize(0:worldrank-1))], HSIZE_T)
   totalShape = int([sum(writeSize)], HSIZE_T)
 
-
 !---------------------------------------------------------------------------------------------------
-! expand phaseAt to consider IPs (is not stored per IP)
+! compound type: label(ID) + entry
-  do i = 1, size(homogenizationAtMaterialpoint,1)
-    homogenizationAtMaterialpoint(i,:) = homogenizationAt
-  enddo
-
-!---------------------------------------------------------------------------------------------------
-! renumber member from my process to all processes
-  do i = 1, size(label)
-    where(homogenizationAtMaterialpoint == i) memberAtGlobal = memberAtLocal + sum(memberOffset(i,0:worldrank-1)) - 1  ! convert to 0-based
-  enddo
-
-!---------------------------------------------------------------------------------------------------
-! compound type: name of phase section + position/index within results array
   call h5tcopy_f(H5T_NATIVE_CHARACTER, dt_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
   call h5tset_size_f(dt_id, int(len(label(1)),SIZE_T), hdferr)
@@ -698,10 +685,10 @@ subroutine results_mapping_homogenization(homogenizationAt,memberAtLocal,label)
   call h5dcreate_f(loc_id, 'homogenization', dtype_id, filespace_id, dset_id, hdferr)
   if(hdferr < 0) error stop 'HDF5 error'
 
-  call h5dwrite_f(dset_id, label_id, reshape(label(pack(homogenizationAtMaterialpoint,.true.)),myShape), &
+  call h5dwrite_f(dset_id, label_id, reshape(label(pack(ID,.true.)),myShape), &
                   myShape, hdferr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
   if(hdferr < 0) error stop 'HDF5 error'
-  call h5dwrite_f(dset_id, entry_id, reshape(pack(memberAtGlobal,.true.),myShape), &
+  call h5dwrite_f(dset_id, entry_id, reshape(pack(entryGlobal,.true.),myShape), &
                   myShape, hdferr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
   if(hdferr < 0) error stop 'HDF5 error'
 
@@ -733,19 +720,18 @@ end subroutine results_mapping_homogenization
 !--------------------------------------------------------------------------------------------------
 subroutine executionStamp(path,description,SIunit)
 
-
   character(len=*), intent(in)           :: path,description
   character(len=*), intent(in), optional :: SIunit
 
 
-  if (HDF5_objectExists(resultsFile,path)) &
-    call HDF5_addAttribute(resultsFile,'description',description,path)
-  if (HDF5_objectExists(resultsFile,path) .and. present(SIunit)) &
-    call HDF5_addAttribute(resultsFile,'unit',SIunit,path)
   if (HDF5_objectExists(resultsFile,path)) &
     call HDF5_addAttribute(resultsFile,'creator','DAMASK '//DAMASKVERSION,path)
   if (HDF5_objectExists(resultsFile,path)) &
     call HDF5_addAttribute(resultsFile,'created',now(),path)
+  if (HDF5_objectExists(resultsFile,path)) &
+    call HDF5_addAttribute(resultsFile,'description',description,path)
+  if (HDF5_objectExists(resultsFile,path) .and. present(SIunit)) &
+    call HDF5_addAttribute(resultsFile,'unit',SIunit,path)
 
 end subroutine executionStamp