reverting back to single level gradually

code/kinematics_cleavage_opening.f90

@@ -0,0 +1,303 @@
+!--------------------------------------------------------------------------------------------------
+! $Id$
+!--------------------------------------------------------------------------------------------------
+!> @author Luv Sharma, Max-Planck-Institut fŸr Eisenforschung GmbH
+!> @author Pratheek Shanthraj, Max-Planck-Institut fŸr Eisenforschung GmbH
+!> @brief material subroutine incorporating kinematics resulting from opening of cleavage planes
+!> @details to be done
+!--------------------------------------------------------------------------------------------------
+module kinematics_cleavage_opening
+ use prec, only: &
+   pReal, &
+   pInt
+
+ implicit none
+ private
+ integer(pInt),                       dimension(:),           allocatable,         public, protected :: &
+   kinematics_cleavage_opening_sizePostResults, &                                                                !< cumulative size of post results
+   kinematics_cleavage_opening_offset, &                                                                         !< which kinematics is my current damage mechanism?
+   kinematics_cleavage_opening_instance                                                                          !< instance of damage kinematics mechanism
+
+ integer(pInt),                       dimension(:,:),         allocatable, target, public  :: &
+   kinematics_cleavage_opening_sizePostResult                                                                    !< size of each post result output
+
+ character(len=64),                   dimension(:,:),         allocatable, target, public  :: &
+   kinematics_cleavage_opening_output                                                                            !< name of each post result output
+   
+ integer(pInt),                       dimension(:),           allocatable, target, public  :: &
+   kinematics_cleavage_opening_Noutput                                                                           !< number of outputs per instance of this damage 
+
+ integer(pInt),                       dimension(:),           allocatable,         private :: &
+   kinematics_cleavage_opening_totalNcleavage                                                                    !< total number of cleavage systems
+   
+ integer(pInt),                       dimension(:,:),         allocatable,         private :: &
+   kinematics_cleavage_opening_Ncleavage                                                                         !< number of cleavage systems per family
+   
+ real(pReal),                         dimension(:),           allocatable,         private :: &
+   kinematics_cleavage_opening_sdot_0, &
+   kinematics_cleavage_opening_N
+
+ real(pReal),                         dimension(:,:),         allocatable,         private :: &
+   kinematics_cleavage_opening_critDisp, &
+   kinematics_cleavage_opening_critLoad
+
+ public :: &
+   kinematics_cleavage_opening_init, &
+   kinematics_cleavage_opening_LiAndItsTangent
+
+contains
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief module initialization
+!> @details reads in material parameters, allocates arrays, and does sanity checks
+!--------------------------------------------------------------------------------------------------
+subroutine kinematics_cleavage_opening_init(fileUnit)
+ use, intrinsic :: iso_fortran_env                                                                  ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
+ use debug, only: &
+   debug_level,&
+   debug_constitutive,&
+   debug_levelBasic
+ use IO, only: &
+   IO_read, &
+   IO_lc, &
+   IO_getTag, &
+   IO_isBlank, &
+   IO_stringPos, &
+   IO_stringValue, &
+   IO_floatValue, &
+   IO_intValue, &
+   IO_warning, &
+   IO_error, &
+   IO_timeStamp, &
+   IO_EOF
+ use material, only: &
+   phase_kinematics, &
+   phase_Nkinematics, &
+   phase_Noutput, &
+   KINEMATICS_cleavage_opening_label, &
+   KINEMATICS_cleavage_opening_ID, &
+   material_Nphase, &
+   MATERIAL_partPhase
+ use numerics,only: &
+   worldrank
+ use lattice, only: &
+   lattice_maxNcleavageFamily, &
+   lattice_NcleavageSystem
+
+ implicit none
+ integer(pInt), intent(in) :: fileUnit
+
+ integer(pInt), allocatable, dimension(:) :: chunkPos
+ integer(pInt) :: maxNinstance,phase,instance,kinematics
+ integer(pInt) :: Nchunks_CleavageFamilies = 0_pInt, j   
+ character(len=65536) :: &
+   tag  = '', &
+   line = ''
+
+ mainProcess: if (worldrank == 0) then 
+   write(6,'(/,a)')   ' <<<+-  kinematics_'//KINEMATICS_cleavage_opening_LABEL//' init  -+>>>'
+   write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
+#include "compilation_info.f90"
+ endif mainProcess
+
+ maxNinstance = int(count(phase_kinematics == KINEMATICS_cleavage_opening_ID),pInt)
+ if (maxNinstance == 0_pInt) return
+ 
+ if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
+   write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
+ 
+ allocate(kinematics_cleavage_opening_offset(material_Nphase), source=0_pInt)
+ allocate(kinematics_cleavage_opening_instance(material_Nphase), source=0_pInt)
+ do phase = 1, material_Nphase
+   kinematics_cleavage_opening_instance(phase) = count(phase_kinematics(:,1:phase) == kinematics_cleavage_opening_ID)
+   do kinematics = 1, phase_Nkinematics(phase)
+     if (phase_kinematics(kinematics,phase) == kinematics_cleavage_opening_ID) &
+       kinematics_cleavage_opening_offset(phase) = kinematics
+   enddo    
+ enddo
+   
+ allocate(kinematics_cleavage_opening_sizePostResults(maxNinstance),                      source=0_pInt)
+ allocate(kinematics_cleavage_opening_sizePostResult(maxval(phase_Noutput),maxNinstance), source=0_pInt)
+ allocate(kinematics_cleavage_opening_output(maxval(phase_Noutput),maxNinstance))
+          kinematics_cleavage_opening_output = ''
+ allocate(kinematics_cleavage_opening_Noutput(maxNinstance),                              source=0_pInt) 
+ allocate(kinematics_cleavage_opening_critDisp(lattice_maxNcleavageFamily,maxNinstance),  source=0.0_pReal) 
+ allocate(kinematics_cleavage_opening_critLoad(lattice_maxNcleavageFamily,maxNinstance),  source=0.0_pReal) 
+ allocate(kinematics_cleavage_opening_Ncleavage(lattice_maxNcleavageFamily,maxNinstance), source=0_pInt)
+ allocate(kinematics_cleavage_opening_totalNcleavage(maxNinstance),                       source=0_pInt)
+ allocate(kinematics_cleavage_opening_sdot_0(maxNinstance),                               source=0.0_pReal) 
+ allocate(kinematics_cleavage_opening_N(maxNinstance),                                    source=0.0_pReal) 
+
+ rewind(fileUnit)
+ phase = 0_pInt
+ do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase)         ! wind forward to <phase>
+   line = IO_read(fileUnit)
+ enddo
+ 
+ parsingFile: do while (trim(line) /= IO_EOF)                                                       ! read through sections of phase part
+   line = IO_read(fileUnit)
+   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
+   if (IO_getTag(line,'<','>') /= '') then                                                          ! stop at next part
+     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
+     exit                                                                                           
+   endif   
+   if (IO_getTag(line,'[',']') /= '') then                                                          ! next phase section
+     phase = phase + 1_pInt                                                                         ! advance phase section counter
+     cycle                                                                                          ! skip to next line
+   endif
+   if (phase > 0_pInt ) then; if (any(phase_kinematics(:,phase) == KINEMATICS_cleavage_opening_ID)) then         ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
+     instance = kinematics_cleavage_opening_instance(phase)                                                         ! which instance of my damage is present phase
+     chunkPos = IO_stringPos(line)
+     tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt))                                             ! extract key
+     select case(tag)
+       case ('anisobrittle_sdot0')
+         kinematics_cleavage_opening_sdot_0(instance) = IO_floatValue(line,chunkPos,2_pInt)
+         
+       case ('anisobrittle_ratesensitivity')
+         kinematics_cleavage_opening_N(instance) = IO_floatValue(line,chunkPos,2_pInt)
+         
+       case ('ncleavage')  !
+         Nchunks_CleavageFamilies = chunkPos(1) - 1_pInt
+         do j = 1_pInt, Nchunks_CleavageFamilies
+           kinematics_cleavage_opening_Ncleavage(j,instance) = IO_intValue(line,chunkPos,1_pInt+j)
+         enddo
+
+       case ('anisobrittle_criticaldisplacement')
+         do j = 1_pInt, Nchunks_CleavageFamilies
+           kinematics_cleavage_opening_critDisp(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j)
+         enddo
+
+       case ('anisobrittle_criticalload')
+         do j = 1_pInt, Nchunks_CleavageFamilies
+           kinematics_cleavage_opening_critLoad(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j)
+         enddo
+
+     end select
+   endif; endif
+ enddo parsingFile
+
+!--------------------------------------------------------------------------------------------------
+!  sanity checks
+ sanityChecks: do phase = 1_pInt, material_Nphase   
+   myPhase: if (any(phase_kinematics(:,phase) == KINEMATICS_cleavage_opening_ID)) then
+     instance = kinematics_cleavage_opening_instance(phase)
+     kinematics_cleavage_opening_Ncleavage(1:lattice_maxNcleavageFamily,instance) = &
+       min(lattice_NcleavageSystem(1:lattice_maxNcleavageFamily,phase),&                            ! limit active cleavage systems per family to min of available and requested
+           kinematics_cleavage_opening_Ncleavage(1:lattice_maxNcleavageFamily,instance))
+     kinematics_cleavage_opening_totalNcleavage(instance)  = sum(kinematics_cleavage_opening_Ncleavage(:,instance)) ! how many cleavage systems altogether
+     if (kinematics_cleavage_opening_sdot_0(instance) <= 0.0_pReal) &
+       call IO_error(211_pInt,el=instance,ext_msg='sdot_0 ('//KINEMATICS_cleavage_opening_LABEL//')')
+     if (any(kinematics_cleavage_opening_critDisp(1:Nchunks_CleavageFamilies,instance) < 0.0_pReal)) &
+       call IO_error(211_pInt,el=instance,ext_msg='critical_displacement ('//KINEMATICS_cleavage_opening_LABEL//')')
+     if (any(kinematics_cleavage_opening_critLoad(1:Nchunks_CleavageFamilies,instance) < 0.0_pReal)) &
+       call IO_error(211_pInt,el=instance,ext_msg='critical_load ('//KINEMATICS_cleavage_opening_LABEL//')')
+     if (kinematics_cleavage_opening_N(instance) <= 0.0_pReal) &
+       call IO_error(211_pInt,el=instance,ext_msg='rate_sensitivity ('//KINEMATICS_cleavage_opening_LABEL//')')
+   endif myPhase
+ enddo sanityChecks
+ 
+end subroutine kinematics_cleavage_opening_init
+
+!--------------------------------------------------------------------------------------------------
+!> @brief  contains the constitutive equation for calculating the velocity gradient  
+!--------------------------------------------------------------------------------------------------
+subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar3333, Tstar_v, ipc, ip, el)
+ use prec, only: &
+   tol_math_check
+ use material, only: &
+   phaseAt, phasememberAt, &
+   material_homog, &
+   damage, &
+   damageMapping
+ use lattice, only: &
+   lattice_Scleavage, &
+   lattice_Scleavage_v, &
+   lattice_maxNcleavageFamily, &
+   lattice_NcleavageSystem
+ 
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< grain number
+   ip, &                                                                                            !< integration point number
+   el                                                                                               !< element number
+ real(pReal),   intent(in),  dimension(6) :: &
+   Tstar_v                                                                                          !< 2nd Piola-Kirchhoff stress
+ real(pReal),   intent(out), dimension(3,3) :: &
+   Ld                                                                                               !< damage velocity gradient
+ real(pReal),   intent(out), dimension(3,3,3,3) :: &
+   dLd_dTstar3333                                                                                   !< derivative of Ld with respect to Tstar (4th-order tensor)
+ integer(pInt) :: &
+   phase, &
+   constituent, &
+   instance, &
+   homog, damageOffset, &
+   f, i, index_myFamily, k, l, m, n
+ real(pReal) :: &
+   traction_d, traction_t, traction_n, traction_crit, &
+   udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt
+
+ phase = phaseAt(ipc,ip,el)
+ constituent = phasememberAt(ipc,ip,el)
+ instance = kinematics_cleavage_opening_instance(phase)
+ homog = material_homog(ip,el)
+ damageOffset = damageMapping(homog)%p(ip,el)
+ 
+ Ld = 0.0_pReal
+ dLd_dTstar3333 = 0.0_pReal
+ do f = 1_pInt,lattice_maxNcleavageFamily
+   index_myFamily = sum(lattice_NcleavageSystem(1:f-1_pInt,phase))                                   ! at which index starts my family
+   do i = 1_pInt,kinematics_cleavage_opening_Ncleavage(f,instance)                                            ! process each (active) cleavage system in family
+     traction_d    = dot_product(Tstar_v,lattice_Scleavage_v(1:6,1,index_myFamily+i,phase))
+     traction_t    = dot_product(Tstar_v,lattice_Scleavage_v(1:6,2,index_myFamily+i,phase))
+     traction_n    = dot_product(Tstar_v,lattice_Scleavage_v(1:6,3,index_myFamily+i,phase))
+     traction_crit = kinematics_cleavage_opening_critLoad(f,instance)* &
+                     damage(homog)%p(damageOffset)*damage(homog)%p(damageOffset)
+     udotd = &
+       sign(1.0_pReal,traction_d)* &
+       kinematics_cleavage_opening_sdot_0(instance)* &
+       (max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**kinematics_cleavage_opening_N(instance)
+     if (abs(udotd) > tol_math_check) then
+       Ld = Ld + udotd*lattice_Scleavage(1:3,1:3,1,index_myFamily+i,phase)
+       dudotd_dt = sign(1.0_pReal,traction_d)*udotd*kinematics_cleavage_opening_N(instance)/ &
+                   max(0.0_pReal, abs(traction_d) - traction_crit)
+       forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
+         dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + &
+           dudotd_dt*lattice_Scleavage(k,l,1,index_myFamily+i,phase)* &
+                     lattice_Scleavage(m,n,1,index_myFamily+i,phase)
+     endif                
+
+     udott = &
+       sign(1.0_pReal,traction_t)* &
+       kinematics_cleavage_opening_sdot_0(instance)* &
+       (max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**kinematics_cleavage_opening_N(instance)
+     if (abs(udott) > tol_math_check) then
+       Ld = Ld + udott*lattice_Scleavage(1:3,1:3,2,index_myFamily+i,phase)
+       dudott_dt = sign(1.0_pReal,traction_t)*udott*kinematics_cleavage_opening_N(instance)/ &
+                   max(0.0_pReal, abs(traction_t) - traction_crit)  
+       forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
+         dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + &
+           dudott_dt*lattice_Scleavage(k,l,2,index_myFamily+i,phase)* &
+                     lattice_Scleavage(m,n,2,index_myFamily+i,phase)
+     endif                
+
+     udotn = &
+       sign(1.0_pReal,traction_n)* &
+       kinematics_cleavage_opening_sdot_0(instance)* &
+       (max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**kinematics_cleavage_opening_N(instance)
+     if (abs(udotn) > tol_math_check) then
+       Ld = Ld + udotn*lattice_Scleavage(1:3,1:3,3,index_myFamily+i,phase)
+       dudotn_dt = sign(1.0_pReal,traction_n)*udotn*kinematics_cleavage_opening_N(instance)/ &
+                   max(0.0_pReal, abs(traction_n) - traction_crit)
+       forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
+         dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + &
+           dudotn_dt*lattice_Scleavage(k,l,3,index_myFamily+i,phase)* &
+                     lattice_Scleavage(m,n,3,index_myFamily+i,phase)
+     endif                
+
+   enddo
+ enddo
+ 
+end subroutine kinematics_cleavage_opening_LiAndItsTangent
+
+end module kinematics_cleavage_opening

code/kinematics_hydrogen_strain.f90

@@ -0,0 +1,264 @@
+!--------------------------------------------------------------------------------------------------
+! $Id$
+!--------------------------------------------------------------------------------------------------
+!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
+!> @brief material subroutine incorporating kinematics resulting from interstitial hydrogen
+!> @details to be done
+!--------------------------------------------------------------------------------------------------
+module kinematics_hydrogen_strain
+ use prec, only: &
+   pReal, &
+   pInt
+
+ implicit none
+ private
+ integer(pInt),                       dimension(:),           allocatable,         public, protected :: &
+   kinematics_hydrogen_strain_sizePostResults, &                                                                !< cumulative size of post results
+   kinematics_hydrogen_strain_offset, &                                                                         !< which kinematics is my current damage mechanism?
+   kinematics_hydrogen_strain_instance                                                                          !< instance of damage kinematics mechanism
+
+ integer(pInt),                       dimension(:,:),         allocatable, target, public  :: &
+   kinematics_hydrogen_strain_sizePostResult                                                                    !< size of each post result output
+
+ character(len=64),                   dimension(:,:),         allocatable, target, public  :: &
+   kinematics_hydrogen_strain_output                                                                            !< name of each post result output
+   
+ integer(pInt),                       dimension(:),           allocatable, target, public  :: &
+   kinematics_hydrogen_strain_Noutput                                                                           !< number of outputs per instance of this damage 
+
+ real(pReal),                         dimension(:),           allocatable,         private :: &
+   kinematics_hydrogen_strain_coeff
+
+ public :: &
+   kinematics_hydrogen_strain_init, &
+   kinematics_hydrogen_strain_initialStrain, &
+   kinematics_hydrogen_strain_LiAndItsTangent, &
+   kinematics_hydrogen_strain_ChemPotAndItsTangent
+
+contains
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief module initialization
+!> @details reads in material parameters, allocates arrays, and does sanity checks
+!--------------------------------------------------------------------------------------------------
+subroutine kinematics_hydrogen_strain_init(fileUnit)
+ use, intrinsic :: iso_fortran_env                                                                  ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
+ use debug, only: &
+   debug_level,&
+   debug_constitutive,&
+   debug_levelBasic
+ use IO, only: &
+   IO_read, &
+   IO_lc, &
+   IO_getTag, &
+   IO_isBlank, &
+   IO_stringPos, &
+   IO_stringValue, &
+   IO_floatValue, &
+   IO_intValue, &
+   IO_warning, &
+   IO_error, &
+   IO_timeStamp, &
+   IO_EOF
+ use material, only: &
+   phase_kinematics, &
+   phase_Nkinematics, &
+   phase_Noutput, &
+   KINEMATICS_hydrogen_strain_label, &
+   KINEMATICS_hydrogen_strain_ID, &
+   material_Nphase, &
+   MATERIAL_partPhase
+ use numerics,only: &
+   worldrank
+
+ implicit none
+ integer(pInt), intent(in) :: fileUnit
+
+ integer(pInt), allocatable, dimension(:) :: chunkPos
+ integer(pInt) :: maxNinstance,phase,instance,kinematics
+ character(len=65536) :: &
+   tag  = '', &
+   line = ''
+
+ mainProcess: if (worldrank == 0) then 
+   write(6,'(/,a)')   ' <<<+-  kinematics_'//KINEMATICS_hydrogen_strain_LABEL//' init  -+>>>'
+   write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
+#include "compilation_info.f90"
+ endif mainProcess
+
+ maxNinstance = int(count(phase_kinematics == KINEMATICS_hydrogen_strain_ID),pInt)
+ if (maxNinstance == 0_pInt) return
+ 
+ if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
+   write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
+ 
+ allocate(kinematics_hydrogen_strain_offset(material_Nphase), source=0_pInt)
+ allocate(kinematics_hydrogen_strain_instance(material_Nphase), source=0_pInt)
+ do phase = 1, material_Nphase
+   kinematics_hydrogen_strain_instance(phase) = count(phase_kinematics(:,1:phase) == kinematics_hydrogen_strain_ID)
+   do kinematics = 1, phase_Nkinematics(phase)
+     if (phase_kinematics(kinematics,phase) == kinematics_hydrogen_strain_ID) &
+       kinematics_hydrogen_strain_offset(phase) = kinematics
+   enddo    
+ enddo
+   
+ allocate(kinematics_hydrogen_strain_sizePostResults(maxNinstance),                     source=0_pInt)
+ allocate(kinematics_hydrogen_strain_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
+ allocate(kinematics_hydrogen_strain_output(maxval(phase_Noutput),maxNinstance))
+          kinematics_hydrogen_strain_output = ''
+ allocate(kinematics_hydrogen_strain_Noutput(maxNinstance),                             source=0_pInt) 
+ allocate(kinematics_hydrogen_strain_coeff(maxNinstance),                               source=0.0_pReal) 
+
+ rewind(fileUnit)
+ phase = 0_pInt
+ do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase)         ! wind forward to <phase>
+   line = IO_read(fileUnit)
+ enddo
+ 
+ parsingFile: do while (trim(line) /= IO_EOF)                                                       ! read through sections of phase part
+   line = IO_read(fileUnit)
+   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
+   if (IO_getTag(line,'<','>') /= '') then                                                          ! stop at next part
+     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
+     exit                                                                                           
+   endif   
+   if (IO_getTag(line,'[',']') /= '') then                                                          ! next phase section
+     phase = phase + 1_pInt                                                                         ! advance phase section counter
+     cycle                                                                                          ! skip to next line
+   endif
+   if (phase > 0_pInt ) then; if (any(phase_kinematics(:,phase) == KINEMATICS_hydrogen_strain_ID)) then         ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
+     instance = kinematics_hydrogen_strain_instance(phase)                                                         ! which instance of my damage is present phase
+     chunkPos = IO_stringPos(line)
+     tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt))                                             ! extract key
+     select case(tag)
+       case ('hydrogen_strain_coeff')
+         kinematics_hydrogen_strain_coeff(instance) = IO_floatValue(line,chunkPos,2_pInt)
+
+     end select
+   endif; endif
+ enddo parsingFile
+
+end subroutine kinematics_hydrogen_strain_init
+
+!--------------------------------------------------------------------------------------------------
+!> @brief  report initial hydrogen strain based on current hydrogen conc deviation from 
+!>         equillibrium (0)
+!--------------------------------------------------------------------------------------------------
+pure function kinematics_hydrogen_strain_initialStrain(ipc, ip, el)
+ use math, only: &
+   math_I3
+ use material, only: &
+   material_phase, &
+   material_homog, &
+   hydrogenConc, &
+   hydrogenfluxMapping
+ use lattice, only: &
+   lattice_equilibriumHydrogenConcentration
+ 
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< grain number
+   ip, &                                                                                            !< integration point number
+   el                                                                                               !< element number
+ real(pReal), dimension(3,3) :: &
+   kinematics_hydrogen_strain_initialStrain                                                       !< initial thermal strain (should be small strain, though)
+ integer(pInt) :: &
+   phase, &
+   homog, offset, instance
+   
+ phase = material_phase(ipc,ip,el)
+ instance = kinematics_hydrogen_strain_instance(phase)
+ homog = material_homog(ip,el)
+ offset = hydrogenfluxMapping(homog)%p(ip,el)
+ 
+ kinematics_hydrogen_strain_initialStrain = &
+   (hydrogenConc(homog)%p(offset) - lattice_equilibriumHydrogenConcentration(phase)) * &
+   kinematics_hydrogen_strain_coeff(instance)* math_I3
+  
+end function kinematics_hydrogen_strain_initialStrain
+
+!--------------------------------------------------------------------------------------------------
+!> @brief  contains the constitutive equation for calculating the velocity gradient  
+!--------------------------------------------------------------------------------------------------
+subroutine kinematics_hydrogen_strain_LiAndItsTangent(Li, dLi_dTstar3333, ipc, ip, el)
+ use material, only: &
+   material_phase, &
+   material_homog, &
+   hydrogenConc, &
+   hydrogenConcRate, &
+   hydrogenfluxMapping
+ use math, only: &
+   math_I3
+ use lattice, only: &
+   lattice_equilibriumHydrogenConcentration
+ 
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< grain number
+   ip, &                                                                                            !< integration point number
+   el                                                                                               !< element number
+ real(pReal),   intent(out), dimension(3,3) :: &
+   Li                                                                                               !< thermal velocity gradient
+ real(pReal),   intent(out), dimension(3,3,3,3) :: &
+   dLi_dTstar3333                                                                                   !< derivative of Li with respect to Tstar (4th-order tensor)
+ integer(pInt) :: &
+   phase, &
+   instance, &
+   homog, offset
+ real(pReal) :: &
+   Ch, ChEq, ChDot  
+   
+ phase = material_phase(ipc,ip,el)
+ instance = kinematics_hydrogen_strain_instance(phase)
+ homog = material_homog(ip,el)
+ offset = hydrogenfluxMapping(homog)%p(ip,el)
+ Ch = hydrogenConc(homog)%p(offset)
+ ChDot = hydrogenConcRate(homog)%p(offset)
+ ChEq = lattice_equilibriumHydrogenConcentration(phase)
+ 
+ Li = ChDot*math_I3* &
+      kinematics_hydrogen_strain_coeff(instance)/ &
+      (1.0_pReal + kinematics_hydrogen_strain_coeff(instance)*(Ch - ChEq))
+ dLi_dTstar3333 = 0.0_pReal
+  
+end subroutine kinematics_hydrogen_strain_LiAndItsTangent
+
+!--------------------------------------------------------------------------------------------------
+!> @brief  contains the kinematic contribution to hydrogen chemical potential  
+!--------------------------------------------------------------------------------------------------
+subroutine kinematics_hydrogen_strain_ChemPotAndItsTangent(ChemPot, dChemPot_dCh, Tstar_v, Fi0, Fi, ipc, ip, el)
+ use material, only: &
+   material_phase
+ use math, only: &
+   math_inv33, &
+   math_mul33x33, &
+   math_Mandel6to33, &
+   math_transpose33
+ 
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< grain number
+   ip, &                                                                                            !< integration point number
+   el                                                                                               !< element number
+ real(pReal),  intent(in), dimension(6) :: &
+   Tstar_v
+ real(pReal),  intent(in), dimension(3,3) :: &
+   Fi0, Fi
+ real(pReal),   intent(out) :: &
+   ChemPot, dChemPot_dCh                                                               
+ integer(pInt) :: &
+   phase, &
+   instance
+   
+ phase = material_phase(ipc,ip,el)
+ instance = kinematics_hydrogen_strain_instance(phase)
+ 
+ ChemPot = -kinematics_hydrogen_strain_coeff(instance)* &
+            sum(math_mul33x33(Fi,math_Mandel6to33(Tstar_v))* &
+                math_mul33x33(math_mul33x33(Fi,math_inv33(Fi0)),Fi))
+ dChemPot_dCh = 0.0_pReal
+   
+end subroutine kinematics_hydrogen_strain_ChemPotAndItsTangent
+
+end module kinematics_hydrogen_strain

code/kinematics_slipplane_opening.f90

@@ -0,0 +1,323 @@
+!--------------------------------------------------------------------------------------------------
+! $Id$
+!--------------------------------------------------------------------------------------------------
+!> @author Luv Sharma, Max-Planck-Institut für Eisenforschung GmbH
+!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
+!> @brief material subroutine incorporating kinematics resulting from opening of slip planes
+!> @details to be done
+!--------------------------------------------------------------------------------------------------
+module kinematics_slipplane_opening
+ use prec, only: &
+   pReal, &
+   pInt
+
+ implicit none
+ private
+ integer(pInt),                       dimension(:),           allocatable,         public, protected :: &
+   kinematics_slipplane_opening_sizePostResults, &                                                            !< cumulative size of post results
+   kinematics_slipplane_opening_offset, &                                                                     !< which kinematics is my current damage mechanism?
+   kinematics_slipplane_opening_instance                                                                      !< instance of damage kinematics mechanism
+
+ integer(pInt),                       dimension(:,:),         allocatable, target, public  :: &
+   kinematics_slipplane_opening_sizePostResult                                                                !< size of each post result output
+
+ character(len=64),                   dimension(:,:),         allocatable, target, public  :: &
+   kinematics_slipplane_opening_output                                                                        !< name of each post result output
+   
+ integer(pInt),                       dimension(:),           allocatable, target, public  :: &
+   kinematics_slipplane_opening_Noutput                                                                       !< number of outputs per instance of this damage 
+   
+ integer(pInt),                       dimension(:),           allocatable,         private :: &
+   kinematics_slipplane_opening_totalNslip                                                                    !< total number of slip systems
+
+ integer(pInt),                       dimension(:,:),         allocatable,         private :: &
+   kinematics_slipplane_opening_Nslip                                                                         !< number of slip systems per family
+   
+ real(pReal),                         dimension(:),           allocatable,         private :: &
+   kinematics_slipplane_opening_sdot_0, &
+   kinematics_slipplane_opening_N
+
+ real(pReal),                         dimension(:,:),         allocatable,         private :: &
+   kinematics_slipplane_opening_critPlasticStrain, &
+   kinematics_slipplane_opening_critLoad
+   
+ public :: &
+   kinematics_slipplane_opening_init, &
+   kinematics_slipplane_opening_LiAndItsTangent
+
+contains
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief module initialization
+!> @details reads in material parameters, allocates arrays, and does sanity checks
+!--------------------------------------------------------------------------------------------------
+subroutine kinematics_slipplane_opening_init(fileUnit)
+ use, intrinsic :: iso_fortran_env                                                                  ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
+ use debug, only: &
+   debug_level,&
+   debug_constitutive,&
+   debug_levelBasic
+ use IO, only: &
+   IO_read, &
+   IO_lc, &
+   IO_getTag, &
+   IO_isBlank, &
+   IO_stringPos, &
+   IO_stringValue, &
+   IO_floatValue, &
+   IO_intValue, &
+   IO_warning, &
+   IO_error, &
+   IO_timeStamp, &
+   IO_EOF
+ use material, only: &
+   phase_kinematics, &
+   phase_Nkinematics, &
+   phase_Noutput, &
+   KINEMATICS_slipplane_opening_label, &
+   KINEMATICS_slipplane_opening_ID, &
+   material_Nphase, &
+   MATERIAL_partPhase
+ use numerics,only: &
+   worldrank
+ use lattice, only: &
+   lattice_maxNslipFamily, &
+   lattice_NslipSystem
+
+ implicit none
+ integer(pInt), intent(in) :: fileUnit
+
+ integer(pInt), allocatable, dimension(:) :: chunkPos
+ integer(pInt) :: maxNinstance,phase,instance,kinematics
+ integer(pInt) :: Nchunks_SlipFamilies = 0_pInt, j   
+ character(len=65536) :: &
+   tag  = '', &
+   line = ''
+
+ mainProcess: if (worldrank == 0) then 
+   write(6,'(/,a)')   ' <<<+-  kinematics_'//KINEMATICS_slipplane_opening_LABEL//' init  -+>>>'
+   write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
+#include "compilation_info.f90"
+ endif mainProcess
+
+ maxNinstance = int(count(phase_kinematics == KINEMATICS_slipplane_opening_ID),pInt)
+ if (maxNinstance == 0_pInt) return
+ 
+ if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
+   write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
+ 
+ allocate(kinematics_slipplane_opening_offset(material_Nphase), source=0_pInt)
+ allocate(kinematics_slipplane_opening_instance(material_Nphase), source=0_pInt)
+ do phase = 1, material_Nphase
+   kinematics_slipplane_opening_instance(phase) = count(phase_kinematics(:,1:phase) == kinematics_slipplane_opening_ID)
+   do kinematics = 1, phase_Nkinematics(phase)
+     if (phase_kinematics(kinematics,phase) == kinematics_slipplane_opening_ID) &
+       kinematics_slipplane_opening_offset(phase) = kinematics
+   enddo    
+ enddo
+   
+ allocate(kinematics_slipplane_opening_sizePostResults(maxNinstance),                     source=0_pInt)
+ allocate(kinematics_slipplane_opening_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
+ allocate(kinematics_slipplane_opening_output(maxval(phase_Noutput),maxNinstance))
+          kinematics_slipplane_opening_output = ''
+ allocate(kinematics_slipplane_opening_Noutput(maxNinstance),                             source=0_pInt) 
+ allocate(kinematics_slipplane_opening_critLoad(lattice_maxNslipFamily,maxNinstance),     source=0.0_pReal) 
+ allocate(kinematics_slipplane_opening_critPlasticStrain(lattice_maxNslipFamily,maxNinstance),source=0.0_pReal) 
+ allocate(kinematics_slipplane_opening_Nslip(lattice_maxNslipFamily,maxNinstance),        source=0_pInt)
+ allocate(kinematics_slipplane_opening_totalNslip(maxNinstance),                          source=0_pInt)
+ allocate(kinematics_slipplane_opening_N(maxNinstance),                                   source=0.0_pReal) 
+ allocate(kinematics_slipplane_opening_sdot_0(maxNinstance),                              source=0.0_pReal) 
+
+ rewind(fileUnit)
+ phase = 0_pInt
+ do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase)         ! wind forward to <phase>
+   line = IO_read(fileUnit)
+ enddo
+ 
+ parsingFile: do while (trim(line) /= IO_EOF)                                                       ! read through sections of phase part
+   line = IO_read(fileUnit)
+   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
+   if (IO_getTag(line,'<','>') /= '') then                                                          ! stop at next part
+     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
+     exit                                                                                           
+   endif   
+   if (IO_getTag(line,'[',']') /= '') then                                                          ! next phase section
+     phase = phase + 1_pInt                                                                         ! advance phase section counter
+     cycle                                                                                          ! skip to next line
+   endif
+   if (phase > 0_pInt ) then; if (any(phase_kinematics(:,phase) == KINEMATICS_slipplane_opening_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
+     instance = kinematics_slipplane_opening_instance(phase)                                        ! which instance of my damage is present phase
+     chunkPos = IO_stringPos(line)
+     tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt))                                              ! extract key
+     select case(tag)
+       case ('nslip')  !
+         Nchunks_SlipFamilies = chunkPos(1) - 1_pInt
+         do j = 1_pInt, Nchunks_SlipFamilies
+           kinematics_slipplane_opening_Nslip(j,instance) = IO_intValue(line,chunkPos,1_pInt+j)
+         enddo
+
+       case ('anisoductile_sdot0')
+         kinematics_slipplane_opening_sdot_0(instance) = IO_floatValue(line,chunkPos,2_pInt)
+         
+       case ('anisoductile_criticalplasticstrain')
+         do j = 1_pInt, Nchunks_SlipFamilies
+           kinematics_slipplane_opening_critPlasticStrain(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j)
+         enddo
+         
+       case ('anisoductile_ratesensitivity')
+         kinematics_slipplane_opening_N(instance) = IO_floatValue(line,chunkPos,2_pInt)
+
+       case ('anisoductile_criticalload')
+         do j = 1_pInt, Nchunks_SlipFamilies
+           kinematics_slipplane_opening_critLoad(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j)
+         enddo
+         
+     end select
+   endif; endif
+ enddo parsingFile
+
+!--------------------------------------------------------------------------------------------------
+!  sanity checks
+ sanityChecks: do phase = 1_pInt, material_Nphase   
+   myPhase: if (any(phase_kinematics(:,phase) == KINEMATICS_slipplane_opening_ID)) then
+     instance = kinematics_slipplane_opening_instance(phase)
+     kinematics_slipplane_opening_Nslip(1:lattice_maxNslipFamily,instance) = &
+       min(lattice_NslipSystem(1:lattice_maxNslipFamily,phase),&                                    ! limit active cleavage systems per family to min of available and requested
+           kinematics_slipplane_opening_Nslip(1:lattice_maxNslipFamily,instance))
+         kinematics_slipplane_opening_totalNslip(instance) = sum(kinematics_slipplane_opening_Nslip(:,instance))
+     if (kinematics_slipplane_opening_sdot_0(instance) <= 0.0_pReal) &
+       call IO_error(211_pInt,el=instance,ext_msg='sdot_0 ('//KINEMATICS_slipplane_opening_LABEL//')')
+     if (any(kinematics_slipplane_opening_critPlasticStrain(:,instance) < 0.0_pReal)) &
+       call IO_error(211_pInt,el=instance,ext_msg='criticaPlasticStrain ('//KINEMATICS_slipplane_opening_LABEL//')')
+     if (kinematics_slipplane_opening_N(instance) <= 0.0_pReal) &
+       call IO_error(211_pInt,el=instance,ext_msg='rate_sensitivity ('//KINEMATICS_slipplane_opening_LABEL//')')
+   endif myPhase
+ enddo sanityChecks
+  
+
+end subroutine kinematics_slipplane_opening_init
+
+!--------------------------------------------------------------------------------------------------
+!> @brief  contains the constitutive equation for calculating the velocity gradient  
+!--------------------------------------------------------------------------------------------------
+subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar3333, Tstar_v, ipc, ip, el)
+ use prec, only: &
+   tol_math_check
+ use lattice, only: &
+   lattice_maxNslipFamily, &
+   lattice_NslipSystem, &
+   lattice_sd, &
+   lattice_st, &
+   lattice_sn
+ use material, only: &
+   phaseAt, phasememberAt, &
+   material_homog, &
+   damage, &
+   damageMapping
+ use math, only: &
+   math_Plain3333to99, &
+   math_I3, &
+   math_identity4th, &
+   math_symmetric33, &
+   math_Mandel33to6, &
+   math_tensorproduct33, &
+   math_det33, &
+   math_mul33x33
+ 
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< grain number
+   ip, &                                                                                            !< integration point number
+   el                                                                                               !< element number
+ real(pReal),   intent(in),  dimension(6) :: &
+   Tstar_v                                                                                          !< 2nd Piola-Kirchhoff stress
+ real(pReal),   intent(out), dimension(3,3) :: &
+   Ld                                                                                               !< damage velocity gradient
+ real(pReal),   intent(out), dimension(3,3,3,3) :: &
+   dLd_dTstar3333                                                                                   !< derivative of Ld with respect to Tstar (4th-order tensor)
+ real(pReal),   dimension(3,3) :: &
+   projection_d, projection_t, projection_n                                                         !< projection modes 3x3 tensor
+ real(pReal),   dimension(6) :: &
+   projection_d_v, projection_t_v, projection_n_v                                                   !< projection modes 3x3 vector
+ integer(pInt) :: &
+   phase, &
+   constituent, &
+   instance, &
+   homog, damageOffset, &
+   f, i, index_myFamily, k, l, m, n
+ real(pReal) :: &
+   traction_d, traction_t, traction_n, traction_crit, &
+   udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt
+   
+ phase = phaseAt(ipc,ip,el)
+ constituent = phasememberAt(ipc,ip,el)
+ instance = kinematics_slipplane_opening_instance(phase)
+ homog = material_homog(ip,el)
+ damageOffset = damageMapping(homog)%p(ip,el)
+ 
+ Ld = 0.0_pReal
+ dLd_dTstar3333 = 0.0_pReal
+ do f = 1_pInt,lattice_maxNslipFamily
+   index_myFamily = sum(lattice_NslipSystem(1:f-1_pInt,phase))                                      ! at which index starts my family
+   do i = 1_pInt,kinematics_slipplane_opening_Nslip(f,instance)                                              ! process each (active) slip system in family
+     projection_d = math_tensorproduct33(lattice_sd(1:3,index_myFamily+i,phase),&
+                                       lattice_sn(1:3,index_myFamily+i,phase))
+     projection_t = math_tensorproduct33(lattice_st(1:3,index_myFamily+i,phase),&
+                                       lattice_sn(1:3,index_myFamily+i,phase))
+     projection_n = math_tensorproduct33(lattice_sn(1:3,index_myFamily+i,phase),&
+                                       lattice_sn(1:3,index_myFamily+i,phase))
+
+     projection_d_v(1:6) = math_Mandel33to6(math_symmetric33(projection_d(1:3,1:3)))
+     projection_t_v(1:6) = math_Mandel33to6(math_symmetric33(projection_t(1:3,1:3)))
+     projection_n_v(1:6) = math_Mandel33to6(math_symmetric33(projection_n(1:3,1:3)))
+   
+     traction_d    = dot_product(Tstar_v,projection_d_v(1:6))
+     traction_t    = dot_product(Tstar_v,projection_t_v(1:6))
+     traction_n    = dot_product(Tstar_v,projection_n_v(1:6))
+     
+     traction_crit = kinematics_slipplane_opening_critLoad(f,instance)* &
+                     damage(homog)%p(damageOffset)                                                        ! degrading critical load carrying capacity by damage 
+
+     udotd = &
+       sign(1.0_pReal,traction_d)* &
+       kinematics_slipplane_opening_sdot_0(instance)* &
+       (abs(traction_d)/traction_crit - &
+        abs(traction_d)/kinematics_slipplane_opening_critLoad(f,instance))**kinematics_slipplane_opening_N(instance)
+     if (abs(udotd) > tol_math_check) then
+       Ld = Ld + udotd*projection_d
+       dudotd_dt = udotd*kinematics_slipplane_opening_N(instance)/traction_d
+       forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
+         dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + &
+           dudotd_dt*projection_d(k,l)*projection_d(m,n)
+     endif                
+ 
+     udott = &
+       sign(1.0_pReal,traction_t)* &
+       kinematics_slipplane_opening_sdot_0(instance)* &
+       (abs(traction_t)/traction_crit - &
+        abs(traction_t)/kinematics_slipplane_opening_critLoad(f,instance))**kinematics_slipplane_opening_N(instance)
+     if (abs(udott) > tol_math_check) then
+       Ld = Ld + udott*projection_t
+       dudott_dt = udott*kinematics_slipplane_opening_N(instance)/traction_t
+       forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
+         dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + &
+           dudott_dt*projection_t(k,l)*projection_t(m,n)
+     endif
+     udotn = &
+       kinematics_slipplane_opening_sdot_0(instance)* &
+       (max(0.0_pReal,traction_n)/traction_crit - &
+        max(0.0_pReal,traction_n)/kinematics_slipplane_opening_critLoad(f,instance))**kinematics_slipplane_opening_N(instance)
+     if (abs(udotn) > tol_math_check) then
+       Ld = Ld + udotn*projection_n
+       dudotn_dt = udotn*kinematics_slipplane_opening_N(instance)/traction_n
+       forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
+         dLd_dTstar3333(k,l,m,n) = dLd_dTstar3333(k,l,m,n) + &
+           dudotn_dt*projection_n(k,l)*projection_n(m,n)
+     endif 
+   enddo
+ enddo
+  
+end subroutine kinematics_slipplane_opening_LiAndItsTangent
+
+end module kinematics_slipplane_opening

code/kinematics_thermal_expansion.f90

@@ -0,0 +1,228 @@
+!--------------------------------------------------------------------------------------------------
+! $Id$
+!--------------------------------------------------------------------------------------------------
+!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
+!> @brief material subroutine incorporating kinematics resulting from thermal expansion
+!> @details to be done
+!--------------------------------------------------------------------------------------------------
+module kinematics_thermal_expansion
+ use prec, only: &
+   pReal, &
+   pInt
+
+ implicit none
+ private
+ integer(pInt),                       dimension(:),           allocatable,         public, protected :: &
+   kinematics_thermal_expansion_sizePostResults, &                                                                !< cumulative size of post results
+   kinematics_thermal_expansion_offset, &                                                                         !< which kinematics is my current damage mechanism?
+   kinematics_thermal_expansion_instance                                                                          !< instance of damage kinematics mechanism
+
+ integer(pInt),                       dimension(:,:),         allocatable, target, public  :: &
+   kinematics_thermal_expansion_sizePostResult                                                                    !< size of each post result output
+
+ character(len=64),                   dimension(:,:),         allocatable, target, public  :: &
+   kinematics_thermal_expansion_output                                                                            !< name of each post result output
+   
+ integer(pInt),                       dimension(:),           allocatable, target, public  :: &
+   kinematics_thermal_expansion_Noutput                                                                           !< number of outputs per instance of this damage 
+
+! enum, bind(c)                                                                                                   ! ToDo kinematics need state machinery to deal with sizePostResult
+!   enumerator :: undefined_ID, &                                                                                 ! possible remedy is to decouple having state vars from having output
+!                 thermalexpansionrate_ID                                                                         ! which means to separate user-defined types tState + tOutput...
+! end enum
+ public :: &
+   kinematics_thermal_expansion_init, &
+   kinematics_thermal_expansion_initialStrain, &
+   kinematics_thermal_expansion_LiAndItsTangent
+
+contains
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief module initialization
+!> @details reads in material parameters, allocates arrays, and does sanity checks
+!--------------------------------------------------------------------------------------------------
+subroutine kinematics_thermal_expansion_init(fileUnit)
+ use, intrinsic :: iso_fortran_env                                                                  ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
+ use debug, only: &
+   debug_level,&
+   debug_constitutive,&
+   debug_levelBasic
+ use IO, only: &
+   IO_read, &
+   IO_lc, &
+   IO_getTag, &
+   IO_isBlank, &
+   IO_stringPos, &
+   IO_stringValue, &
+   IO_floatValue, &
+   IO_intValue, &
+   IO_warning, &
+   IO_error, &
+   IO_timeStamp, &
+   IO_EOF
+ use material, only: &
+   phase_kinematics, &
+   phase_Nkinematics, &
+   phase_Noutput, &
+   KINEMATICS_thermal_expansion_label, &
+   KINEMATICS_thermal_expansion_ID, &
+   material_Nphase, &
+   MATERIAL_partPhase
+ use numerics,only: &
+   worldrank
+
+ implicit none
+ integer(pInt), intent(in) :: fileUnit
+
+ integer(pInt), allocatable, dimension(:) :: chunkPos
+ integer(pInt) :: maxNinstance,phase,instance,kinematics
+ character(len=65536) :: &
+   tag     = '', &
+   output  = '', &
+   line    = ''
+
+ mainProcess: if (worldrank == 0) then 
+   write(6,'(/,a)')   ' <<<+-  kinematics_'//KINEMATICS_thermal_expansion_LABEL//' init  -+>>>'
+   write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
+#include "compilation_info.f90"
+ endif mainProcess
+
+ maxNinstance = int(count(phase_kinematics == KINEMATICS_thermal_expansion_ID),pInt)
+ if (maxNinstance == 0_pInt) return
+ 
+ if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
+   write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
+ 
+ allocate(kinematics_thermal_expansion_offset(material_Nphase), source=0_pInt)
+ allocate(kinematics_thermal_expansion_instance(material_Nphase), source=0_pInt)
+ do phase = 1, material_Nphase
+   kinematics_thermal_expansion_instance(phase) = count(phase_kinematics(:,1:phase) == kinematics_thermal_expansion_ID)
+   do kinematics = 1, phase_Nkinematics(phase)
+     if (phase_kinematics(kinematics,phase) == kinematics_thermal_expansion_ID) &
+       kinematics_thermal_expansion_offset(phase) = kinematics
+   enddo    
+ enddo
+   
+ allocate(kinematics_thermal_expansion_sizePostResults(maxNinstance),                     source=0_pInt)
+ allocate(kinematics_thermal_expansion_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
+ allocate(kinematics_thermal_expansion_output(maxval(phase_Noutput),maxNinstance))
+          kinematics_thermal_expansion_output = ''
+ allocate(kinematics_thermal_expansion_Noutput(maxNinstance),                             source=0_pInt) 
+
+ rewind(fileUnit)
+ phase = 0_pInt
+ do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase)         ! wind forward to <phase>
+   line = IO_read(fileUnit)
+ enddo
+ 
+ parsingFile: do while (trim(line) /= IO_EOF)                                                       ! read through sections of phase part
+   line = IO_read(fileUnit)
+   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
+   if (IO_getTag(line,'<','>') /= '') then                                                          ! stop at next part
+     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
+     exit                                                                                           
+   endif   
+   if (IO_getTag(line,'[',']') /= '') then                                                          ! next phase section
+     phase = phase + 1_pInt                                                                         ! advance phase section counter
+     cycle                                                                                          ! skip to next line
+   endif
+   if (phase > 0_pInt ) then; if (any(phase_kinematics(:,phase) == KINEMATICS_thermal_expansion_ID)) then         ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
+     instance = kinematics_thermal_expansion_instance(phase)                                                         ! which instance of my damage is present phase
+     chunkPos = IO_stringPos(line)
+     tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt))                                             ! extract key...
+     select case(tag)
+!       case ('(output)')
+!         output = IO_lc(IO_stringValue(line,chunkPos,2_pInt))                                       ! ...and corresponding output
+!         select case(output)
+!           case ('thermalexpansionrate')
+!             kinematics_thermal_expansion_Noutput(instance) = kinematics_thermal_expansion_Noutput(instance) + 1_pInt
+!             kinematics_thermal_expansion_outputID(kinematics_thermal_expansion_Noutput(instance),instance) = &
+!               thermalexpansionrate_ID
+!             kinematics_thermal_expansion_output(kinematics_thermal_expansion_Noutput(instance),instance) = output
+! ToDo add sizePostResult loop afterwards...
+
+     end select
+   endif; endif
+ enddo parsingFile
+
+end subroutine kinematics_thermal_expansion_init
+
+!--------------------------------------------------------------------------------------------------
+!> @brief  report initial thermal strain based on current temperature deviation from reference
+!--------------------------------------------------------------------------------------------------
+pure function kinematics_thermal_expansion_initialStrain(ipc, ip, el)
+ use material, only: &
+   material_phase, &
+   material_homog, &
+   temperature, &
+   thermalMapping
+ use lattice, only: &
+   lattice_thermalExpansion33, &
+   lattice_referenceTemperature
+ 
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< grain number
+   ip, &                                                                                            !< integration point number
+   el                                                                                               !< element number
+ real(pReal), dimension(3,3) :: &
+   kinematics_thermal_expansion_initialStrain                                                       !< initial thermal strain (should be small strain, though)
+ integer(pInt) :: &
+   phase, &
+   homog, offset
+   
+ phase = material_phase(ipc,ip,el)
+ homog = material_homog(ip,el)
+ offset = thermalMapping(homog)%p(ip,el)
+ 
+ kinematics_thermal_expansion_initialStrain = &
+   (temperature(homog)%p(offset) - lattice_referenceTemperature(phase)) * &
+   lattice_thermalExpansion33(1:3,1:3,phase)
+  
+end function kinematics_thermal_expansion_initialStrain
+
+!--------------------------------------------------------------------------------------------------
+!> @brief  contains the constitutive equation for calculating the velocity gradient  
+!--------------------------------------------------------------------------------------------------
+subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar3333, ipc, ip, el)
+ use material, only: &
+   material_phase, &
+   material_homog, &
+   temperature, &
+   temperatureRate, &
+   thermalMapping
+ use lattice, only: &
+   lattice_thermalExpansion33, &
+   lattice_referenceTemperature
+ 
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< grain number
+   ip, &                                                                                            !< integration point number
+   el                                                                                               !< element number
+ real(pReal),   intent(out), dimension(3,3) :: &
+   Li                                                                                               !< thermal velocity gradient
+ real(pReal),   intent(out), dimension(3,3,3,3) :: &
+   dLi_dTstar3333                                                                                   !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero)
+ integer(pInt) :: &
+   phase, &
+   homog, offset
+ real(pReal) :: &
+   T, TRef, TDot  
+   
+ phase = material_phase(ipc,ip,el)
+ homog = material_homog(ip,el)
+ offset = thermalMapping(homog)%p(ip,el)
+ T = temperature(homog)%p(offset)
+ TDot = temperatureRate(homog)%p(offset)
+ TRef = lattice_referenceTemperature(phase)
+ 
+ Li = TDot* &
+      lattice_thermalExpansion33(1:3,1:3,phase)/ &
+      (1.0_pReal + lattice_thermalExpansion33(1:3,1:3,phase)*(T - TRef))
+ dLi_dTstar3333 = 0.0_pReal 
+  
+end subroutine kinematics_thermal_expansion_LiAndItsTangent
+
+end module kinematics_thermal_expansion

code/kinematics_vacancy_strain.f90

@@ -0,0 +1,265 @@
+!--------------------------------------------------------------------------------------------------
+! $Id$
+!--------------------------------------------------------------------------------------------------
+!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
+!> @brief material subroutine incorporating kinematics resulting from vacancy point defects
+!> @details to be done
+!--------------------------------------------------------------------------------------------------
+module kinematics_vacancy_strain
+ use prec, only: &
+   pReal, &
+   pInt
+
+ implicit none
+ private
+ integer(pInt),                       dimension(:),           allocatable,         public, protected :: &
+   kinematics_vacancy_strain_sizePostResults, &                                                                !< cumulative size of post results
+   kinematics_vacancy_strain_offset, &                                                                         !< which kinematics is my current damage mechanism?
+   kinematics_vacancy_strain_instance                                                                          !< instance of damage kinematics mechanism
+
+ integer(pInt),                       dimension(:,:),         allocatable, target, public  :: &
+   kinematics_vacancy_strain_sizePostResult                                                                    !< size of each post result output
+
+ character(len=64),                   dimension(:,:),         allocatable, target, public  :: &
+   kinematics_vacancy_strain_output                                                                            !< name of each post result output
+   
+ integer(pInt),                       dimension(:),           allocatable, target, public  :: &
+   kinematics_vacancy_strain_Noutput                                                                           !< number of outputs per instance of this damage 
+
+ real(pReal),                         dimension(:),           allocatable,         private :: &
+   kinematics_vacancy_strain_coeff
+
+ public :: &
+   kinematics_vacancy_strain_init, &
+   kinematics_vacancy_strain_initialStrain, &
+   kinematics_vacancy_strain_LiAndItsTangent, &
+   kinematics_vacancy_strain_ChemPotAndItsTangent
+
+contains
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief module initialization
+!> @details reads in material parameters, allocates arrays, and does sanity checks
+!--------------------------------------------------------------------------------------------------
+subroutine kinematics_vacancy_strain_init(fileUnit)
+ use, intrinsic :: iso_fortran_env                                                                  ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
+ use debug, only: &
+   debug_level,&
+   debug_constitutive,&
+   debug_levelBasic
+ use IO, only: &
+   IO_read, &
+   IO_lc, &
+   IO_getTag, &
+   IO_isBlank, &
+   IO_stringPos, &
+   IO_stringValue, &
+   IO_floatValue, &
+   IO_intValue, &
+   IO_warning, &
+   IO_error, &
+   IO_timeStamp, &
+   IO_EOF
+ use material, only: &
+   phase_kinematics, &
+   phase_Nkinematics, &
+   phase_Noutput, &
+   KINEMATICS_vacancy_strain_label, &
+   KINEMATICS_vacancy_strain_ID, &
+   material_Nphase, &
+   MATERIAL_partPhase
+ use numerics,only: &
+   worldrank
+
+ implicit none
+ integer(pInt), intent(in) :: fileUnit
+
+ integer(pInt), allocatable, dimension(:) :: chunkPos
+ integer(pInt) :: maxNinstance,phase,instance,kinematics
+ character(len=65536) :: &
+   tag  = '', &
+   line = ''
+
+ mainProcess: if (worldrank == 0) then 
+   write(6,'(/,a)')   ' <<<+-  kinematics_'//KINEMATICS_vacancy_strain_LABEL//' init  -+>>>'
+   write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
+#include "compilation_info.f90"
+ endif mainProcess
+
+ maxNinstance = int(count(phase_kinematics == KINEMATICS_vacancy_strain_ID),pInt)
+ if (maxNinstance == 0_pInt) return
+ 
+ if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
+   write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
+ 
+ allocate(kinematics_vacancy_strain_offset(material_Nphase), source=0_pInt)
+ allocate(kinematics_vacancy_strain_instance(material_Nphase), source=0_pInt)
+ do phase = 1, material_Nphase
+   kinematics_vacancy_strain_instance(phase) = count(phase_kinematics(:,1:phase) == kinematics_vacancy_strain_ID)
+   do kinematics = 1, phase_Nkinematics(phase)
+     if (phase_kinematics(kinematics,phase) == kinematics_vacancy_strain_ID) &
+       kinematics_vacancy_strain_offset(phase) = kinematics
+   enddo    
+ enddo
+   
+ allocate(kinematics_vacancy_strain_sizePostResults(maxNinstance),                     source=0_pInt)
+ allocate(kinematics_vacancy_strain_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
+ allocate(kinematics_vacancy_strain_output(maxval(phase_Noutput),maxNinstance))
+          kinematics_vacancy_strain_output = ''
+ allocate(kinematics_vacancy_strain_Noutput(maxNinstance),                             source=0_pInt) 
+ allocate(kinematics_vacancy_strain_coeff(maxNinstance),                               source=0.0_pReal) 
+
+ rewind(fileUnit)
+ phase = 0_pInt
+ do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase)         ! wind forward to <phase>
+   line = IO_read(fileUnit)
+ enddo
+ 
+ parsingFile: do while (trim(line) /= IO_EOF)                                                       ! read through sections of phase part
+   line = IO_read(fileUnit)
+   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
+   if (IO_getTag(line,'<','>') /= '') then                                                          ! stop at next part
+     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
+     exit                                                                                           
+   endif   
+   if (IO_getTag(line,'[',']') /= '') then                                                          ! next phase section
+     phase = phase + 1_pInt                                                                         ! advance phase section counter
+     cycle                                                                                          ! skip to next line
+   endif
+   if (phase > 0_pInt ) then; if (any(phase_kinematics(:,phase) == KINEMATICS_vacancy_strain_ID)) then         ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
+     instance = kinematics_vacancy_strain_instance(phase)                                                         ! which instance of my damage is present phase
+     chunkPos = IO_stringPos(line)
+     tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt))                                             ! extract key
+     select case(tag)
+       case ('vacancy_strain_coeff')
+         kinematics_vacancy_strain_coeff(instance) = IO_floatValue(line,chunkPos,2_pInt)
+
+     end select
+   endif; endif
+ enddo parsingFile
+
+end subroutine kinematics_vacancy_strain_init
+
+!--------------------------------------------------------------------------------------------------
+!> @brief  report initial vacancy strain based on current vacancy conc deviation from equillibrium
+!--------------------------------------------------------------------------------------------------
+pure function kinematics_vacancy_strain_initialStrain(ipc, ip, el)
+ use math, only: &
+   math_I3
+ use material, only: &
+   material_phase, &
+   material_homog, &
+   vacancyConc, &
+   vacancyfluxMapping
+ use lattice, only: &
+   lattice_equilibriumVacancyConcentration
+ 
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< grain number
+   ip, &                                                                                            !< integration point number
+   el                                                                                               !< element number
+ real(pReal), dimension(3,3) :: &
+   kinematics_vacancy_strain_initialStrain                                                       !< initial thermal strain (should be small strain, though)
+ integer(pInt) :: &
+   phase, &
+   homog, offset, instance
+   
+ phase = material_phase(ipc,ip,el)
+ instance = kinematics_vacancy_strain_instance(phase)
+ homog = material_homog(ip,el)
+ offset = vacancyfluxMapping(homog)%p(ip,el)
+ 
+ kinematics_vacancy_strain_initialStrain = &
+   (vacancyConc(homog)%p(offset) - lattice_equilibriumVacancyConcentration(phase)) * &
+   kinematics_vacancy_strain_coeff(instance)* math_I3
+  
+end function kinematics_vacancy_strain_initialStrain
+
+!--------------------------------------------------------------------------------------------------
+!> @brief  contains the constitutive equation for calculating the velocity gradient  
+!--------------------------------------------------------------------------------------------------
+subroutine kinematics_vacancy_strain_LiAndItsTangent(Li, dLi_dTstar3333, ipc, ip, el)
+ use material, only: &
+   material_phase, &
+   material_homog, &
+   vacancyConc, &
+   vacancyConcRate, &
+   vacancyfluxMapping
+ use math, only: &
+   math_I3
+ use lattice, only: &
+   lattice_equilibriumVacancyConcentration
+ 
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< grain number
+   ip, &                                                                                            !< integration point number
+   el                                                                                               !< element number
+ real(pReal),   intent(out), dimension(3,3) :: &
+   Li                                                                                               !< thermal velocity gradient
+ real(pReal),   intent(out), dimension(3,3,3,3) :: &
+   dLi_dTstar3333                                                                                   !< derivative of Li with respect to Tstar (4th-order tensor)
+ integer(pInt) :: &
+   phase, &
+   instance, &
+   homog, offset
+ real(pReal) :: &
+   Cv, CvEq, CvDot  
+   
+ phase = material_phase(ipc,ip,el)
+ instance = kinematics_vacancy_strain_instance(phase)
+ homog = material_homog(ip,el)
+ offset = vacancyfluxMapping(homog)%p(ip,el)
+ 
+ Cv = vacancyConc(homog)%p(offset)
+ CvDot = vacancyConcRate(homog)%p(offset)
+ CvEq = lattice_equilibriumvacancyConcentration(phase)
+ 
+ Li = CvDot*math_I3* &
+      kinematics_vacancy_strain_coeff(instance)/ &
+      (1.0_pReal + kinematics_vacancy_strain_coeff(instance)*(Cv - CvEq))
+      
+ dLi_dTstar3333 = 0.0_pReal
+  
+end subroutine kinematics_vacancy_strain_LiAndItsTangent
+
+!--------------------------------------------------------------------------------------------------
+!> @brief  contains the kinematic contribution to vacancy chemical potential  
+!--------------------------------------------------------------------------------------------------
+subroutine kinematics_vacancy_strain_ChemPotAndItsTangent(ChemPot, dChemPot_dCv, Tstar_v, Fi0, Fi, ipc, ip, el)
+ use material, only: &
+   material_phase
+ use math, only: &
+   math_inv33, &
+   math_mul33x33, &
+   math_Mandel6to33, &
+   math_transpose33
+ 
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< grain number
+   ip, &                                                                                            !< integration point number
+   el                                                                                               !< element number
+ real(pReal),  intent(in), dimension(6) :: &
+   Tstar_v
+ real(pReal),  intent(in), dimension(3,3) :: &
+   Fi0, Fi
+ real(pReal),   intent(out) :: &
+   ChemPot, dChemPot_dCv                                                               
+ integer(pInt) :: &
+   phase, &
+   instance
+   
+ phase = material_phase(ipc,ip,el)
+ instance = kinematics_vacancy_strain_instance(phase)
+ 
+ ChemPot = -kinematics_vacancy_strain_coeff(instance)* &
+            sum(math_mul33x33(Fi,math_Mandel6to33(Tstar_v))* &
+                math_mul33x33(math_mul33x33(Fi,math_inv33(Fi0)),Fi))
+ dChemPot_dCv = 0.0_pReal
+   
+end subroutine kinematics_vacancy_strain_ChemPotAndItsTangent
+
+end module kinematics_vacancy_strain

code/plastic_isotropic.f90

@@ -0,0 +1,678 @@
+!--------------------------------------------------------------------------------------------------
+! $Id$
+!--------------------------------------------------------------------------------------------------
+!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
+!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
+!> @brief material subroutine for isotropic (ISOTROPIC) plasticity
+!> @details Isotropic (ISOTROPIC) Plasticity which resembles the phenopowerlaw plasticity without
+!! resolving the stress on the slip systems. Will give the response of phenopowerlaw for an
+!! untextured polycrystal
+!--------------------------------------------------------------------------------------------------
+module plastic_isotropic
+#ifdef HDF
+ use hdf5, only: &
+   HID_T
+#endif
+
+ use prec, only: &
+   pReal,&
+   pInt
+ 
+ implicit none
+ private
+ integer(pInt),                       dimension(:),     allocatable,         public, protected :: &
+   plastic_isotropic_sizePostResults                                                                  !< cumulative size of post results
+   
+ integer(pInt),                       dimension(:,:),   allocatable, target, public :: &
+   plastic_isotropic_sizePostResult                                                                   !< size of each post result output
+   
+ character(len=64),                   dimension(:,:),   allocatable, target, public :: &
+   plastic_isotropic_output                                                                           !< name of each post result output
+ 
+ integer(pInt),                       dimension(:),     allocatable, target, public :: &
+   plastic_isotropic_Noutput                                                                          !< number of outputs per instance
+ 
+ enum, bind(c) 
+   enumerator :: undefined_ID, &
+                 flowstress_ID, &
+                 strainrate_ID
+ end enum
+
+ type, private :: tParameters                                                                         !< container type for internal constitutive parameters
+   integer(kind(undefined_ID)), allocatable, dimension(:) :: & 
+     outputID
+  real(pReal) :: &
+     fTaylor, &
+     tau0, &
+     gdot0, &
+     n, &
+     h0, &
+     h0_slopeLnRate, &
+     tausat, &
+     a, &
+     aTolFlowstress, &
+     aTolShear     , &
+     tausat_SinhFitA, &
+     tausat_SinhFitB, &
+     tausat_SinhFitC, &
+     tausat_SinhFitD
+  logical :: &
+     dilatation
+ end type
+
+ type(tParameters), dimension(:), allocatable, private :: param                                       !< containers of constitutive parameters (len Ninstance)
+ 
+ type, private :: tIsotropicState                                                                     !< internal state aliases
+   real(pReal), pointer,     dimension(:) :: &                                                        ! scalars along NipcMyInstance
+     flowstress, &
+     accumulatedShear
+ end type
+ type, private :: tIsotropicAbsTol                                                                  !< internal alias for abs tolerance in state
+   real(pReal), pointer :: &                                                                        ! scalars along NipcMyInstance
+     flowstress, &
+     accumulatedShear
+ end type
+ type(tIsotropicState), allocatable, dimension(:), private :: &                                       !< state aliases per instance
+   state, &
+   state0, &
+   dotState
+ type(tIsotropicAbsTol), allocatable, dimension(:), private :: &                                       !< state aliases per instance
+   stateAbsTol
+
+ public  :: &
+   plastic_isotropic_init, &
+   plastic_isotropic_LpAndItsTangent, &
+   plastic_isotropic_LiAndItsTangent, &
+   plastic_isotropic_dotState, &
+   plastic_isotropic_postResults
+
+contains
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief module initialization
+!> @details reads in material parameters, allocates arrays, and does sanity checks
+!--------------------------------------------------------------------------------------------------
+subroutine plastic_isotropic_init(fileUnit)
+ use, intrinsic :: iso_fortran_env                                                                  ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
+ use debug, only: &
+   debug_level, &
+   debug_constitutive, &
+   debug_levelBasic
+ use numerics, only: &
+   analyticJaco, &
+   worldrank, &
+   numerics_integrator
+ use math, only: &
+   math_Mandel3333to66, &
+   math_Voigt66to3333
+ use IO, only: &
+   IO_read, &
+   IO_lc, &
+   IO_getTag, &
+   IO_isBlank, &
+   IO_stringPos, &
+   IO_stringValue, &
+   IO_floatValue, &
+   IO_error, &
+   IO_timeStamp, &
+   IO_EOF
+ use material, only: &
+   phase_plasticity, &
+   phase_plasticityInstance, &
+   phase_Noutput, &
+   PLASTICITY_ISOTROPIC_label, &
+   PLASTICITY_ISOTROPIC_ID, &
+   material_phase, &
+   plasticState, &
+   MATERIAL_partPhase
+   
+ use lattice  
+
+ implicit none
+ integer(pInt), intent(in) :: fileUnit
+ 
+ 
+ integer(pInt), allocatable, dimension(:) :: chunkPos
+ integer(pInt) :: &
+   o, &
+   phase, & 
+   instance, &
+   maxNinstance, &
+   mySize, &
+   sizeDotState, &
+   sizeState, &
+   sizeDeltaState
+ character(len=65536) :: &
+   tag       = '', &
+   outputtag = '', &
+   line      = '', &
+   extmsg    = ''
+  integer(pInt) :: NipcMyPhase
+
+ mainProcess: if (worldrank == 0) then 
+   write(6,'(/,a)')   ' <<<+-  constitutive_'//PLASTICITY_ISOTROPIC_label//' init  -+>>>'
+   write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
+#include "compilation_info.f90"
+ endif mainProcess
+ 
+ maxNinstance = int(count(phase_plasticity == PLASTICITY_ISOTROPIC_ID),pInt)
+ if (maxNinstance == 0_pInt) return
+
+ if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
+   write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
+
+ allocate(plastic_isotropic_sizePostResults(maxNinstance),                      source=0_pInt)
+ allocate(plastic_isotropic_sizePostResult(maxval(phase_Noutput), maxNinstance),source=0_pInt)
+ allocate(plastic_isotropic_output(maxval(phase_Noutput), maxNinstance))
+          plastic_isotropic_output = ''
+ allocate(plastic_isotropic_Noutput(maxNinstance),                              source=0_pInt)
+
+ allocate(param(maxNinstance))                                                                      ! one container of parameters per instance
+ 
+ rewind(fileUnit)
+ phase = 0_pInt
+ do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partPhase)         ! wind forward to <phase>
+   line = IO_read(fileUnit)
+ enddo
+ 
+ parsingFile: do while (trim(line) /= IO_EOF)                                                       ! read through sections of phase part
+   line = IO_read(fileUnit)
+   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
+   if (IO_getTag(line,'<','>') /= '') then                                                          ! stop at next part
+     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
+     exit                                                                                           
+   endif
+   if (IO_getTag(line,'[',']') /= '') then                                                          ! next section
+     phase = phase + 1_pInt                                                                         ! advance section counter
+     if (phase_plasticity(phase) == PLASTICITY_ISOTROPIC_ID) then
+       instance = phase_plasticityInstance(phase)
+
+     endif
+     cycle                                                                                          ! skip to next line
+   endif
+   if (phase > 0_pInt) then; if (phase_plasticity(phase) == PLASTICITY_ISOTROPIC_ID) then           ! one of my phases. Do not short-circuit here (.and. between if-statements), it's not safe in Fortran
+     instance = phase_plasticityInstance(phase)                                                     ! which instance of my plasticity is present phase
+     allocate(param(instance)%outputID(phase_Noutput(phase)))                                       ! allocate space for IDs of every requested output
+     chunkPos = IO_stringPos(line) 
+     tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt))                                              ! extract key
+     extmsg = trim(tag)//' ('//PLASTICITY_ISOTROPIC_label//')'                                      ! prepare error message identifier
+
+     select case(tag)
+       case ('(output)')
+         outputtag = IO_lc(IO_stringValue(line,chunkPos,2_pInt))
+         select case(outputtag)
+           case ('flowstress')
+             plastic_isotropic_Noutput(instance) = plastic_isotropic_Noutput(instance) + 1_pInt
+             param(instance)%outputID (plastic_isotropic_Noutput(instance)) = flowstress_ID
+             plastic_isotropic_output(plastic_isotropic_Noutput(instance),instance) = outputtag
+           case ('strainrate')
+             plastic_isotropic_Noutput(instance) = plastic_isotropic_Noutput(instance) + 1_pInt
+             param(instance)%outputID (plastic_isotropic_Noutput(instance)) = strainrate_ID
+             plastic_isotropic_output(plastic_isotropic_Noutput(instance),instance) = outputtag
+
+         end select
+
+       case ('/dilatation/')
+         param(instance)%dilatation      = .true.
+
+       case ('tau0')
+         param(instance)%tau0            = IO_floatValue(line,chunkPos,2_pInt)
+         if (param(instance)%tau0 < 0.0_pReal)                  call IO_error(211_pInt,ext_msg=extmsg)
+
+       case ('gdot0')
+         param(instance)%gdot0           = IO_floatValue(line,chunkPos,2_pInt)
+         if (param(instance)%gdot0 <= 0.0_pReal)                call IO_error(211_pInt,ext_msg=extmsg)
+
+       case ('n')
+         param(instance)%n               = IO_floatValue(line,chunkPos,2_pInt)
+         if (param(instance)%n <= 0.0_pReal)                    call IO_error(211_pInt,ext_msg=extmsg)
+
+       case ('h0')
+         param(instance)%h0              = IO_floatValue(line,chunkPos,2_pInt)
+
+       case ('h0_slope','slopelnrate')
+         param(instance)%h0_slopeLnRate  = IO_floatValue(line,chunkPos,2_pInt)
+
+       case ('tausat')
+         param(instance)%tausat          = IO_floatValue(line,chunkPos,2_pInt)
+         if (param(instance)%tausat <= 0.0_pReal)               call IO_error(211_pInt,ext_msg=extmsg)
+
+       case ('tausat_sinhfita')
+         param(instance)%tausat_SinhFitA = IO_floatValue(line,chunkPos,2_pInt)
+
+       case ('tausat_sinhfitb')
+         param(instance)%tausat_SinhFitB = IO_floatValue(line,chunkPos,2_pInt)
+
+       case ('tausat_sinhfitc')
+         param(instance)%tausat_SinhFitC = IO_floatValue(line,chunkPos,2_pInt)
+
+       case ('tausat_sinhfitd')
+         param(instance)%tausat_SinhFitD = IO_floatValue(line,chunkPos,2_pInt)
+
+       case ('a', 'w0')
+         param(instance)%a               = IO_floatValue(line,chunkPos,2_pInt)
+         if (param(instance)%a <= 0.0_pReal)                    call IO_error(211_pInt,ext_msg=extmsg)
+
+       case ('taylorfactor')
+         param(instance)%fTaylor         = IO_floatValue(line,chunkPos,2_pInt)
+         if (param(instance)%fTaylor <= 0.0_pReal)              call IO_error(211_pInt,ext_msg=extmsg)
+
+       case ('atol_flowstress')
+         param(instance)%aTolFlowstress  = IO_floatValue(line,chunkPos,2_pInt)
+         if (param(instance)%aTolFlowstress <= 0.0_pReal)       call IO_error(211_pInt,ext_msg=extmsg)
+
+       case ('atol_shear')
+         param(instance)%aTolShear       = IO_floatValue(line,chunkPos,2_pInt)
+
+       case default
+
+     end select
+   endif; endif
+ enddo parsingFile
+
+ allocate(state(maxNinstance))                                                                      ! internal state aliases
+ allocate(state0(maxNinstance))
+ allocate(dotState(maxNinstance))
+ allocate(stateAbsTol(maxNinstance))
+
+ initializeInstances: do phase = 1_pInt, size(phase_plasticity)                                     ! loop over every plasticity
+   myPhase: if (phase_plasticity(phase) == PLASTICITY_isotropic_ID) then                            ! isolate instances of own constitutive description
+     NipcMyPhase = count(material_phase == phase)                                                   ! number of own material points (including point components ipc)
+     instance = phase_plasticityInstance(phase)
+!--------------------------------------------------------------------------------------------------
+!  sanity checks
+     if (param(instance)%aTolShear <= 0.0_pReal) &
+       param(instance)%aTolShear = 1.0e-6_pReal                                         ! default absolute tolerance 1e-6
+
+!--------------------------------------------------------------------------------------------------
+!  Determine size of postResults array
+     outputsLoop: do o = 1_pInt,plastic_isotropic_Noutput(instance)
+       select case(param(instance)%outputID(o))
+         case(flowstress_ID,strainrate_ID)
+           mySize = 1_pInt
+         case default
+       end select
+  
+       outputFound: if (mySize > 0_pInt) then
+         plastic_isotropic_sizePostResult(o,instance) = mySize
+         plastic_isotropic_sizePostResults(instance) = &
+         plastic_isotropic_sizePostResults(instance) + mySize
+       endif outputFound
+     enddo outputsLoop
+
+!--------------------------------------------------------------------------------------------------
+! allocate state arrays
+     sizeState    = 2_pInt                                                                           ! flowstress, accumulated_shear
+     sizeDotState = sizeState                                                                        ! both evolve
+     sizeDeltaState = 0_pInt                                                                         ! no sudden jumps in state
+     plasticState(phase)%sizeState = sizeState
+     plasticState(phase)%sizeDotState = sizeDotState
+     plasticState(phase)%sizeDeltaState = sizeDeltaState
+     plasticState(phase)%sizePostResults = plastic_isotropic_sizePostResults(instance)
+     plasticState(phase)%nSlip = 1
+     plasticState(phase)%nTwin = 0
+     plasticState(phase)%nTrans= 0
+     allocate(plasticState(phase)%aTolState          (   sizeState))
+
+     allocate(plasticState(phase)%state0             (   sizeState,NipcMyPhase),source=0.0_pReal)
+
+     allocate(plasticState(phase)%partionedState0    (   sizeState,NipcMyPhase),source=0.0_pReal)
+     allocate(plasticState(phase)%subState0          (   sizeState,NipcMyPhase),source=0.0_pReal)
+     allocate(plasticState(phase)%state              (   sizeState,NipcMyPhase),source=0.0_pReal)
+     allocate(plasticState(phase)%dotState           (sizeDotState,NipcMyPhase),source=0.0_pReal)
+     allocate(plasticState(phase)%deltaState       (sizeDeltaState,NipcMyPhase),source=0.0_pReal)
+     if (.not. analyticJaco) then
+       allocate(plasticState(phase)%state_backup     (   sizeState,NipcMyPhase),source=0.0_pReal)
+       allocate(plasticState(phase)%dotState_backup  (sizeDotState,NipcMyPhase),source=0.0_pReal)
+     endif
+     if (any(numerics_integrator == 1_pInt)) then
+       allocate(plasticState(phase)%previousDotState (sizeDotState,NipcMyPhase),source=0.0_pReal)
+       allocate(plasticState(phase)%previousDotState2(sizeDotState,NipcMyPhase),source=0.0_pReal)
+     endif
+     if (any(numerics_integrator == 4_pInt)) &
+       allocate(plasticState(phase)%RK4dotState      (sizeDotState,NipcMyPhase),source=0.0_pReal)
+     if (any(numerics_integrator == 5_pInt)) &
+       allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NipcMyPhase),source=0.0_pReal)
+
+!--------------------------------------------------------------------------------------------------
+! globally required state aliases
+     plasticState(phase)%slipRate           => plasticState(phase)%dotState(2:2,1:NipcMyPhase)
+     plasticState(phase)%accumulatedSlip    => plasticState(phase)%state   (2:2,1:NipcMyPhase)
+
+!--------------------------------------------------------------------------------------------------
+! locally defined state aliases
+     state(instance)%flowstress             => plasticState(phase)%state    (1,1:NipcMyPhase)
+     state0(instance)%flowstress            => plasticState(phase)%state0   (1,1:NipcMyPhase)
+     dotState(instance)%flowstress          => plasticState(phase)%dotState (1,1:NipcMyPhase)
+     stateAbsTol(instance)%flowstress       => plasticState(phase)%aTolState(1)
+
+     state(instance)%accumulatedShear       => plasticState(phase)%state    (2,1:NipcMyPhase)
+     state0(instance)%accumulatedShear      => plasticState(phase)%state0   (2,1:NipcMyPhase)
+     dotState(instance)%accumulatedShear    => plasticState(phase)%dotState (2,1:NipcMyPhase)
+     stateAbsTol(instance)%accumulatedShear => plasticState(phase)%aTolState(2)
+
+!--------------------------------------------------------------------------------------------------
+! init state
+     state0(instance)%flowstress       = param(instance)%tau0
+     state0(instance)%accumulatedShear = 0.0_pReal
+
+!--------------------------------------------------------------------------------------------------
+! init absolute state tolerances
+     stateAbsTol(instance)%flowstress       = param(instance)%aTolFlowstress
+     stateAbsTol(instance)%accumulatedShear = param(instance)%aTolShear
+
+   endif myPhase
+ enddo initializeInstances
+
+end subroutine plastic_isotropic_init
+
+!--------------------------------------------------------------------------------------------------
+!> @brief calculates plastic velocity gradient and its tangent
+!--------------------------------------------------------------------------------------------------
+subroutine plastic_isotropic_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
+ use debug, only: &
+   debug_level, &
+   debug_constitutive, &
+   debug_levelBasic, &
+   debug_levelExtensive, &
+   debug_levelSelective, &
+   debug_e, &
+   debug_i, &
+   debug_g
+ use math, only: &
+   math_mul6x6, &
+   math_Mandel6to33, &
+   math_Plain3333to99, &
+   math_deviatoric33, &
+   math_mul33xx33, &
+   math_transpose33
+ use material, only: &
+   phaseAt, phasememberAt, &
+   plasticState, &
+   material_phase, &
+   phase_plasticityInstance
+
+ implicit none
+ real(pReal), dimension(3,3), intent(out) :: &
+   Lp                                                                                               !< plastic velocity gradient
+ real(pReal), dimension(9,9), intent(out) :: &
+   dLp_dTstar99                                                                                     !< derivative of Lp with respect to 2nd Piola Kirchhoff stress
+
+ real(pReal), dimension(6),   intent(in) :: &
+   Tstar_v                                                                                          !< 2nd Piola Kirchhoff stress tensor in Mandel notation
+ integer(pInt),               intent(in) :: &
+   ipc, &                                                                                           !< component-ID of integration point
+   ip, &                                                                                            !< integration point
+   el                                                                                               !< element
+
+ real(pReal), dimension(3,3) :: &
+   Tstar_dev_33                                                                                     !< deviatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor
+ real(pReal), dimension(3,3,3,3) :: &
+   dLp_dTstar_3333                                                                                  !< derivative of Lp with respect to Tstar as 4th order tensor
+ real(pReal) :: &
+   gamma_dot, &                                                                                     !< strainrate
+   norm_Tstar_dev, &                                                                                !< euclidean norm of Tstar_dev
+   squarenorm_Tstar_dev                                                                             !< square of the euclidean norm of Tstar_dev
+ integer(pInt) :: &
+   instance, of, &
+   k, l, m, n
+
+ of = phasememberAt(ipc,ip,el)                                                                      ! phasememberAt should be tackled by material and be renamed to material_phasemember
+ instance = phase_plasticityInstance(phaseAt(ipc,ip,el))                                            ! "phaseAt" equivalent to "material_phase" !!
+
+ Tstar_dev_33 = math_deviatoric33(math_Mandel6to33(Tstar_v))                                        ! deviatoric part of 2nd Piola-Kirchhoff stress
+ squarenorm_Tstar_dev = math_mul33xx33(Tstar_dev_33,Tstar_dev_33)
+ norm_Tstar_dev = sqrt(squarenorm_Tstar_dev) 
+
+ if (norm_Tstar_dev <= 0.0_pReal) then                                                              ! Tstar == 0 --> both Lp and dLp_dTstar are zero
+   Lp = 0.0_pReal
+   dLp_dTstar99 = 0.0_pReal
+ else
+   gamma_dot = param(instance)%gdot0 &
+             * ( sqrt(1.5_pReal) * norm_Tstar_dev / param(instance)%fTaylor / state(instance)%flowstress(of) ) &
+             **param(instance)%n
+
+   Lp = Tstar_dev_33/norm_Tstar_dev * gamma_dot/param(instance)%fTaylor
+
+   if (iand(debug_level(debug_constitutive), debug_levelExtensive) /= 0_pInt &
+       .and. ((el == debug_e .and. ip == debug_i .and. ipc == debug_g) &
+              .or. .not. iand(debug_level(debug_constitutive),debug_levelSelective) /= 0_pInt)) then
+     write(6,'(a,i8,1x,i2,1x,i3)') '<< CONST isotropic >> at el ip g ',el,ip,ipc
+     write(6,'(/,a,/,3(12x,3(f12.4,1x)/))') '<< CONST isotropic >> Tstar (dev) / MPa', &
+                                      math_transpose33(Tstar_dev_33(1:3,1:3))*1.0e-6_pReal
+     write(6,'(/,a,/,f12.5)') '<< CONST isotropic >> norm Tstar / MPa', norm_Tstar_dev*1.0e-6_pReal
+     write(6,'(/,a,/,f12.5)') '<< CONST isotropic >> gdot', gamma_dot
+   end if
+!--------------------------------------------------------------------------------------------------
+! Calculation of the tangent of Lp
+   forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
+     dLp_dTstar_3333(k,l,m,n) = (param(instance)%n-1.0_pReal) * &
+                                      Tstar_dev_33(k,l)*Tstar_dev_33(m,n) / squarenorm_Tstar_dev
+   forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt) &
+     dLp_dTstar_3333(k,l,k,l) = dLp_dTstar_3333(k,l,k,l) + 1.0_pReal
+   forall (k=1_pInt:3_pInt,m=1_pInt:3_pInt) &
+     dLp_dTstar_3333(k,k,m,m) = dLp_dTstar_3333(k,k,m,m) - 1.0_pReal/3.0_pReal
+   dLp_dTstar99 = math_Plain3333to99(gamma_dot / param(instance)%fTaylor * &
+                                      dLp_dTstar_3333 / norm_Tstar_dev)
+ end if
+end subroutine plastic_isotropic_LpAndItsTangent
+
+!--------------------------------------------------------------------------------------------------
+!> @brief calculates plastic velocity gradient and its tangent
+!--------------------------------------------------------------------------------------------------
+subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dTstar_3333,Tstar_v,ipc,ip,el)
+ use math, only: &
+   math_mul6x6, &
+   math_Mandel6to33, &
+   math_Plain3333to99, &
+   math_spherical33, &
+   math_mul33xx33
+ use material, only: &
+   phaseAt, phasememberAt, &
+   plasticState, &
+   material_phase, &
+   phase_plasticityInstance
+
+ implicit none
+ real(pReal), dimension(3,3), intent(out) :: &
+   Li                                                                                               !< plastic velocity gradient
+
+ real(pReal), dimension(6),   intent(in) :: &
+   Tstar_v                                                                                          !< 2nd Piola Kirchhoff stress tensor in Mandel notation
+ integer(pInt),               intent(in) :: &
+   ipc, &                                                                                           !< component-ID of integration point
+   ip, &                                                                                            !< integration point
+   el                                                                                               !< element
+
+ real(pReal), dimension(3,3) :: &
+   Tstar_sph_33                                                                                     !< sphiatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor
+ real(pReal), dimension(3,3,3,3), intent(out)  :: &
+   dLi_dTstar_3333                                                                                  !< derivative of Li with respect to Tstar as 4th order tensor
+ real(pReal) :: &
+   gamma_dot, &                                                                                     !< strainrate
+   norm_Tstar_sph, &                                                                                !< euclidean norm of Tstar_sph
+   squarenorm_Tstar_sph                                                                             !< square of the euclidean norm of Tstar_sph
+ integer(pInt) :: &
+   instance, of, &
+   k, l, m, n
+
+ of = phasememberAt(ipc,ip,el)                                                                      ! phasememberAt should be tackled by material and be renamed to material_phasemember
+ instance = phase_plasticityInstance(phaseAt(ipc,ip,el))                                            ! "phaseAt" equivalent to "material_phase" !!
+
+ Tstar_sph_33 = math_spherical33(math_Mandel6to33(Tstar_v))                                         ! spherical part of 2nd Piola-Kirchhoff stress
+ squarenorm_Tstar_sph = math_mul33xx33(Tstar_sph_33,Tstar_sph_33)
+ norm_Tstar_sph = sqrt(squarenorm_Tstar_sph) 
+
+ if (param(instance)%dilatation) then
+     if (norm_Tstar_sph <= 0.0_pReal) then                                                          ! Tstar == 0 --> both Li and dLi_dTstar are zero
+       Li = 0.0_pReal
+       dLi_dTstar_3333 = 0.0_pReal
+     else
+       gamma_dot = param(instance)%gdot0 &
+                   * (sqrt(1.5_pReal) * norm_Tstar_sph / param(instance)%fTaylor / state(instance)%flowstress(of) ) &
+                   **param(instance)%n
+
+       Li = Tstar_sph_33/norm_Tstar_sph * gamma_dot/param(instance)%fTaylor
+
+       !--------------------------------------------------------------------------------------------------
+       ! Calculation of the tangent of Li
+       forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
+         dLi_dTstar_3333(k,l,m,n) = (param(instance)%n-1.0_pReal) * &
+                                          Tstar_sph_33(k,l)*Tstar_sph_33(m,n) / squarenorm_Tstar_sph
+       forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt) &
+         dLi_dTstar_3333(k,l,k,l) = dLi_dTstar_3333(k,l,k,l) + 1.0_pReal
+
+       dLi_dTstar_3333 = gamma_dot / param(instance)%fTaylor * &
+                                          dLi_dTstar_3333 / norm_Tstar_sph
+     endif
+ endif
+ 
+end subroutine plastic_isotropic_LiAndItsTangent
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief calculates the rate of change of microstructure
+!--------------------------------------------------------------------------------------------------
+subroutine plastic_isotropic_dotState(Tstar_v,ipc,ip,el)
+ use math, only: &
+   math_mul6x6
+ use material, only: &
+   phaseAt, phasememberAt, &
+   plasticState, &
+   material_phase, &
+   phase_plasticityInstance
+ 
+ implicit none
+ real(pReal), dimension(6), intent(in):: &
+   Tstar_v                                                                                          !< 2nd Piola Kirchhoff stress tensor in Mandel notation
+ integer(pInt),             intent(in) :: &
+   ipc, &                                                                                           !< component-ID of integration point
+   ip, &                                                                                            !< integration point
+   el                                                                                               !< element
+ real(pReal), dimension(6) :: &
+   Tstar_dev_v                                                                                      !< deviatoric 2nd Piola Kirchhoff stress tensor in Mandel notation
+ real(pReal) :: &
+   gamma_dot, &                                                                                     !< strainrate
+   hardening, &                                                                                     !< hardening coefficient
+   saturation, &                                                                                    !< saturation flowstress
+   norm_Tstar_v                                                                                     !< euclidean norm of Tstar_dev
+ integer(pInt) :: &
+   instance, &                                                                                      !< instance of my instance (unique number of my constitutive model)
+   of                                                                                               !< shortcut notation for offset position in state array
+
+ of = phasememberAt(ipc,ip,el)                                                                      ! phasememberAt should be tackled by material and be renamed to material_phasemember
+ instance = phase_plasticityInstance(phaseAt(ipc,ip,el))                                            ! "phaseAt" equivalent to "material_phase" !!
+
+!--------------------------------------------------------------------------------------------------
+! norm of (deviatoric) 2nd Piola-Kirchhoff stress
+ if (param(instance)%dilatation) then
+   norm_Tstar_v = sqrt(math_mul6x6(Tstar_v,Tstar_v))
+ else
+   Tstar_dev_v(1:3) = Tstar_v(1:3) - sum(Tstar_v(1:3))/3.0_pReal
+   Tstar_dev_v(4:6) = Tstar_v(4:6)
+   norm_Tstar_v = sqrt(math_mul6x6(Tstar_dev_v,Tstar_dev_v))
+ end if
+!--------------------------------------------------------------------------------------------------
+! strain rate 
+ gamma_dot = param(instance)%gdot0 * ( sqrt(1.5_pReal) * norm_Tstar_v & 
+            / &!-----------------------------------------------------------------------------------
+           (param(instance)%fTaylor*state(instance)%flowstress(of) ))**param(instance)%n
+ 
+!--------------------------------------------------------------------------------------------------
+! hardening coefficient
+ if (abs(gamma_dot) > 1e-12_pReal) then
+   if (abs(param(instance)%tausat_SinhFitA) <= tiny(0.0_pReal)) then
+     saturation = param(instance)%tausat
+   else
+     saturation = (  param(instance)%tausat &
+                   + ( log(  ( gamma_dot / param(instance)%tausat_SinhFitA&
+                               )**(1.0_pReal / param(instance)%tausat_SinhFitD)&
+                            + sqrt(  ( gamma_dot / param(instance)%tausat_SinhFitA &
+                                      )**(2.0_pReal / param(instance)%tausat_SinhFitD) &
+                                   + 1.0_pReal ) &
+                            ) & ! asinh(K) = ln(K + sqrt(K^2 +1))
+                       )**(1.0_pReal / param(instance)%tausat_SinhFitC) &
+                   / (  param(instance)%tausat_SinhFitB &
+                      * (gamma_dot / param(instance)%gdot0)**(1.0_pReal / param(instance)%n) &
+                      ) &
+                   )
+   endif
+   hardening = ( param(instance)%h0 + param(instance)%h0_slopeLnRate * log(gamma_dot) ) &
+               * abs( 1.0_pReal - state(instance)%flowstress(of)/saturation )**param(instance)%a &
+               * sign(1.0_pReal, 1.0_pReal - state(instance)%flowstress(of)/saturation)
+ else
+   hardening = 0.0_pReal
+ endif
+
+ dotState(instance)%flowstress      (of) = hardening * gamma_dot
+ dotState(instance)%accumulatedShear(of) =             gamma_dot
+
+end subroutine plastic_isotropic_dotState
+
+!--------------------------------------------------------------------------------------------------
+!> @brief return array of constitutive results
+!--------------------------------------------------------------------------------------------------
+function plastic_isotropic_postResults(Tstar_v,ipc,ip,el)
+ use math, only: &
+   math_mul6x6
+ use material, only: &
+   material_phase, &
+   plasticState, &
+   phaseAt, phasememberAt, &
+   phase_plasticityInstance
+
+ implicit none
+ real(pReal), dimension(6),  intent(in) :: &
+   Tstar_v                                                                                          !< 2nd Piola Kirchhoff stress tensor in Mandel notation
+ integer(pInt),              intent(in) :: &
+   ipc, &                                                                                           !< component-ID of integration point
+   ip, &                                                                                            !< integration point
+   el                                                                                               !< element
+ real(pReal), dimension(plastic_isotropic_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
+                                           plastic_isotropic_postResults
+
+ real(pReal), dimension(6) :: &
+   Tstar_dev_v                                                                                      !< deviatoric 2nd Piola Kirchhoff stress tensor in Mandel notation
+ real(pReal) :: &
+   norm_Tstar_v                                                                                     ! euclidean norm of Tstar_dev
+ integer(pInt) :: &
+   instance, &                                                                                      !< instance of my instance (unique number of my constitutive model)
+   of, &                                                                                            !< shortcut notation for offset position in state array
+   c, &
+   o
+
+ of = phasememberAt(ipc,ip,el)                                                                      ! phasememberAt should be tackled by material and be renamed to material_phasemember
+ instance = phase_plasticityInstance(phaseAt(ipc,ip,el))                                            ! "phaseAt" equivalent to "material_phase" !!
+ 
+!--------------------------------------------------------------------------------------------------
+! norm of (deviatoric) 2nd Piola-Kirchhoff stress
+ if (param(instance)%dilatation) then
+   norm_Tstar_v = sqrt(math_mul6x6(Tstar_v,Tstar_v))
+ else
+   Tstar_dev_v(1:3) = Tstar_v(1:3) - sum(Tstar_v(1:3))/3.0_pReal
+   Tstar_dev_v(4:6) = Tstar_v(4:6)
+   norm_Tstar_v = sqrt(math_mul6x6(Tstar_dev_v,Tstar_dev_v))
+ end if
+ 
+ c = 0_pInt
+ plastic_isotropic_postResults = 0.0_pReal
+
+ outputsLoop: do o = 1_pInt,plastic_isotropic_Noutput(instance)
+   select case(param(instance)%outputID(o))
+     case (flowstress_ID)
+       plastic_isotropic_postResults(c+1_pInt) = state(instance)%flowstress(of)
+       c = c + 1_pInt
+     case (strainrate_ID)
+       plastic_isotropic_postResults(c+1_pInt) = &
+                param(instance)%gdot0 * (            sqrt(1.5_pReal) * norm_Tstar_v & 
+             / &!----------------------------------------------------------------------------------
+              (param(instance)%fTaylor * state(instance)%flowstress(of)) ) ** param(instance)%n
+       c = c + 1_pInt
+   end select
+ enddo outputsLoop
+
+end function plastic_isotropic_postResults
+
+
+end module plastic_isotropic

code/plastic_j2.f90

@@ -0,0 +1,579 @@
+!--------------------------------------------------------------------------------------------------
+! $Id$
+!--------------------------------------------------------------------------------------------------
+!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
+!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
+!> @brief material subroutine for isotropic (J2) plasticity
+!> @details Isotropic (J2) Plasticity which resembles the phenopowerlaw plasticity without
+!! resolving the stress on the slip systems. Will give the response of phenopowerlaw for an
+!! untextured polycrystal
+!--------------------------------------------------------------------------------------------------
+module plastic_j2
+#ifdef HDF
+ use hdf5, only: &
+   HID_T
+#endif
+
+ use prec, only: &
+   pReal,&
+   pInt
+ 
+ implicit none
+ private
+ integer(pInt),                       dimension(:),     allocatable,         public, protected :: &
+   plastic_j2_sizePostResults                                                                  !< cumulative size of post results
+   
+ integer(pInt),                       dimension(:,:),   allocatable, target, public :: &
+   plastic_j2_sizePostResult                                                                   !< size of each post result output
+   
+ character(len=64),                   dimension(:,:),   allocatable, target, public :: &
+   plastic_j2_output                                                                           !< name of each post result output
+ 
+ integer(pInt),                       dimension(:),     allocatable, target, public :: &
+   plastic_j2_Noutput                                                                          !< number of outputs per instance
+ real(pReal),                         dimension(:),     allocatable,         private :: &
+   plastic_j2_fTaylor, &                                                                       !< Taylor factor
+   plastic_j2_tau0, &                                                                          !< initial plastic stress
+   plastic_j2_gdot0, &                                                                         !< reference velocity
+   plastic_j2_n, &                                                                             !< Visco-plastic parameter
+!--------------------------------------------------------------------------------------------------
+! h0 as function of h0 = A + B log (gammadot) 
+   plastic_j2_h0, &
+   plastic_j2_h0_slopeLnRate, &
+   plastic_j2_tausat, &                                                                        !< final plastic stress
+   plastic_j2_a, &
+   plastic_j2_aTolResistance, &
+   plastic_j2_aTolShear, &
+!--------------------------------------------------------------------------------------------------
+! tausat += (asinh((gammadot / SinhFitA)**(1 / SinhFitD)))**(1 / SinhFitC) / (SinhFitB * (gammadot / gammadot0)**(1/n))
+   plastic_j2_tausat_SinhFitA, &                                                               !< fitting parameter for normalized strain rate vs. stress function
+   plastic_j2_tausat_SinhFitB, &                                                               !< fitting parameter for normalized strain rate vs. stress function
+   plastic_j2_tausat_SinhFitC, &                                                               !< fitting parameter for normalized strain rate vs. stress function
+   plastic_j2_tausat_SinhFitD                                                                  !< fitting parameter for normalized strain rate vs. stress function
+
+ enum, bind(c) 
+   enumerator :: undefined_ID, &
+                 flowstress_ID, &
+                 strainrate_ID
+ end enum
+ integer(kind(undefined_ID)),           dimension(:,:),   allocatable,        private :: & 
+   plastic_j2_outputID                                                                         !< ID of each post result output
+ 
+
+#ifdef HDF 
+ type plastic_j2_tOutput
+   real(pReal),                         dimension(:),     allocatable,         private :: &
+     flowstress, &
+     strainrate
+ logical :: flowstressActive = .false., strainrateActive = .false.                ! if we can write the output block wise, this is not needed anymore because we can do an if(allocated(xxx))                                 
+ end type plastic_j2_tOutput
+ type(plastic_j2_tOutput), allocatable, dimension(:) :: plastic_j2_Output2
+integer(HID_T), allocatable, dimension(:) :: outID
+#endif
+
+
+ public  :: &
+   plastic_j2_init, &
+   plastic_j2_LpAndItsTangent, &
+   plastic_j2_dotState, &
+   plastic_j2_postResults
+
+contains
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief module initialization
+!> @details reads in material parameters, allocates arrays, and does sanity checks
+!--------------------------------------------------------------------------------------------------
+subroutine plastic_j2_init(fileUnit)
+ use, intrinsic :: iso_fortran_env                                                                  ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
+#ifdef HDF 
+ use hdf5
+#endif 
+ use debug, only: &
+   debug_level, &
+   debug_constitutive, &
+   debug_levelBasic
+ use numerics, only: &
+   analyticJaco, &
+   worldrank, &
+   numerics_integrator
+ use math, only: &
+   math_Mandel3333to66, &
+   math_Voigt66to3333
+ use IO, only: &
+   IO_read, &
+   IO_lc, &
+   IO_getTag, &
+   IO_isBlank, &
+   IO_stringPos, &
+   IO_stringValue, &
+   IO_floatValue, &
+   IO_error, &
+   IO_timeStamp, &
+#ifdef HDF 
+   tempResults, &
+   HDF5_addGroup, &
+   HDF5_addScalarDataset,&
+#endif
+   IO_EOF
+ use material, only: &
+   phase_plasticity, &
+   phase_plasticityInstance, &
+   phase_Noutput, &
+   PLASTICITY_J2_label, &
+   PLASTICITY_J2_ID, &
+   material_phase, &
+   plasticState, &
+   MATERIAL_partPhase
+   
+ use lattice  
+
+ implicit none
+ integer(pInt), intent(in) :: fileUnit
+ 
+ 
+ integer(pInt), allocatable, dimension(:) :: chunkPos
+ integer(pInt) :: &
+   o, &
+   phase, & 
+   maxNinstance, &
+   instance, &
+   mySize, &
+   sizeDotState, &
+   sizeState, &
+   sizeDeltaState
+ character(len=65536) :: &
+   tag  = '', &
+   line = ''
+  integer(pInt) :: NofMyPhase
+
+#ifdef HDF 
+ character(len=5) :: &
+   str1
+ integer(HID_T) :: ID,ID2,ID4
+#endif
+
+ mainProcess: if (worldrank == 0) then 
+   write(6,'(/,a)')   ' <<<+-  constitutive_'//PLASTICITY_J2_label//' init  -+>>>'
+   write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
+#include "compilation_info.f90"
+ endif mainProcess
+ 
+ maxNinstance = int(count(phase_plasticity == PLASTICITY_J2_ID),pInt)
+ if (maxNinstance == 0_pInt) return
+
+ if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
+   write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
+
+#ifdef HDF 
+  allocate(plastic_j2_Output2(maxNinstance))
+  allocate(outID(maxNinstance))
+#endif
+
+ allocate(plastic_j2_sizePostResults(maxNinstance),                      source=0_pInt)
+ allocate(plastic_j2_sizePostResult(maxval(phase_Noutput), maxNinstance),source=0_pInt)
+ allocate(plastic_j2_output(maxval(phase_Noutput), maxNinstance))
+          plastic_j2_output = ''
+ allocate(plastic_j2_outputID(maxval(phase_Noutput),maxNinstance),       source=undefined_ID)
+ allocate(plastic_j2_Noutput(maxNinstance),                              source=0_pInt)
+ allocate(plastic_j2_fTaylor(maxNinstance),                              source=0.0_pReal)
+ allocate(plastic_j2_tau0(maxNinstance),                                 source=0.0_pReal)
+ allocate(plastic_j2_gdot0(maxNinstance),                                source=0.0_pReal)
+ allocate(plastic_j2_n(maxNinstance),                                    source=0.0_pReal)
+ allocate(plastic_j2_h0(maxNinstance),                                   source=0.0_pReal)
+ allocate(plastic_j2_h0_slopeLnRate(maxNinstance),                       source=0.0_pReal)
+ allocate(plastic_j2_tausat(maxNinstance),                               source=0.0_pReal)
+ allocate(plastic_j2_a(maxNinstance),                                    source=0.0_pReal)
+ allocate(plastic_j2_aTolResistance(maxNinstance),                       source=0.0_pReal)
+ allocate(plastic_j2_aTolShear     (maxNinstance),                       source=0.0_pReal)
+ allocate(plastic_j2_tausat_SinhFitA(maxNinstance),                      source=0.0_pReal)
+ allocate(plastic_j2_tausat_SinhFitB(maxNinstance),                      source=0.0_pReal)
+ allocate(plastic_j2_tausat_SinhFitC(maxNinstance),                      source=0.0_pReal)
+ allocate(plastic_j2_tausat_SinhFitD(maxNinstance),                      source=0.0_pReal)
+ 
+ rewind(fileUnit)
+ phase = 0_pInt
+ do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= material_partPhase)         ! wind forward to <phase>
+   line = IO_read(fileUnit)
+ enddo
+ 
+ parsingFile: do while (trim(line) /= IO_EOF)                                                       ! read through sections of phase part
+   line = IO_read(fileUnit)
+   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
+   if (IO_getTag(line,'<','>') /= '') then                                                          ! stop at next part
+     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
+     exit                                                                                           
+   endif
+   if (IO_getTag(line,'[',']') /= '') then                                                          ! next section
+     phase = phase + 1_pInt                                                                         ! advance section counter
+     if (phase_plasticity(phase) == PLASTICITY_J2_ID) then
+       instance = phase_plasticityInstance(phase)
+#ifdef HDF
+       outID(instance)=HDF5_addGroup(str1,tempResults)
+#endif
+     endif
+     cycle                                                                                          ! skip to next line
+   endif
+   if (phase > 0_pInt ) then; if (phase_plasticity(phase) == PLASTICITY_J2_ID) then                 ! one of my phases. Do not short-circuit here (.and. between if-statements), it's not safe in Fortran
+     instance = phase_plasticityInstance(phase)                                                     ! which instance of my plasticity is present phase
+     chunkPos = IO_stringPos(line) 
+     tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt))                                             ! extract key
+
+     select case(tag)
+       case ('(output)')
+         select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt)))
+           case ('flowstress')
+             plastic_j2_Noutput(instance) = plastic_j2_Noutput(instance) + 1_pInt
+             plastic_j2_outputID(plastic_j2_Noutput(instance),instance) = flowstress_ID
+             plastic_j2_output(plastic_j2_Noutput(instance),instance) = &
+                                                IO_lc(IO_stringValue(line,chunkPos,2_pInt))
+#ifdef HDF 
+             call HDF5_addScalarDataset(outID(instance),myConstituents,'flowstress','MPa')
+             allocate(plastic_j2_Output2(instance)%flowstress(myConstituents))
+             plastic_j2_Output2(instance)%flowstressActive = .true.
+#endif
+           case ('strainrate')
+             plastic_j2_Noutput(instance) = plastic_j2_Noutput(instance) + 1_pInt
+             plastic_j2_outputID(plastic_j2_Noutput(instance),instance) = strainrate_ID
+             plastic_j2_output(plastic_j2_Noutput(instance),instance) = &
+                                                IO_lc(IO_stringValue(line,chunkPos,2_pInt))
+#ifdef HDF 
+             call HDF5_addScalarDataset(outID(instance),myConstituents,'strainrate','1/s')
+             allocate(plastic_j2_Output2(instance)%strainrate(myConstituents))
+             plastic_j2_Output2(instance)%strainrateActive = .true.
+#endif
+           case default
+
+         end select
+       case ('tau0')
+         plastic_j2_tau0(instance)         = IO_floatValue(line,chunkPos,2_pInt)
+         if (plastic_j2_tau0(instance) < 0.0_pReal) &
+           call IO_error(211_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_J2_label//')')
+       case ('gdot0')
+         plastic_j2_gdot0(instance)        = IO_floatValue(line,chunkPos,2_pInt)
+         if (plastic_j2_gdot0(instance) <= 0.0_pReal) &
+           call IO_error(211_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_J2_label//')')
+       case ('n')
+         plastic_j2_n(instance)            = IO_floatValue(line,chunkPos,2_pInt)
+         if (plastic_j2_n(instance) <= 0.0_pReal) &
+           call IO_error(211_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_J2_label//')')
+       case ('h0')
+         plastic_j2_h0(instance)           = IO_floatValue(line,chunkPos,2_pInt)
+       case ('h0_slope','slopelnrate')
+         plastic_j2_h0_slopeLnRate(instance)  = IO_floatValue(line,chunkPos,2_pInt)
+       case ('tausat')
+         plastic_j2_tausat(instance)          = IO_floatValue(line,chunkPos,2_pInt)
+         if (plastic_j2_tausat(instance) <= 0.0_pReal) &
+           call IO_error(211_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_J2_label//')')
+       case ('tausat_sinhfita')
+         plastic_j2_tausat_SinhFitA(instance) = IO_floatValue(line,chunkPos,2_pInt)
+       case ('tausat_sinhfitb')
+         plastic_j2_tausat_SinhFitB(instance) = IO_floatValue(line,chunkPos,2_pInt)
+       case ('tausat_sinhfitc')
+         plastic_j2_tausat_SinhFitC(instance) = IO_floatValue(line,chunkPos,2_pInt)
+       case ('tausat_sinhfitd')
+         plastic_j2_tausat_SinhFitD(instance) = IO_floatValue(line,chunkPos,2_pInt)
+       case ('a', 'w0')
+         plastic_j2_a(instance)               = IO_floatValue(line,chunkPos,2_pInt)
+         if (plastic_j2_a(instance) <= 0.0_pReal) &
+           call IO_error(211_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_J2_label//')')
+       case ('taylorfactor')
+         plastic_j2_fTaylor(instance)         = IO_floatValue(line,chunkPos,2_pInt)
+         if (plastic_j2_fTaylor(instance) <= 0.0_pReal) &
+           call IO_error(211_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_J2_label//')')
+       case ('atol_resistance')
+         plastic_j2_aTolResistance(instance)  = IO_floatValue(line,chunkPos,2_pInt)
+         if (plastic_j2_aTolResistance(instance) <= 0.0_pReal) &
+           call IO_error(211_pInt,ext_msg=trim(tag)//' ('//PLASTICITY_J2_label//')')
+       case ('atol_shear')
+         plastic_j2_aTolShear(instance)  = IO_floatValue(line,chunkPos,2_pInt)
+
+       case default
+
+     end select
+   endif; endif
+ enddo parsingFile
+
+ initializeInstances: do phase = 1_pInt, size(phase_plasticity)
+   myPhase: if (phase_plasticity(phase) == PLASTICITY_j2_ID) then
+     NofMyPhase=count(material_phase==phase)
+     instance = phase_plasticityInstance(phase)
+!--------------------------------------------------------------------------------------------------
+!  sanity checks
+     if (plastic_j2_aTolShear(instance) <= 0.0_pReal) &
+       plastic_j2_aTolShear(instance) = 1.0e-6_pReal                                         ! default absolute tolerance 1e-6
+
+!--------------------------------------------------------------------------------------------------
+!  Determine size of postResults array
+     outputsLoop: do o = 1_pInt,plastic_j2_Noutput(instance)
+       select case(plastic_j2_outputID(o,instance))
+         case(flowstress_ID,strainrate_ID)
+           mySize = 1_pInt
+         case default
+       end select
+  
+       outputFound: if (mySize > 0_pInt) then
+         plastic_j2_sizePostResult(o,instance) = mySize
+         plastic_j2_sizePostResults(instance) = &
+         plastic_j2_sizePostResults(instance) + mySize
+       endif outputFound
+     enddo outputsLoop
+
+!--------------------------------------------------------------------------------------------------
+! allocate state arrays
+     sizeState    = 2_pInt
+     sizeDotState = sizeState
+     sizeDeltaState = 0_pInt
+     plasticState(phase)%sizeState = sizeState
+     plasticState(phase)%sizeDotState = sizeDotState
+     plasticState(phase)%sizeDeltaState = sizeDeltaState
+     plasticState(phase)%sizePostResults = plastic_j2_sizePostResults(instance)
+     plasticState(phase)%nSlip = 1
+     plasticState(phase)%nTwin = 0
+     plasticState(phase)%nTrans= 0
+     allocate(plasticState(phase)%aTolState          (   sizeState))
+     plasticState(phase)%aTolState(1) = plastic_j2_aTolResistance(instance)
+     plasticState(phase)%aTolState(2) = plastic_j2_aTolShear(instance)
+     allocate(plasticState(phase)%state0             (   sizeState,NofMyPhase))
+     plasticState(phase)%state0(1,1:NofMyPhase) = plastic_j2_tau0(instance)
+     plasticState(phase)%state0(2,1:NofMyPhase) = 0.0_pReal
+     allocate(plasticState(phase)%partionedState0    (   sizeState,NofMyPhase),source=0.0_pReal)
+     allocate(plasticState(phase)%subState0          (   sizeState,NofMyPhase),source=0.0_pReal)
+     allocate(plasticState(phase)%state              (   sizeState,NofMyPhase),source=0.0_pReal)
+     allocate(plasticState(phase)%dotState           (sizeDotState,NofMyPhase),source=0.0_pReal)
+     allocate(plasticState(phase)%deltaState       (sizeDeltaState,NofMyPhase),source=0.0_pReal)
+     if (.not. analyticJaco) then
+       allocate(plasticState(phase)%state_backup     (   sizeState,NofMyPhase),source=0.0_pReal)
+       allocate(plasticState(phase)%dotState_backup  (sizeDotState,NofMyPhase),source=0.0_pReal)
+     endif
+     if (any(numerics_integrator == 1_pInt)) then
+       allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase),source=0.0_pReal)
+       allocate(plasticState(phase)%previousDotState2(sizeDotState,NofMyPhase),source=0.0_pReal)
+     endif
+     if (any(numerics_integrator == 4_pInt)) &
+       allocate(plasticState(phase)%RK4dotState      (sizeDotState,NofMyPhase),source=0.0_pReal)
+     if (any(numerics_integrator == 5_pInt)) &
+       allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase),source=0.0_pReal)
+     plasticState(phase)%slipRate        => plasticState(phase)%dotState(2:2,1:NofMyPhase)
+     plasticState(phase)%accumulatedSlip => plasticState(phase)%state   (2:2,1:NofMyPhase)
+   endif myPhase
+ enddo initializeInstances
+
+end subroutine plastic_j2_init
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief calculates plastic velocity gradient and its tangent
+!--------------------------------------------------------------------------------------------------
+subroutine plastic_j2_LpAndItsTangent(Lp,dLp_dTstar99,Tstar_v,ipc,ip,el)
+ use math, only: &
+   math_mul6x6, &
+   math_Mandel6to33, &
+   math_Plain3333to99, &
+   math_deviatoric33, &
+   math_mul33xx33
+ use material, only: &
+   phaseAt, phasememberAt, &
+   plasticState, &
+   material_phase, &
+   phase_plasticityInstance
+
+ implicit none
+ real(pReal), dimension(3,3), intent(out) :: &
+   Lp                                                                                               !< plastic velocity gradient
+ real(pReal), dimension(9,9), intent(out) :: &
+   dLp_dTstar99                                                                                     !< derivative of Lp with respect to 2nd Piola Kirchhoff stress
+
+ real(pReal), dimension(6),   intent(in) :: &
+   Tstar_v                                                                                          !< 2nd Piola Kirchhoff stress tensor in Mandel notation
+ integer(pInt),               intent(in) :: &
+   ipc, &                                                                                           !< component-ID of integration point
+   ip, &                                                                                            !< integration point
+   el                                                                                               !< element
+
+ real(pReal), dimension(3,3) :: &
+   Tstar_dev_33                                                                                     !< deviatoric part of the 2nd Piola Kirchhoff stress tensor as 2nd order tensor
+ real(pReal), dimension(3,3,3,3) :: &
+   dLp_dTstar_3333                                                                                  !< derivative of Lp with respect to Tstar as 4th order tensor
+ real(pReal) :: &
+   gamma_dot, &                                                                                     !< strainrate
+   norm_Tstar_dev, &                                                                                !< euclidean norm of Tstar_dev
+   squarenorm_Tstar_dev                                                                             !< square of the euclidean norm of Tstar_dev
+ integer(pInt) :: &
+   instance, &
+   k, l, m, n
+
+ instance = phase_plasticityInstance(material_phase(ipc,ip,el))
+ Tstar_dev_33 = math_deviatoric33(math_Mandel6to33(Tstar_v))                                        ! deviatoric part of 2nd Piola-Kirchhoff stress
+ squarenorm_Tstar_dev = math_mul33xx33(Tstar_dev_33,Tstar_dev_33)
+ norm_Tstar_dev = sqrt(squarenorm_Tstar_dev) 
+
+ if (norm_Tstar_dev <= 0.0_pReal) then                                                              ! Tstar == 0 --> both Lp and dLp_dTstar are zero
+   Lp = 0.0_pReal
+   dLp_dTstar99 = 0.0_pReal
+ else
+   gamma_dot = plastic_j2_gdot0(instance) &
+             * (sqrt(1.5_pReal) * norm_Tstar_dev / (plastic_j2_fTaylor(instance) * &
+               plasticState(phaseAt(ipc,ip,el))%state(1,phasememberAt(ipc,ip,el)))) &
+                                                  **plastic_j2_n(instance)
+
+   Lp = Tstar_dev_33/norm_Tstar_dev * gamma_dot/plastic_j2_fTaylor(instance)
+
+!--------------------------------------------------------------------------------------------------
+! Calculation of the tangent of Lp
+   forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt,m=1_pInt:3_pInt,n=1_pInt:3_pInt) &
+     dLp_dTstar_3333(k,l,m,n) = (plastic_j2_n(instance)-1.0_pReal) * &
+                                      Tstar_dev_33(k,l)*Tstar_dev_33(m,n) / squarenorm_Tstar_dev
+   forall (k=1_pInt:3_pInt,l=1_pInt:3_pInt) &
+     dLp_dTstar_3333(k,l,k,l) = dLp_dTstar_3333(k,l,k,l) + 1.0_pReal
+   forall (k=1_pInt:3_pInt,m=1_pInt:3_pInt) &
+     dLp_dTstar_3333(k,k,m,m) = dLp_dTstar_3333(k,k,m,m) - 1.0_pReal/3.0_pReal
+   dLp_dTstar99 = math_Plain3333to99(gamma_dot / plastic_j2_fTaylor(instance) * &
+                                      dLp_dTstar_3333 / norm_Tstar_dev)
+ end if
+end subroutine plastic_j2_LpAndItsTangent
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief calculates the rate of change of microstructure
+!--------------------------------------------------------------------------------------------------
+subroutine plastic_j2_dotState(Tstar_v,ipc,ip,el)
+ use math, only: &
+   math_mul6x6
+ use material, only: &
+   phaseAt, phasememberAt, &
+   plasticState, &
+   material_phase, &
+   phase_plasticityInstance
+ 
+ implicit none
+ real(pReal), dimension(6), intent(in):: &
+   Tstar_v                                                                                          !< 2nd Piola Kirchhoff stress tensor in Mandel notation
+ integer(pInt),             intent(in) :: &
+   ipc, &                                                                                           !< component-ID of integration point
+   ip, &                                                                                            !< integration point
+   el                                                                                               !< element
+ real(pReal), dimension(6) :: &
+   Tstar_dev_v                                                                                      !< deviatoric part of the 2nd Piola Kirchhoff stress tensor in Mandel notation
+ real(pReal) :: &
+   gamma_dot, &                                                                                     !< strainrate
+   hardening, &                                                                                     !< hardening coefficient
+   saturation, &                                                                                    !< saturation resistance
+   norm_Tstar_dev                                                                                   !< euclidean norm of Tstar_dev
+ integer(pInt) :: &
+   instance, &                                                                                      !< instance of my instance (unique number of my constitutive model)
+   of, &                                                                                            !< shortcut notation for offset position in state array
+   ph                                                                                               !< shortcut notation for phase ID (unique number of all phases, regardless of constitutive model)
+
+ of = phasememberAt(ipc,ip,el)
+ ph = phaseAt(ipc,ip,el)
+ instance = phase_plasticityInstance(material_phase(ipc,ip,el))
+
+!--------------------------------------------------------------------------------------------------
+! norm of deviatoric part of 2nd Piola-Kirchhoff stress
+ Tstar_dev_v(1:3) = Tstar_v(1:3) - sum(Tstar_v(1:3))/3.0_pReal
+ Tstar_dev_v(4:6) = Tstar_v(4:6)
+ norm_Tstar_dev = sqrt(math_mul6x6(Tstar_dev_v,Tstar_dev_v))
+
+!--------------------------------------------------------------------------------------------------
+! strain rate 
+ gamma_dot = plastic_j2_gdot0(instance) * ( sqrt(1.5_pReal) * norm_Tstar_dev & 
+            / &!-----------------------------------------------------------------------------------
+           (plastic_j2_fTaylor(instance)*plasticState(ph)%state(1,of)) )**plastic_j2_n(instance)
+ 
+!--------------------------------------------------------------------------------------------------
+! hardening coefficient
+ if (abs(gamma_dot) > 1e-12_pReal) then
+   if (abs(plastic_j2_tausat_SinhFitA(instance)) <= tiny(0.0_pReal)) then
+     saturation = plastic_j2_tausat(instance)
+   else
+     saturation = (  plastic_j2_tausat(instance) &
+                   + ( log(  ( gamma_dot / plastic_j2_tausat_SinhFitA(instance)&
+                               )**(1.0_pReal / plastic_j2_tausat_SinhFitD(instance))&
+                            + sqrt(  ( gamma_dot / plastic_j2_tausat_SinhFitA(instance) &
+                                      )**(2.0_pReal / plastic_j2_tausat_SinhFitD(instance)) &
+                                   + 1.0_pReal ) &
+                            ) & ! asinh(K) = ln(K + sqrt(K^2 +1))
+                       )**(1.0_pReal / plastic_j2_tausat_SinhFitC(instance)) &
+                   / (  plastic_j2_tausat_SinhFitB(instance) &
+                      * (gamma_dot / plastic_j2_gdot0(instance))**(1.0_pReal / plastic_j2_n(instance)) &
+                      ) &
+                   )
+   endif
+   hardening = ( plastic_j2_h0(instance) + plastic_j2_h0_slopeLnRate(instance) * log(gamma_dot) ) &
+               * abs( 1.0_pReal - plasticState(ph)%state(1,of)/saturation )**plastic_j2_a(instance) &
+               * sign(1.0_pReal, 1.0_pReal - plasticState(ph)%state(1,of)/saturation)
+ else
+   hardening = 0.0_pReal
+ endif
+
+  plasticState(ph)%dotState(1,of) = hardening * gamma_dot
+  plasticState(ph)%dotState(2,of) =             gamma_dot
+
+end subroutine plastic_j2_dotState
+
+!--------------------------------------------------------------------------------------------------
+!> @brief return array of constitutive results
+!--------------------------------------------------------------------------------------------------
+function plastic_j2_postResults(Tstar_v,ipc,ip,el)
+ use math, only: &
+   math_mul6x6
+ use material, only: &
+   material_phase, &
+   plasticState, &
+   phaseAt, phasememberAt, &
+   phase_plasticityInstance
+
+ implicit none
+ real(pReal), dimension(6),  intent(in) :: &
+   Tstar_v                                                                                          !< 2nd Piola Kirchhoff stress tensor in Mandel notation
+ integer(pInt),              intent(in) :: &
+   ipc, &                                                                                           !< component-ID of integration point
+   ip, &                                                                                            !< integration point
+   el                                                                                               !< element
+ real(pReal), dimension(plastic_j2_sizePostResults(phase_plasticityInstance(material_phase(ipc,ip,el)))) :: &
+                                           plastic_j2_postResults
+
+ real(pReal), dimension(6) :: &
+   Tstar_dev_v                                                                                      ! deviatoric part of the 2nd Piola Kirchhoff stress tensor in Mandel notation
+ real(pReal) :: &
+   norm_Tstar_dev                                                                                   ! euclidean norm of Tstar_dev
+ integer(pInt) :: &
+   instance, &                                                                                      !< instance of my instance (unique number of my constitutive model)
+   of, &                                                                                            !< shortcut notation for offset position in state array
+   ph, &                                                                                            !< shortcut notation for phase ID (unique number of all phases, regardless of constitutive model)
+   c, &
+   o
+
+ of = phasememberAt(ipc,ip,el)
+ ph = phaseAt(ipc,ip,el)
+ instance = phase_plasticityInstance(material_phase(ipc,ip,el))
+ 
+!--------------------------------------------------------------------------------------------------
+! calculate deviatoric part of 2nd Piola-Kirchhoff stress and its norm
+ Tstar_dev_v(1:3) = Tstar_v(1:3) - sum(Tstar_v(1:3))/3.0_pReal
+ Tstar_dev_v(4:6) = Tstar_v(4:6)
+ norm_Tstar_dev = sqrt(math_mul6x6(Tstar_dev_v,Tstar_dev_v))
+ 
+ c = 0_pInt
+ plastic_j2_postResults = 0.0_pReal
+
+ outputsLoop: do o = 1_pInt,plastic_j2_Noutput(instance)
+   select case(plastic_j2_outputID(o,instance))
+     case (flowstress_ID)
+       plastic_j2_postResults(c+1_pInt) = plasticState(ph)%state(1,of)
+       c = c + 1_pInt
+     case (strainrate_ID)
+       plastic_j2_postResults(c+1_pInt) = &
+                plastic_j2_gdot0(instance) * (            sqrt(1.5_pReal) * norm_Tstar_dev & 
+             / &!----------------------------------------------------------------------------------
+              (plastic_j2_fTaylor(instance) * plasticState(ph)%state(1,of)) ) ** plastic_j2_n(instance)
+       c = c + 1_pInt
+   end select
+ enddo outputsLoop
+
+end function plastic_j2_postResults
+
+
+end module plastic_j2

code/plastic_none.f90

@@ -0,0 +1,109 @@
+!--------------------------------------------------------------------------------------------------
+! $Id$
+!--------------------------------------------------------------------------------------------------
+!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
+!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
+!> @brief material subroutine for purely elastic material
+!--------------------------------------------------------------------------------------------------
+module plastic_none
+ use prec, only: &
+   pInt
+
+ implicit none
+ private
+ integer(pInt),                       dimension(:),     allocatable,          public, protected :: &
+   plastic_none_sizePostResults
+
+ integer(pInt),                       dimension(:,:),   allocatable, target,  public :: &
+   plastic_none_sizePostResult                                                                 !< size of each post result output
+
+ public :: &
+   plastic_none_init
+
+contains
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief module initialization
+!> @details reads in material parameters, allocates arrays, and does sanity checks
+!--------------------------------------------------------------------------------------------------
+subroutine plastic_none_init
+ use, intrinsic :: iso_fortran_env                                                                  ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
+ use debug, only: &
+   debug_level, &
+   debug_constitutive, &
+   debug_levelBasic
+ use IO, only: &
+   IO_timeStamp
+ use numerics, only: &
+   worldrank, &
+   numerics_integrator
+ use material, only: &
+   phase_plasticity, &
+   PLASTICITY_NONE_label, &
+   material_phase, &
+   plasticState, &
+   PLASTICITY_none_ID
+
+ implicit none
+
+ integer(pInt) :: &
+   maxNinstance, &
+   phase, &
+   NofMyPhase, &
+   sizeState, &
+   sizeDotState, &
+   sizeDeltaState
+ 
+ mainProcess: if (worldrank == 0) then 
+   write(6,'(/,a)')   ' <<<+-  constitutive_'//PLASTICITY_NONE_label//' init  -+>>>'
+   write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
+#include "compilation_info.f90"
+ endif mainProcess
+ 
+ maxNinstance = int(count(phase_plasticity == PLASTICITY_none_ID),pInt)
+ if (maxNinstance == 0_pInt) return
+
+ if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
+   write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
+
+ initializeInstances: do phase = 1_pInt, size(phase_plasticity)
+   if (phase_plasticity(phase) == PLASTICITY_none_ID) then
+   NofMyPhase=count(material_phase==phase)
+
+     sizeState    = 0_pInt
+     plasticState(phase)%sizeState = sizeState
+     sizeDotState = sizeState
+     plasticState(phase)%sizeDotState = sizeDotState
+     sizeDeltaState = 0_pInt
+     plasticState(phase)%sizeDeltaState = sizeDeltaState
+     plasticState(phase)%sizePostResults = 0_pInt
+     plasticState(phase)%nSlip  = 0_pInt
+     plasticState(phase)%nTwin  = 0_pInt
+     plasticState(phase)%nTrans = 0_pInt
+     allocate(plasticState(phase)%aTolState          (sizeState))
+     allocate(plasticState(phase)%state0             (sizeState,NofMyPhase))
+     allocate(plasticState(phase)%partionedState0    (sizeState,NofMyPhase))
+     allocate(plasticState(phase)%subState0          (sizeState,NofMyPhase))
+     allocate(plasticState(phase)%state              (sizeState,NofMyPhase))
+     allocate(plasticState(phase)%state_backup       (sizeState,NofMyPhase))
+
+     allocate(plasticState(phase)%dotState           (sizeDotState,NofMyPhase))
+     allocate(plasticState(phase)%deltaState        (sizeDeltaState,NofMyPhase))
+     allocate(plasticState(phase)%dotState_backup    (sizeDotState,NofMyPhase))
+     if (any(numerics_integrator == 1_pInt)) then
+       allocate(plasticState(phase)%previousDotState (sizeDotState,NofMyPhase))
+       allocate(plasticState(phase)%previousDotState2(sizeDotState,NofMyPhase))
+     endif
+     if (any(numerics_integrator == 4_pInt)) &
+       allocate(plasticState(phase)%RK4dotState      (sizeDotState,NofMyPhase))
+     if (any(numerics_integrator == 5_pInt)) &
+       allocate(plasticState(phase)%RKCK45dotState (6,sizeDotState,NofMyPhase))
+   endif
+ enddo initializeInstances
+
+ allocate(plastic_none_sizePostResults(maxNinstance), source=0_pInt)
+
+end subroutine plastic_none_init
+
+end module plastic_none

code/source_damage_anisoBrittle.f90

@@ -0,0 +1,425 @@
+!--------------------------------------------------------------------------------------------------
+! $Id$
+!--------------------------------------------------------------------------------------------------
+!> @author Luv Sharma, Max-Planck-Institut fŸr Eisenforschung GmbH
+!> @author Pratheek Shanthraj, Max-Planck-Institut fŸr Eisenforschung GmbH
+!> @brief material subroutine incorporating anisotropic brittle damage source mechanism
+!> @details to be done
+!--------------------------------------------------------------------------------------------------
+module source_damage_anisoBrittle
+ use prec, only: &
+   pReal, &
+   pInt
+
+ implicit none
+ private
+ integer(pInt),                       dimension(:),           allocatable,         public, protected :: &
+   source_damage_anisoBrittle_sizePostResults, &                                                                !< cumulative size of post results
+   source_damage_anisoBrittle_offset, &                                                                         !< which source is my current source mechanism?
+   source_damage_anisoBrittle_instance                                                                          !< instance of source mechanism
+
+ integer(pInt),                       dimension(:,:),         allocatable, target, public  :: &
+   source_damage_anisoBrittle_sizePostResult                                                                    !< size of each post result output
+
+ character(len=64),                   dimension(:,:),         allocatable, target, public  :: &
+   source_damage_anisoBrittle_output                                                                            !< name of each post result output
+   
+ integer(pInt),                       dimension(:),           allocatable, target, public  :: &
+   source_damage_anisoBrittle_Noutput                                                                           !< number of outputs per instance of this source 
+
+ integer(pInt),                       dimension(:),           allocatable,         private :: &
+   source_damage_anisoBrittle_totalNcleavage                                                                    !< total number of cleavage systems
+   
+ integer(pInt),                       dimension(:,:),         allocatable,         private :: &
+   source_damage_anisoBrittle_Ncleavage                                                                         !< number of cleavage systems per family
+   
+ real(pReal),                         dimension(:),           allocatable,         private :: &
+   source_damage_anisoBrittle_aTol, &
+   source_damage_anisoBrittle_sdot_0, &
+   source_damage_anisoBrittle_N
+
+ real(pReal),                         dimension(:,:),         allocatable,         private :: &
+   source_damage_anisoBrittle_critDisp, &
+   source_damage_anisoBrittle_critLoad
+
+ enum, bind(c) 
+   enumerator :: undefined_ID, &
+                 damage_drivingforce_ID
+ end enum                                                
+ 
+ integer(kind(undefined_ID)),         dimension(:,:),         allocatable,          private :: & 
+   source_damage_anisoBrittle_outputID                                                                  !< ID of each post result output
+
+
+ public :: &
+   source_damage_anisoBrittle_init, &
+   source_damage_anisoBrittle_dotState, &
+   source_damage_anisobrittle_getRateAndItsTangent, &
+   source_damage_anisoBrittle_postResults
+
+contains
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief module initialization
+!> @details reads in material parameters, allocates arrays, and does sanity checks
+!--------------------------------------------------------------------------------------------------
+subroutine source_damage_anisoBrittle_init(fileUnit)
+ use, intrinsic :: iso_fortran_env                                                                  ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
+ use debug, only: &
+   debug_level,&
+   debug_constitutive,&
+   debug_levelBasic
+ use IO, only: &
+   IO_read, &
+   IO_lc, &
+   IO_getTag, &
+   IO_isBlank, &
+   IO_stringPos, &
+   IO_stringValue, &
+   IO_floatValue, &
+   IO_intValue, &
+   IO_warning, &
+   IO_error, &
+   IO_timeStamp, &
+   IO_EOF
+ use material, only: &
+   phase_source, &
+   phase_Nsources, &
+   phase_Noutput, &
+   SOURCE_damage_anisoBrittle_label, &
+   SOURCE_damage_anisoBrittle_ID, &
+   material_Nphase, &
+   material_phase, &  
+   sourceState, &
+   MATERIAL_partPhase
+ use numerics,only: &
+   analyticJaco, &
+   worldrank, &
+   numerics_integrator
+ use lattice, only: &
+   lattice_maxNcleavageFamily, &
+   lattice_NcleavageSystem
+
+ implicit none
+ integer(pInt), intent(in) :: fileUnit
+
+ integer(pInt), allocatable, dimension(:) :: chunkPos
+ integer(pInt) :: maxNinstance,mySize=0_pInt,phase,instance,source,sourceOffset,o
+ integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
+ integer(pInt) :: NofMyPhase   
+ integer(pInt) :: Nchunks_CleavageFamilies = 0_pInt, j   
+ character(len=65536) :: &
+   tag  = '', &
+   line = ''
+
+ mainProcess: if (worldrank == 0) then 
+   write(6,'(/,a)')   ' <<<+-  source_'//SOURCE_damage_anisoBrittle_LABEL//' init  -+>>>'
+   write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
+#include "compilation_info.f90"
+ endif mainProcess
+
+ maxNinstance = int(count(phase_source == SOURCE_damage_anisoBrittle_ID),pInt)
+ if (maxNinstance == 0_pInt) return
+ 
+ if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
+   write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
+ 
+ allocate(source_damage_anisoBrittle_offset(material_Nphase), source=0_pInt)
+ allocate(source_damage_anisoBrittle_instance(material_Nphase), source=0_pInt)
+ do phase = 1, material_Nphase
+   source_damage_anisoBrittle_instance(phase) = count(phase_source(:,1:phase) == source_damage_anisoBrittle_ID)
+   do source = 1, phase_Nsources(phase)
+     if (phase_source(source,phase) == source_damage_anisoBrittle_ID) &
+       source_damage_anisoBrittle_offset(phase) = source
+   enddo    
+ enddo
+   
+ allocate(source_damage_anisoBrittle_sizePostResults(maxNinstance),                      source=0_pInt)
+ allocate(source_damage_anisoBrittle_sizePostResult(maxval(phase_Noutput),maxNinstance), source=0_pInt)
+ allocate(source_damage_anisoBrittle_output(maxval(phase_Noutput),maxNinstance))
+          source_damage_anisoBrittle_output = ''
+ allocate(source_damage_anisoBrittle_outputID(maxval(phase_Noutput),maxNinstance),       source=undefined_ID)
+ allocate(source_damage_anisoBrittle_Noutput(maxNinstance),                              source=0_pInt) 
+ allocate(source_damage_anisoBrittle_critDisp(lattice_maxNcleavageFamily,maxNinstance),  source=0.0_pReal) 
+ allocate(source_damage_anisoBrittle_critLoad(lattice_maxNcleavageFamily,maxNinstance),  source=0.0_pReal) 
+ allocate(source_damage_anisoBrittle_Ncleavage(lattice_maxNcleavageFamily,maxNinstance), source=0_pInt)
+ allocate(source_damage_anisoBrittle_totalNcleavage(maxNinstance),                       source=0_pInt)
+ allocate(source_damage_anisoBrittle_aTol(maxNinstance),                                 source=0.0_pReal) 
+ allocate(source_damage_anisoBrittle_sdot_0(maxNinstance),                               source=0.0_pReal) 
+ allocate(source_damage_anisoBrittle_N(maxNinstance),                                    source=0.0_pReal) 
+
+ rewind(fileUnit)
+ phase = 0_pInt
+ do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase)         ! wind forward to <phase>
+   line = IO_read(fileUnit)
+ enddo
+ 
+ parsingFile: do while (trim(line) /= IO_EOF)                                                       ! read through sections of phase part
+   line = IO_read(fileUnit)
+   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
+   if (IO_getTag(line,'<','>') /= '') then                                                          ! stop at next part
+     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
+     exit                                                                                           
+   endif   
+   if (IO_getTag(line,'[',']') /= '') then                                                          ! next phase section
+     phase = phase + 1_pInt                                                                         ! advance phase section counter
+     cycle                                                                                          ! skip to next line
+   endif
+   if (phase > 0_pInt ) then; if (any(phase_source(:,phase) == SOURCE_damage_anisoBrittle_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
+     instance = source_damage_anisoBrittle_instance(phase)                                          ! which instance of my damage is present phase
+     chunkPos = IO_stringPos(line)
+     tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt))                                             ! extract key
+     select case(tag)
+       case ('(output)')
+         select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt)))
+           case ('anisobrittle_drivingforce')
+             source_damage_anisoBrittle_Noutput(instance) = source_damage_anisoBrittle_Noutput(instance) + 1_pInt
+             source_damage_anisoBrittle_outputID(source_damage_anisoBrittle_Noutput(instance),instance) = damage_drivingforce_ID
+             source_damage_anisoBrittle_output(source_damage_anisoBrittle_Noutput(instance),instance) = &
+                                                       IO_lc(IO_stringValue(line,chunkPos,2_pInt))
+          end select
+
+       case ('anisobrittle_atol')
+         source_damage_anisoBrittle_aTol(instance) = IO_floatValue(line,chunkPos,2_pInt)
+         
+       case ('anisobrittle_sdot0')
+         source_damage_anisoBrittle_sdot_0(instance) = IO_floatValue(line,chunkPos,2_pInt)
+         
+       case ('anisobrittle_ratesensitivity')
+         source_damage_anisoBrittle_N(instance) = IO_floatValue(line,chunkPos,2_pInt)
+         
+       case ('ncleavage')  !
+         Nchunks_CleavageFamilies = chunkPos(1) - 1_pInt
+         do j = 1_pInt, Nchunks_CleavageFamilies
+           source_damage_anisoBrittle_Ncleavage(j,instance) = IO_intValue(line,chunkPos,1_pInt+j)
+         enddo
+
+       case ('anisobrittle_criticaldisplacement')
+         do j = 1_pInt, Nchunks_CleavageFamilies
+           source_damage_anisoBrittle_critDisp(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j)
+         enddo
+
+       case ('anisobrittle_criticalload')
+         do j = 1_pInt, Nchunks_CleavageFamilies
+           source_damage_anisoBrittle_critLoad(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j)
+         enddo
+
+     end select
+   endif; endif
+ enddo parsingFile
+
+!--------------------------------------------------------------------------------------------------
+!  sanity checks
+ sanityChecks: do phase = 1_pInt, material_Nphase  
+   myPhase: if (any(phase_source(:,phase) == SOURCE_damage_anisoBrittle_ID)) then
+     instance = source_damage_anisoBrittle_instance(phase)
+     source_damage_anisoBrittle_Ncleavage(1:lattice_maxNcleavageFamily,instance) = &
+       min(lattice_NcleavageSystem(1:lattice_maxNcleavageFamily,phase),&                            ! limit active cleavage systems per family to min of available and requested
+           source_damage_anisoBrittle_Ncleavage(1:lattice_maxNcleavageFamily,instance))
+     source_damage_anisoBrittle_totalNcleavage(instance)  = sum(source_damage_anisoBrittle_Ncleavage(:,instance)) ! how many cleavage systems altogether
+     if (source_damage_anisoBrittle_aTol(instance) < 0.0_pReal) &
+       source_damage_anisoBrittle_aTol(instance) = 1.0e-3_pReal                                     ! default absolute tolerance 1e-3
+     if (source_damage_anisoBrittle_sdot_0(instance) <= 0.0_pReal) &
+       call IO_error(211_pInt,el=instance,ext_msg='sdot_0 ('//SOURCE_damage_anisoBrittle_LABEL//')')
+     if (any(source_damage_anisoBrittle_critDisp(1:Nchunks_CleavageFamilies,instance) < 0.0_pReal)) &
+       call IO_error(211_pInt,el=instance,ext_msg='critical_displacement ('//SOURCE_damage_anisoBrittle_LABEL//')')
+     if (any(source_damage_anisoBrittle_critLoad(1:Nchunks_CleavageFamilies,instance) < 0.0_pReal)) &
+       call IO_error(211_pInt,el=instance,ext_msg='critical_load ('//SOURCE_damage_anisoBrittle_LABEL//')')
+     if (source_damage_anisoBrittle_N(instance) <= 0.0_pReal) &
+       call IO_error(211_pInt,el=instance,ext_msg='rate_sensitivity ('//SOURCE_damage_anisoBrittle_LABEL//')')
+   endif myPhase
+ enddo sanityChecks
+ 
+ initializeInstances: do phase = 1_pInt, material_Nphase
+   if (any(phase_source(:,phase) == SOURCE_damage_anisoBrittle_ID)) then
+     NofMyPhase=count(material_phase==phase)
+     instance = source_damage_anisoBrittle_instance(phase)
+     sourceOffset = source_damage_anisoBrittle_offset(phase)
+
+!--------------------------------------------------------------------------------------------------
+!  Determine size of postResults array
+     outputsLoop: do o = 1_pInt,source_damage_anisoBrittle_Noutput(instance)
+       select case(source_damage_anisoBrittle_outputID(o,instance))
+         case(damage_drivingforce_ID)
+           mySize = 1_pInt
+       end select
+ 
+       if (mySize > 0_pInt) then  ! any meaningful output found
+          source_damage_anisoBrittle_sizePostResult(o,instance) = mySize
+          source_damage_anisoBrittle_sizePostResults(instance)  = source_damage_anisoBrittle_sizePostResults(instance) + mySize
+       endif
+     enddo outputsLoop
+
+!--------------------------------------------------------------------------------------------------
+! Determine size of state array
+     sizeDotState              = 1_pInt
+     sizeDeltaState            = 0_pInt
+     sizeState                 = 1_pInt
+
+     sourceState(phase)%p(sourceOffset)%sizeState      = sizeState
+     sourceState(phase)%p(sourceOffset)%sizeDotState   = sizeDotState
+     sourceState(phase)%p(sourceOffset)%sizeDeltaState = sizeDeltaState
+     sourceState(phase)%p(sourceOffset)%sizePostResults = source_damage_anisoBrittle_sizePostResults(instance)
+     allocate(sourceState(phase)%p(sourceOffset)%aTolState           (sizeState),                &
+              source=source_damage_anisoBrittle_aTol(instance))
+     allocate(sourceState(phase)%p(sourceOffset)%state0              (sizeState,NofMyPhase),     source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%partionedState0     (sizeState,NofMyPhase),     source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%subState0           (sizeState,NofMyPhase),     source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%state               (sizeState,NofMyPhase),     source=0.0_pReal)
+
+     allocate(sourceState(phase)%p(sourceOffset)%dotState            (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%deltaState        (sizeDeltaState,NofMyPhase),  source=0.0_pReal)
+     if (.not. analyticJaco) then
+       allocate(sourceState(phase)%p(sourceOffset)%state_backup      (sizeState,NofMyPhase),     source=0.0_pReal)
+       allocate(sourceState(phase)%p(sourceOffset)%dotState_backup   (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     endif
+     if (any(numerics_integrator == 1_pInt)) then
+       allocate(sourceState(phase)%p(sourceOffset)%previousDotState  (sizeDotState,NofMyPhase),  source=0.0_pReal)
+       allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     endif
+     if (any(numerics_integrator == 4_pInt)) &
+       allocate(sourceState(phase)%p(sourceOffset)%RK4dotState       (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     if (any(numerics_integrator == 5_pInt)) &
+       allocate(sourceState(phase)%p(sourceOffset)%RKCK45dotState    (6,sizeDotState,NofMyPhase),source=0.0_pReal)
+
+   endif
+ 
+ enddo initializeInstances
+end subroutine source_damage_anisoBrittle_init
+
+!--------------------------------------------------------------------------------------------------
+!> @brief calculates derived quantities from state
+!--------------------------------------------------------------------------------------------------
+subroutine source_damage_anisoBrittle_dotState(Tstar_v, ipc, ip, el)
+ use material, only: &
+   phaseAt, phasememberAt, &
+   sourceState, &
+   material_homog, &
+   damage, &
+   damageMapping
+ use lattice, only: &
+   lattice_Scleavage_v, &
+   lattice_maxNcleavageFamily, &
+   lattice_NcleavageSystem
+
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< component-ID of integration point
+   ip, &                                                                                            !< integration point
+   el                                                                                               !< element
+ real(pReal),  intent(in), dimension(6) :: &
+   Tstar_v                                                                                          !< 2nd Piola Kirchhoff stress tensor (Mandel)
+ integer(pInt) :: &
+   phase, &
+   constituent, &
+   instance, &
+   sourceOffset, &
+   damageOffset, &
+   homog, &
+   f, i, index_myFamily
+ real(pReal) :: &
+   traction_d, traction_t, traction_n, traction_crit
+
+ phase = phaseAt(ipc,ip,el)
+ constituent = phasememberAt(ipc,ip,el)
+ instance = source_damage_anisoBrittle_instance(phase)
+ sourceOffset = source_damage_anisoBrittle_offset(phase)
+ homog = material_homog(ip,el)
+ damageOffset = damageMapping(homog)%p(ip,el)
+ 
+ sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = 0.0_pReal
+ do f = 1_pInt,lattice_maxNcleavageFamily
+   index_myFamily = sum(lattice_NcleavageSystem(1:f-1_pInt,phase))                                    ! at which index starts my family
+   do i = 1_pInt,source_damage_anisoBrittle_Ncleavage(f,instance)                                     ! process each (active) cleavage system in family
+     traction_d    = dot_product(Tstar_v,lattice_Scleavage_v(1:6,1,index_myFamily+i,phase))
+     traction_t    = dot_product(Tstar_v,lattice_Scleavage_v(1:6,2,index_myFamily+i,phase))
+     traction_n    = dot_product(Tstar_v,lattice_Scleavage_v(1:6,3,index_myFamily+i,phase))
+     
+     traction_crit = source_damage_anisoBrittle_critLoad(f,instance)* &
+                     damage(homog)%p(damageOffset)*damage(homog)%p(damageOffset)
+     sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = &
+       sourceState(phase)%p(sourceOffset)%dotState(1,constituent) + &
+       source_damage_anisoBrittle_sdot_0(instance)* &
+       ((max(0.0_pReal, abs(traction_d) - traction_crit)/traction_crit)**source_damage_anisoBrittle_N(instance) + &
+        (max(0.0_pReal, abs(traction_t) - traction_crit)/traction_crit)**source_damage_anisoBrittle_N(instance) + &
+        (max(0.0_pReal, abs(traction_n) - traction_crit)/traction_crit)**source_damage_anisoBrittle_N(instance))/ &
+       source_damage_anisoBrittle_critDisp(f,instance)
+
+   enddo
+ enddo
+
+end subroutine source_damage_anisoBrittle_dotState
+
+!--------------------------------------------------------------------------------------------------
+!> @brief returns local part of nonlocal damage driving force
+!--------------------------------------------------------------------------------------------------
+subroutine source_damage_anisobrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ipc, ip,  el)
+ use material, only: &
+   phaseAt, phasememberAt, &
+   sourceState
+
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< component-ID of integration point
+   ip, &                                                                                            !< integration point
+   el                                                                                               !< element
+ real(pReal),  intent(in) :: &
+   phi
+ real(pReal),  intent(out) :: &
+   localphiDot, &
+   dLocalphiDot_dPhi
+ integer(pInt) :: &
+   phase, constituent, sourceOffset
+
+ phase = phaseAt(ipc,ip,el)
+ constituent = phasememberAt(ipc,ip,el)
+ sourceOffset = source_damage_anisoBrittle_offset(phase)
+ 
+ localphiDot = 1.0_pReal - &
+               sourceState(phase)%p(sourceOffset)%state(1,constituent)*phi
+ 
+ dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent)
+ 
+end subroutine source_damage_anisobrittle_getRateAndItsTangent
+ 
+!--------------------------------------------------------------------------------------------------
+!> @brief return array of local damage results
+!--------------------------------------------------------------------------------------------------
+function source_damage_anisoBrittle_postResults(ipc,ip,el)
+ use material, only: &
+   phaseAt, phasememberAt, &
+   sourceState
+
+ implicit none
+ integer(pInt),              intent(in) :: &
+   ipc, &                                                                                           !< component-ID of integration point
+   ip, &                                                                                            !< integration point
+   el                                                                                               !< element
+ real(pReal), dimension(source_damage_anisoBrittle_sizePostResults( &
+                          source_damage_anisoBrittle_instance(phaseAt(ipc,ip,el)))) :: &
+   source_damage_anisoBrittle_postResults
+
+ integer(pInt) :: &
+   instance, phase, constituent, sourceOffset, o, c
+   
+ phase = phaseAt(ipc,ip,el)
+ constituent = phasememberAt(ipc,ip,el)
+ instance = source_damage_anisoBrittle_instance(phase)
+ sourceOffset = source_damage_anisoBrittle_offset(phase)
+
+ c = 0_pInt
+ source_damage_anisoBrittle_postResults = 0.0_pReal
+
+ do o = 1_pInt,source_damage_anisoBrittle_Noutput(instance)
+    select case(source_damage_anisoBrittle_outputID(o,instance))
+      case (damage_drivingforce_ID)
+        source_damage_anisoBrittle_postResults(c+1_pInt) = &
+          sourceState(phase)%p(sourceOffset)%state(1,constituent)
+        c = c + 1_pInt
+
+    end select
+ enddo
+end function source_damage_anisoBrittle_postResults
+
+end module source_damage_anisoBrittle

code/source_damage_anisoDuctile.f90

@@ -0,0 +1,415 @@
+!--------------------------------------------------------------------------------------------------
+! $Id$
+!--------------------------------------------------------------------------------------------------
+!> @author Luv Sharma, Max-Planck-Institut für Eisenforschung GmbH
+!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
+!> @brief material subroutine incorporating anisotropic ductile damage source mechanism
+!> @details to be done
+!--------------------------------------------------------------------------------------------------
+module source_damage_anisoDuctile
+ use prec, only: &
+   pReal, &
+   pInt
+
+ implicit none
+ private
+ integer(pInt),                       dimension(:),           allocatable,         public, protected :: &
+   source_damage_anisoDuctile_sizePostResults, &                                                                !< cumulative size of post results
+   source_damage_anisoDuctile_offset, &                                                                         !< which source is my current damage mechanism?
+   source_damage_anisoDuctile_instance                                                                          !< instance of damage source mechanism
+
+ integer(pInt),                       dimension(:,:),         allocatable, target, public  :: &
+   source_damage_anisoDuctile_sizePostResult                                                                    !< size of each post result output
+
+ character(len=64),                   dimension(:,:),         allocatable, target, public  :: &
+   source_damage_anisoDuctile_output                                                                            !< name of each post result output
+   
+ integer(pInt),                       dimension(:),           allocatable, target, public  :: &
+   source_damage_anisoDuctile_Noutput                                                                           !< number of outputs per instance of this damage 
+   
+ integer(pInt),                       dimension(:),           allocatable,         private :: &
+   source_damage_anisoDuctile_totalNslip                                                                    !< total number of slip systems
+
+ integer(pInt),                       dimension(:,:),         allocatable,         private :: &
+   source_damage_anisoDuctile_Nslip                                                                         !< number of slip systems per family
+   
+ real(pReal),                         dimension(:),           allocatable,         private :: &
+   source_damage_anisoDuctile_aTol
+
+ real(pReal),                         dimension(:,:),         allocatable,         private :: &
+   source_damage_anisoDuctile_critPlasticStrain
+
+ real(pReal),                         dimension(:),           allocatable,         private :: &
+   source_damage_anisoDuctile_sdot_0, &
+   source_damage_anisoDuctile_N
+
+ real(pReal),                         dimension(:,:),         allocatable,         private :: &
+   source_damage_anisoDuctile_critLoad
+   
+ enum, bind(c) 
+   enumerator :: undefined_ID, &
+                 damage_drivingforce_ID
+ end enum 
+ 
+ integer(kind(undefined_ID)),         dimension(:,:),         allocatable,          private :: & 
+   source_damage_anisoDuctile_outputID                                                                  !< ID of each post result output
+
+
+ public :: &
+   source_damage_anisoDuctile_init, &
+   source_damage_anisoDuctile_dotState, &
+   source_damage_anisoDuctile_getRateAndItsTangent, &
+   source_damage_anisoDuctile_postResults
+
+contains
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief module initialization
+!> @details reads in material parameters, allocates arrays, and does sanity checks
+!--------------------------------------------------------------------------------------------------
+subroutine source_damage_anisoDuctile_init(fileUnit)
+ use, intrinsic :: iso_fortran_env                                                                  ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
+ use debug, only: &
+   debug_level,&
+   debug_constitutive,&
+   debug_levelBasic
+ use IO, only: &
+   IO_read, &
+   IO_lc, &
+   IO_getTag, &
+   IO_isBlank, &
+   IO_stringPos, &
+   IO_stringValue, &
+   IO_floatValue, &
+   IO_intValue, &
+   IO_warning, &
+   IO_error, &
+   IO_timeStamp, &
+   IO_EOF
+ use material, only: &
+   phase_source, &
+   phase_Nsources, &
+   phase_Noutput, &
+   SOURCE_damage_anisoDuctile_label, &
+   SOURCE_damage_anisoDuctile_ID, &
+   material_Nphase, &
+   material_phase, &  
+   sourceState, &
+   MATERIAL_partPhase
+ use numerics,only: &
+   analyticJaco, &
+   worldrank, &
+   numerics_integrator
+ use lattice, only: &
+   lattice_maxNslipFamily, &
+   lattice_NslipSystem
+
+ implicit none
+ integer(pInt), intent(in) :: fileUnit
+
+ integer(pInt), allocatable, dimension(:) :: chunkPos
+ integer(pInt) :: maxNinstance,mySize=0_pInt,phase,instance,source,sourceOffset,o
+ integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
+ integer(pInt) :: NofMyPhase   
+ integer(pInt) :: Nchunks_SlipFamilies = 0_pInt, j   
+ character(len=65536) :: &
+   tag  = '', &
+   line = ''
+
+ mainProcess: if (worldrank == 0) then 
+   write(6,'(/,a)')   ' <<<+-  source_'//SOURCE_damage_anisoDuctile_LABEL//' init  -+>>>'
+   write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
+#include "compilation_info.f90"
+ endif mainProcess
+
+ maxNinstance = int(count(phase_source == SOURCE_damage_anisoDuctile_ID),pInt)
+ if (maxNinstance == 0_pInt) return
+ 
+ if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
+   write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
+ 
+ allocate(source_damage_anisoDuctile_offset(material_Nphase), source=0_pInt)
+ allocate(source_damage_anisoDuctile_instance(material_Nphase), source=0_pInt)
+ do phase = 1, material_Nphase
+   source_damage_anisoDuctile_instance(phase) = count(phase_source(:,1:phase) == source_damage_anisoDuctile_ID)
+   do source = 1, phase_Nsources(phase)
+     if (phase_source(source,phase) == source_damage_anisoDuctile_ID) &
+       source_damage_anisoDuctile_offset(phase) = source
+   enddo    
+ enddo
+   
+ allocate(source_damage_anisoDuctile_sizePostResults(maxNinstance),                     source=0_pInt)
+ allocate(source_damage_anisoDuctile_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
+ allocate(source_damage_anisoDuctile_output(maxval(phase_Noutput),maxNinstance))
+          source_damage_anisoDuctile_output = ''
+ allocate(source_damage_anisoDuctile_outputID(maxval(phase_Noutput),maxNinstance),      source=undefined_ID)
+ allocate(source_damage_anisoDuctile_Noutput(maxNinstance),                             source=0_pInt) 
+ allocate(source_damage_anisoDuctile_critLoad(lattice_maxNslipFamily,maxNinstance), source=0.0_pReal) 
+ allocate(source_damage_anisoDuctile_critPlasticStrain(lattice_maxNslipFamily,maxNinstance),source=0.0_pReal) 
+ allocate(source_damage_anisoDuctile_Nslip(lattice_maxNslipFamily,maxNinstance),        source=0_pInt)
+ allocate(source_damage_anisoDuctile_totalNslip(maxNinstance),                          source=0_pInt)
+ allocate(source_damage_anisoDuctile_N(maxNinstance),                                   source=0.0_pReal) 
+ allocate(source_damage_anisoDuctile_sdot_0(maxNinstance),                              source=0.0_pReal) 
+ allocate(source_damage_anisoDuctile_aTol(maxNinstance),                                source=0.0_pReal) 
+
+ rewind(fileUnit)
+ phase = 0_pInt
+ do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase)         ! wind forward to <phase>
+   line = IO_read(fileUnit)
+ enddo
+ 
+ parsingFile: do while (trim(line) /= IO_EOF)                                                       ! read through sections of phase part
+   line = IO_read(fileUnit)
+   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
+   if (IO_getTag(line,'<','>') /= '') then                                                          ! stop at next part
+     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
+     exit                                                                                           
+   endif   
+   if (IO_getTag(line,'[',']') /= '') then                                                          ! next phase section
+     phase = phase + 1_pInt                                                                         ! advance phase section counter
+     cycle                                                                                          ! skip to next line
+   endif
+   if (phase > 0_pInt ) then; if (any(phase_source(:,phase) == SOURCE_damage_anisoDuctile_ID)) then         ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
+     instance = source_damage_anisoDuctile_instance(phase)                                                         ! which instance of my damage is present phase
+     chunkPos = IO_stringPos(line)
+     tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt))                                             ! extract key
+     select case(tag)
+       case ('(output)')
+         select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt)))
+           case ('anisoductile_drivingforce')
+             source_damage_anisoDuctile_Noutput(instance) = source_damage_anisoDuctile_Noutput(instance) + 1_pInt
+             source_damage_anisoDuctile_outputID(source_damage_anisoDuctile_Noutput(instance),instance) = damage_drivingforce_ID
+             source_damage_anisoDuctile_output(source_damage_anisoDuctile_Noutput(instance),instance) = &
+                                                       IO_lc(IO_stringValue(line,chunkPos,2_pInt))
+          end select
+
+       case ('anisoductile_atol')
+         source_damage_anisoDuctile_aTol(instance) = IO_floatValue(line,chunkPos,2_pInt)
+         
+       case ('nslip')  !
+         Nchunks_SlipFamilies = chunkPos(1) - 1_pInt
+         do j = 1_pInt, Nchunks_SlipFamilies
+           source_damage_anisoDuctile_Nslip(j,instance) = IO_intValue(line,chunkPos,1_pInt+j)
+         enddo
+
+       case ('anisoductile_sdot0')
+         source_damage_anisoDuctile_sdot_0(instance) = IO_floatValue(line,chunkPos,2_pInt)
+         
+       case ('anisoductile_criticalplasticstrain')
+         do j = 1_pInt, Nchunks_SlipFamilies
+           source_damage_anisoDuctile_critPlasticStrain(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j)
+         enddo
+         
+       case ('anisoductile_ratesensitivity')
+         source_damage_anisoDuctile_N(instance) = IO_floatValue(line,chunkPos,2_pInt)
+
+       case ('anisoductile_criticalload')
+         do j = 1_pInt, Nchunks_SlipFamilies
+           source_damage_anisoDuctile_critLoad(j,instance) = IO_floatValue(line,chunkPos,1_pInt+j)
+         enddo
+         
+     end select
+   endif; endif
+ enddo parsingFile
+
+!--------------------------------------------------------------------------------------------------
+!  sanity checks
+ sanityChecks: do phase = 1_pInt, size(phase_source)   
+   myPhase: if (any(phase_source(:,phase) == SOURCE_damage_anisoDuctile_ID)) then
+     instance = source_damage_anisoDuctile_instance(phase)
+     source_damage_anisoDuctile_Nslip(1:lattice_maxNslipFamily,instance) = &
+       min(lattice_NslipSystem(1:lattice_maxNslipFamily,phase),&                                    ! limit active cleavage systems per family to min of available and requested
+           source_damage_anisoDuctile_Nslip(1:lattice_maxNslipFamily,instance))
+         source_damage_anisoDuctile_totalNslip(instance) = sum(source_damage_anisoDuctile_Nslip(:,instance))
+     if (source_damage_anisoDuctile_aTol(instance) < 0.0_pReal) &
+       source_damage_anisoDuctile_aTol(instance) = 1.0e-3_pReal                                     ! default absolute tolerance 1e-3
+     if (source_damage_anisoDuctile_sdot_0(instance) <= 0.0_pReal) &
+       call IO_error(211_pInt,el=instance,ext_msg='sdot_0 ('//SOURCE_damage_anisoDuctile_LABEL//')')
+     if (any(source_damage_anisoDuctile_critPlasticStrain(:,instance) < 0.0_pReal)) &
+       call IO_error(211_pInt,el=instance,ext_msg='criticaPlasticStrain ('//SOURCE_damage_anisoDuctile_LABEL//')')
+     if (source_damage_anisoDuctile_N(instance) <= 0.0_pReal) &
+       call IO_error(211_pInt,el=instance,ext_msg='rate_sensitivity ('//SOURCE_damage_anisoDuctile_LABEL//')')
+   endif myPhase
+ enddo sanityChecks
+  
+
+ initializeInstances: do phase = 1_pInt, material_Nphase
+   if (any(phase_source(:,phase) == SOURCE_damage_anisoDuctile_ID)) then
+     NofMyPhase=count(material_phase==phase)
+     instance = source_damage_anisoDuctile_instance(phase)
+     sourceOffset = source_damage_anisoDuctile_offset(phase)
+
+!--------------------------------------------------------------------------------------------------
+!  Determine size of postResults array
+     outputsLoop: do o = 1_pInt,source_damage_anisoDuctile_Noutput(instance)
+       select case(source_damage_anisoDuctile_outputID(o,instance))
+         case(damage_drivingforce_ID)
+           mySize = 1_pInt
+       end select
+ 
+       if (mySize > 0_pInt) then  ! any meaningful output found
+          source_damage_anisoDuctile_sizePostResult(o,instance) = mySize
+          source_damage_anisoDuctile_sizePostResults(instance)  = source_damage_anisoDuctile_sizePostResults(instance) + mySize
+       endif
+     enddo outputsLoop
+
+!--------------------------------------------------------------------------------------------------
+! Determine size of state array
+     sizeDotState              = 1_pInt
+     sizeDeltaState            = 0_pInt
+     sizeState                 = 1_pInt
+     sourceState(phase)%p(sourceOffset)%sizeState      = sizeState
+     sourceState(phase)%p(sourceOffset)%sizeDotState   = sizeDotState
+     sourceState(phase)%p(sourceOffset)%sizeDeltaState = sizeDeltaState
+     sourceState(phase)%p(sourceOffset)%sizePostResults = source_damage_anisoDuctile_sizePostResults(instance)
+     allocate(sourceState(phase)%p(sourceOffset)%aTolState           (sizeState),                &
+              source=source_damage_anisoDuctile_aTol(instance))
+     allocate(sourceState(phase)%p(sourceOffset)%state0              (sizeState,NofMyPhase),     source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%partionedState0     (sizeState,NofMyPhase),     source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%subState0           (sizeState,NofMyPhase),     source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%state               (sizeState,NofMyPhase),     source=0.0_pReal)
+
+     allocate(sourceState(phase)%p(sourceOffset)%dotState            (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%deltaState        (sizeDeltaState,NofMyPhase),  source=0.0_pReal)
+     if (.not. analyticJaco) then
+       allocate(sourceState(phase)%p(sourceOffset)%state_backup      (sizeState,NofMyPhase),     source=0.0_pReal)
+       allocate(sourceState(phase)%p(sourceOffset)%dotState_backup   (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     endif
+     if (any(numerics_integrator == 1_pInt)) then
+       allocate(sourceState(phase)%p(sourceOffset)%previousDotState  (sizeDotState,NofMyPhase),  source=0.0_pReal)
+       allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     endif
+     if (any(numerics_integrator == 4_pInt)) &
+       allocate(sourceState(phase)%p(sourceOffset)%RK4dotState       (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     if (any(numerics_integrator == 5_pInt)) &
+       allocate(sourceState(phase)%p(sourceOffset)%RKCK45dotState    (6,sizeDotState,NofMyPhase),source=0.0_pReal)
+
+   endif
+ 
+ enddo initializeInstances
+end subroutine source_damage_anisoDuctile_init
+
+!--------------------------------------------------------------------------------------------------
+!> @brief calculates derived quantities from state
+!--------------------------------------------------------------------------------------------------
+subroutine source_damage_anisoDuctile_dotState(ipc, ip, el)
+ use material, only: &
+   phaseAt, phasememberAt, &
+   plasticState, &
+   sourceState, &
+   material_homog, &
+   damage, &
+   damageMapping
+ use lattice, only: &
+   lattice_maxNslipFamily
+
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< component-ID of integration point
+   ip, &                                                                                            !< integration point
+   el                                                                                               !< element
+ integer(pInt) :: &
+   phase, &
+   constituent, &
+   sourceOffset, &
+   homog, damageOffset, &
+   instance, &
+   index, f, i
+
+ phase = phaseAt(ipc,ip,el)
+ constituent = phasememberAt(ipc,ip,el)
+ instance = source_damage_anisoDuctile_instance(phase)
+ sourceOffset = source_damage_anisoDuctile_offset(phase)
+ homog = material_homog(ip,el)
+ damageOffset = damageMapping(homog)%p(ip,el)
+
+ index = 1_pInt
+ sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = 0.0_pReal
+ do f = 1_pInt,lattice_maxNslipFamily
+   do i = 1_pInt,source_damage_anisoDuctile_Nslip(f,instance)                                              ! process each (active) slip system in family
+     sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = &
+       sourceState(phase)%p(sourceOffset)%dotState(1,constituent) + &
+       plasticState(phase)%slipRate(index,constituent)/ &
+       ((damage(homog)%p(damageOffset))**source_damage_anisoDuctile_N(instance))/ & 
+       source_damage_anisoDuctile_critPlasticStrain(f,instance) 
+     
+     index = index + 1_pInt
+   enddo
+ enddo
+ 
+end subroutine source_damage_anisoDuctile_dotState
+
+!--------------------------------------------------------------------------------------------------
+!> @brief returns local part of nonlocal damage driving force
+!--------------------------------------------------------------------------------------------------
+subroutine source_damage_anisoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ipc, ip,  el)
+ use material, only: &
+   phaseAt, phasememberAt, &
+   sourceState
+
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< component-ID of integration point
+   ip, &                                                                                            !< integration point
+   el                                                                                               !< element
+ real(pReal),  intent(in) :: &
+   phi
+ real(pReal),  intent(out) :: &
+   localphiDot, &
+   dLocalphiDot_dPhi
+ integer(pInt) :: &
+   phase, constituent, sourceOffset
+
+ phase = phaseAt(ipc,ip,el)
+ constituent = phasememberAt(ipc,ip,el)
+ sourceOffset = source_damage_anisoDuctile_offset(phase)
+ 
+ localphiDot = 1.0_pReal - &
+               sourceState(phase)%p(sourceOffset)%state(1,constituent)* &
+               phi
+ 
+ dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent)
+ 
+end subroutine source_damage_anisoDuctile_getRateAndItsTangent
+ 
+!--------------------------------------------------------------------------------------------------
+!> @brief return array of local damage results
+!--------------------------------------------------------------------------------------------------
+function source_damage_anisoDuctile_postResults(ipc,ip,el)
+ use material, only: &
+   phaseAt, phasememberAt, &
+   sourceState
+
+ implicit none
+ integer(pInt),              intent(in) :: &
+   ipc, &                                                                                           !< component-ID of integration point
+   ip, &                                                                                            !< integration point
+   el                                                                                               !< element
+ real(pReal), dimension(source_damage_anisoDuctile_sizePostResults( &
+                          source_damage_anisoDuctile_instance(phaseAt(ipc,ip,el)))) :: &
+   source_damage_anisoDuctile_postResults
+
+ integer(pInt) :: &
+   instance, phase, constituent, sourceOffset, o, c
+   
+ phase = phaseAt(ipc,ip,el)
+ constituent = phasememberAt(ipc,ip,el)
+ instance = source_damage_anisoDuctile_instance(phase)
+ sourceOffset = source_damage_anisoDuctile_offset(phase)
+
+ c = 0_pInt
+ source_damage_anisoDuctile_postResults = 0.0_pReal
+
+ do o = 1_pInt,source_damage_anisoDuctile_Noutput(instance)
+    select case(source_damage_anisoDuctile_outputID(o,instance))
+      case (damage_drivingforce_ID)
+        source_damage_anisoDuctile_postResults(c+1_pInt) = &
+          sourceState(phase)%p(sourceOffset)%state(1,constituent)
+        c = c + 1_pInt
+
+    end select
+ enddo
+end function source_damage_anisoDuctile_postResults
+
+end module source_damage_anisoDuctile

code/source_damage_isoBrittle.f90

@@ -0,0 +1,383 @@
+!--------------------------------------------------------------------------------------------------
+! $Id$
+!--------------------------------------------------------------------------------------------------
+!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
+!> @author Luv Sharma, Max-Planck-Institut für Eisenforschung GmbH
+!> @brief material subroutine incoprorating isotropic brittle damage source mechanism
+!> @details to be done
+!--------------------------------------------------------------------------------------------------
+module source_damage_isoBrittle
+ use prec, only: &
+   pReal, &
+   pInt
+
+ implicit none
+ private
+ integer(pInt),                       dimension(:),           allocatable,         public, protected :: &
+   source_damage_isoBrittle_sizePostResults, &                                                        !< cumulative size of post results
+   source_damage_isoBrittle_offset, &                                                                 !< which source is my current damage mechanism?
+   source_damage_isoBrittle_instance                                                                  !< instance of damage source mechanism
+
+ integer(pInt),                       dimension(:,:),         allocatable, target, public :: &
+   source_damage_isoBrittle_sizePostResult                                                            !< size of each post result output
+
+ character(len=64),                   dimension(:,:),         allocatable, target, public :: &
+   source_damage_isoBrittle_output                                                                    !< name of each post result output
+   
+ integer(pInt),                       dimension(:),           allocatable, target, public :: &
+   source_damage_isoBrittle_Noutput                                                                   !< number of outputs per instance of this damage 
+
+ real(pReal),                         dimension(:),     allocatable,         private :: &
+   source_damage_isoBrittle_aTol, &
+   source_damage_isoBrittle_N, &
+   source_damage_isoBrittle_critStrainEnergy
+
+ enum, bind(c) 
+   enumerator :: undefined_ID, &
+                 damage_drivingforce_ID
+ end enum                                                 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11 ToDo
+ 
+ integer(kind(undefined_ID)),         dimension(:,:),         allocatable,          private :: & 
+   source_damage_isoBrittle_outputID                                                                  !< ID of each post result output
+
+
+ public :: &
+   source_damage_isoBrittle_init, &
+   source_damage_isoBrittle_deltaState, &
+   source_damage_isoBrittle_getRateAndItsTangent, &
+   source_damage_isoBrittle_postResults
+
+contains
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief module initialization
+!> @details reads in material parameters, allocates arrays, and does sanity checks
+!--------------------------------------------------------------------------------------------------
+subroutine source_damage_isoBrittle_init(fileUnit)
+ use, intrinsic :: iso_fortran_env                                                                  ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
+ use debug, only: &
+   debug_level,&
+   debug_constitutive,&
+   debug_levelBasic
+ use IO, only: &
+   IO_read, &
+   IO_lc, &
+   IO_getTag, &
+   IO_isBlank, &
+   IO_stringPos, &
+   IO_stringValue, &
+   IO_floatValue, &
+   IO_intValue, &
+   IO_warning, &
+   IO_error, &
+   IO_timeStamp, &
+   IO_EOF
+ use material, only: &
+   phase_source, &
+   phase_Nsources, &
+   phase_Noutput, &
+   SOURCE_damage_isoBrittle_label, &
+   SOURCE_damage_isoBrittle_ID, &
+   material_Nphase, &
+   material_phase, &  
+   sourceState, &
+   MATERIAL_partPhase
+ use numerics,only: &
+   analyticJaco, &
+   worldrank, &
+   numerics_integrator
+
+ implicit none
+ integer(pInt), intent(in) :: fileUnit
+
+ integer(pInt), allocatable, dimension(:) :: chunkPos
+ integer(pInt) :: maxNinstance,mySize=0_pInt,phase,instance,source,sourceOffset,o
+ integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
+ integer(pInt) :: NofMyPhase   
+ character(len=65536) :: &
+   tag  = '', &
+   line = ''
+
+ mainProcess: if (worldrank == 0) then 
+   write(6,'(/,a)')   ' <<<+-  source_'//SOURCE_damage_isoBrittle_label//' init  -+>>>'
+   write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
+#include "compilation_info.f90"
+ endif mainProcess
+
+ maxNinstance = int(count(phase_source == SOURCE_damage_isoBrittle_ID),pInt)
+ if (maxNinstance == 0_pInt) return
+ 
+ if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
+   write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
+ 
+ allocate(source_damage_isoBrittle_offset(material_Nphase), source=0_pInt)
+ allocate(source_damage_isoBrittle_instance(material_Nphase), source=0_pInt)
+ do phase = 1, material_Nphase
+   source_damage_isoBrittle_instance(phase) = count(phase_source(:,1:phase) == source_damage_isoBrittle_ID)
+   do source = 1, phase_Nsources(phase)
+     if (phase_source(source,phase) == source_damage_isoBrittle_ID) &
+       source_damage_isoBrittle_offset(phase) = source
+   enddo    
+ enddo
+   
+ allocate(source_damage_isoBrittle_sizePostResults(maxNinstance),                     source=0_pInt)
+ allocate(source_damage_isoBrittle_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
+ allocate(source_damage_isoBrittle_output(maxval(phase_Noutput),maxNinstance))
+          source_damage_isoBrittle_output = ''
+ allocate(source_damage_isoBrittle_outputID(maxval(phase_Noutput),maxNinstance),      source=undefined_ID)
+ allocate(source_damage_isoBrittle_Noutput(maxNinstance),                             source=0_pInt) 
+ allocate(source_damage_isoBrittle_critStrainEnergy(maxNinstance),                    source=0.0_pReal) 
+ allocate(source_damage_isoBrittle_N(maxNinstance),                                   source=1.0_pReal) 
+ allocate(source_damage_isoBrittle_aTol(maxNinstance),                                source=0.0_pReal) 
+
+ rewind(fileUnit)
+ phase = 0_pInt
+ do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase)         ! wind forward to <phase>
+   line = IO_read(fileUnit)
+ enddo
+ 
+ parsingFile: do while (trim(line) /= IO_EOF)                                                       ! read through sections of phase part
+   line = IO_read(fileUnit)
+   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
+   if (IO_getTag(line,'<','>') /= '') then                                                          ! stop at next part
+     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
+     exit                                                                                           
+   endif   
+   if (IO_getTag(line,'[',']') /= '') then                                                          ! next phase section
+     phase = phase + 1_pInt                                                                         ! advance phase section counter
+     cycle                                                                                          ! skip to next line
+   endif
+   if (phase > 0_pInt ) then; if (any(phase_source(:,phase) == SOURCE_damage_isoBrittle_ID)) then   ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
+     instance = source_damage_isoBrittle_instance(phase)                                            ! which instance of my damage is present phase
+     chunkPos = IO_stringPos(line)
+     tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt))                                             ! extract key
+     select case(tag)
+       case ('(output)')
+         select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt)))
+           case ('isobrittle_drivingforce')
+             source_damage_isoBrittle_Noutput(instance) = source_damage_isoBrittle_Noutput(instance) + 1_pInt
+             source_damage_isoBrittle_outputID(source_damage_isoBrittle_Noutput(instance),instance) = damage_drivingforce_ID
+             source_damage_isoBrittle_output(source_damage_isoBrittle_Noutput(instance),instance) = &
+                                                       IO_lc(IO_stringValue(line,chunkPos,2_pInt))
+          end select
+
+       case ('isobrittle_criticalstrainenergy')
+         source_damage_isoBrittle_critStrainEnergy(instance) = IO_floatValue(line,chunkPos,2_pInt)
+
+       case ('isobrittle_n')
+         source_damage_isoBrittle_N(instance) = IO_floatValue(line,chunkPos,2_pInt)
+
+       case ('isobrittle_atol')
+         source_damage_isoBrittle_aTol(instance) = IO_floatValue(line,chunkPos,2_pInt)
+
+     end select
+   endif; endif
+ enddo parsingFile
+
+
+!--------------------------------------------------------------------------------------------------
+!  sanity checks
+ sanityChecks: do phase = 1_pInt, material_Nphase 
+   myPhase: if (any(phase_source(:,phase) == SOURCE_damage_isoBrittle_ID)) then
+     instance = source_damage_isoBrittle_instance(phase)
+     if (source_damage_isoBrittle_aTol(instance) < 0.0_pReal) &
+       source_damage_isoBrittle_aTol(instance) = 1.0e-3_pReal                                              ! default absolute tolerance 1e-3
+     if (source_damage_isoBrittle_critStrainEnergy(instance) <= 0.0_pReal) &
+       call IO_error(211_pInt,el=instance,ext_msg='criticalStrainEnergy ('//SOURCE_damage_isoBrittle_LABEL//')')
+   endif myPhase
+ enddo sanityChecks
+
+ initializeInstances: do phase = 1_pInt, material_Nphase
+   if (any(phase_source(:,phase) == SOURCE_damage_isoBrittle_ID)) then
+     NofMyPhase=count(material_phase==phase)
+     instance = source_damage_isoBrittle_instance(phase)
+     sourceOffset = source_damage_isoBrittle_offset(phase)
+!--------------------------------------------------------------------------------------------------
+!  Determine size of postResults array
+     outputsLoop: do o = 1_pInt,source_damage_isoBrittle_Noutput(instance)
+       select case(source_damage_isoBrittle_outputID(o,instance))
+         case(damage_drivingforce_ID)
+           mySize = 1_pInt
+       end select
+ 
+       if (mySize > 0_pInt) then  ! any meaningful output found
+          source_damage_isoBrittle_sizePostResult(o,instance) = mySize
+          source_damage_isoBrittle_sizePostResults(instance)  = source_damage_isoBrittle_sizePostResults(instance) + mySize
+       endif
+     enddo outputsLoop
+! Determine size of state array
+     sizeDotState              =   1_pInt
+     sizeDeltaState            =   1_pInt
+     sizeState                 =   1_pInt
+                
+     sourceState(phase)%p(sourceOffset)%sizeState      = sizeState
+     sourceState(phase)%p(sourceOffset)%sizeDotState   = sizeDotState
+     sourceState(phase)%p(sourceOffset)%sizeDeltaState = sizeDeltaState
+     sourceState(phase)%p(sourceOffset)%sizePostResults = source_damage_isoBrittle_sizePostResults(instance)
+     allocate(sourceState(phase)%p(sourceOffset)%aTolState           (sizeState),                &
+              source=source_damage_isoBrittle_aTol(instance))
+     allocate(sourceState(phase)%p(sourceOffset)%state0              (sizeState,NofMyPhase),     source=.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%partionedState0     (sizeState,NofMyPhase),     source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%subState0           (sizeState,NofMyPhase),     source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%state               (sizeState,NofMyPhase),     source=0.0_pReal)
+
+     allocate(sourceState(phase)%p(sourceOffset)%dotState            (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%deltaState        (sizeDeltaState,NofMyPhase),  source=0.0_pReal)
+     if (.not. analyticJaco) then
+       allocate(sourceState(phase)%p(sourceOffset)%state_backup      (sizeState,NofMyPhase),     source=0.0_pReal)
+       allocate(sourceState(phase)%p(sourceOffset)%dotState_backup   (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     endif
+     if (any(numerics_integrator == 1_pInt)) then
+       allocate(sourceState(phase)%p(sourceOffset)%previousDotState  (sizeDotState,NofMyPhase),  source=0.0_pReal)
+       allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     endif
+     if (any(numerics_integrator == 4_pInt)) &
+       allocate(sourceState(phase)%p(sourceOffset)%RK4dotState       (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     if (any(numerics_integrator == 5_pInt)) &
+       allocate(sourceState(phase)%p(sourceOffset)%RKCK45dotState    (6,sizeDotState,NofMyPhase),source=0.0_pReal)
+
+   endif
+ 
+ enddo initializeInstances
+end subroutine source_damage_isoBrittle_init
+
+!--------------------------------------------------------------------------------------------------
+!> @brief calculates derived quantities from state
+!--------------------------------------------------------------------------------------------------
+subroutine source_damage_isoBrittle_deltaState(C, Fe, ipc, ip, el)
+ use material, only: &
+   phaseAt, phasememberAt, &
+   sourceState, &
+   material_homog, &
+   phase_NstiffnessDegradations, &
+   phase_stiffnessDegradation, &
+   porosity, &
+   porosityMapping, &  
+   STIFFNESS_DEGRADATION_porosity_ID
+ use math, only : &
+   math_mul33x33, &
+   math_mul66x6, &
+   math_Mandel33to6, &
+   math_transpose33, &
+   math_I3
+
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< component-ID of integration point
+   ip, &                                                                                            !< integration point
+   el                                                                                               !< element
+ real(pReal),  intent(in), dimension(3,3) :: &
+   Fe
+ real(pReal),  intent(in), dimension(6,6) :: &
+   C
+ integer(pInt) :: &
+   phase, constituent, instance, sourceOffset, mech
+ real(pReal) :: &
+   strain(6), &
+   stiffness(6,6), &
+   strainenergy
+
+ phase = phaseAt(ipc,ip,el)                                                            !< phase ID at ipc,ip,el
+ constituent = phasememberAt(ipc,ip,el)                                                      !< state array offset for phase ID at ipc,ip,el
+ ! ToDo: capability for multiple instances of SAME source within given phase. Needs Ninstance loop from here on!
+ instance = source_damage_isoBrittle_instance(phase)                                                 !< instance of damage_isoBrittle source
+ sourceOffset = source_damage_isoBrittle_offset(phase)
+
+ stiffness = C
+ do mech = 1_pInt, phase_NstiffnessDegradations(phase)
+   select case(phase_stiffnessDegradation(mech,phase))
+     case (STIFFNESS_DEGRADATION_porosity_ID)
+       stiffness = porosity(material_homog(ip,el))%p(porosityMapping(material_homog(ip,el))%p(ip,el))* &
+                   porosity(material_homog(ip,el))%p(porosityMapping(material_homog(ip,el))%p(ip,el))* &
+                   stiffness 
+   end select
+ enddo                                                 
+ strain = 0.5_pReal*math_Mandel33to6(math_mul33x33(math_transpose33(Fe),Fe)-math_I3)
+
+ strainenergy = 2.0_pReal*sum(strain*math_mul66x6(stiffness,strain))/ &
+                source_damage_isoBrittle_critStrainEnergy(instance) 
+ if (strainenergy > sourceState(phase)%p(sourceOffset)%subState0(1,constituent)) then
+   sourceState(phase)%p(sourceOffset)%deltaState(1,constituent) = &
+     strainenergy - sourceState(phase)%p(sourceOffset)%state(1,constituent)
+ else
+   sourceState(phase)%p(sourceOffset)%deltaState(1,constituent) = &
+     sourceState(phase)%p(sourceOffset)%subState0(1,constituent) - &
+     sourceState(phase)%p(sourceOffset)%state(1,constituent)
+ endif
+ 
+end subroutine source_damage_isoBrittle_deltaState
+ 
+!--------------------------------------------------------------------------------------------------
+!> @brief returns local part of nonlocal damage driving force
+!--------------------------------------------------------------------------------------------------
+subroutine source_damage_isoBrittle_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ipc, ip, el)
+ use material, only: &
+   phaseAt, phasememberAt, &
+   sourceState
+
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< component-ID of integration point
+   ip, &                                                                                            !< integration point
+   el                                                                                               !< element
+ real(pReal),  intent(in) :: &
+   phi
+ real(pReal),  intent(out) :: &
+   localphiDot, &
+   dLocalphiDot_dPhi
+ integer(pInt) :: &
+   phase, constituent, instance, sourceOffset
+
+ phase = phaseAt(ipc,ip,el)
+ constituent = phasememberAt(ipc,ip,el)
+ instance = source_damage_isoBrittle_instance(phase)
+ sourceOffset = source_damage_isoBrittle_offset(phase)
+ 
+ localphiDot = (1.0_pReal - phi)**(source_damage_isoBrittle_N(instance) - 1.0_pReal) - &
+               phi*sourceState(phase)%p(sourceOffset)%state(1,constituent)
+ dLocalphiDot_dPhi = - (source_damage_isoBrittle_N(instance) - 1.0_pReal)* &
+                       (1.0_pReal - phi)**max(0.0_pReal,source_damage_isoBrittle_N(instance) - 2.0_pReal) &
+                     - sourceState(phase)%p(sourceOffset)%state(1,constituent)
+ 
+end subroutine source_damage_isoBrittle_getRateAndItsTangent
+ 
+!--------------------------------------------------------------------------------------------------
+!> @brief return array of local damage results
+!--------------------------------------------------------------------------------------------------
+function source_damage_isoBrittle_postResults(ipc,ip,el)
+ use material, only: &
+   phaseAt, phasememberAt, &
+   sourceState
+
+ implicit none
+ integer(pInt),              intent(in) :: &
+   ipc, &                                                                                           !< component-ID of integration point
+   ip, &                                                                                            !< integration point
+   el                                                                                               !< element
+ real(pReal), dimension(source_damage_isoBrittle_sizePostResults( &
+                         source_damage_isoBrittle_instance(phaseAt(ipc,ip,el)))) :: &
+   source_damage_isoBrittle_postResults
+
+ integer(pInt) :: &
+   instance, phase, constituent, sourceOffset, o, c
+   
+ phase = phaseAt(ipc,ip,el)
+ constituent = phasememberAt(ipc,ip,el)
+ instance = source_damage_isoBrittle_instance(phase)
+ sourceOffset = source_damage_isoBrittle_offset(phase)
+
+ c = 0_pInt
+ source_damage_isoBrittle_postResults = 0.0_pReal
+
+ do o = 1_pInt,source_damage_isoBrittle_Noutput(instance)
+    select case(source_damage_isoBrittle_outputID(o,instance))
+      case (damage_drivingforce_ID)
+        source_damage_isoBrittle_postResults(c+1_pInt) = sourceState(phase)%p(sourceOffset)%state(1,constituent)
+        c = c + 1
+
+    end select
+ enddo
+end function source_damage_isoBrittle_postResults
+
+end module source_damage_isoBrittle

code/source_damage_isoDuctile.f90

@@ -0,0 +1,350 @@
+!--------------------------------------------------------------------------------------------------
+! $Id$
+!--------------------------------------------------------------------------------------------------
+!> @author Pratheek Shanthraj, Max-Planck-Institut fŸr Eisenforschung GmbH
+!> @author Luv Sharma, Max-Planck-Institut fŸr Eisenforschung GmbH
+!> @brief material subroutine incoprorating isotropic ductile damage source mechanism
+!> @details to be done
+!--------------------------------------------------------------------------------------------------
+module source_damage_isoDuctile
+ use prec, only: &
+   pReal, &
+   pInt
+
+ implicit none
+ private
+ integer(pInt),                       dimension(:),           allocatable,         public, protected :: &
+   source_damage_isoDuctile_sizePostResults, &                                                        !< cumulative size of post results
+   source_damage_isoDuctile_offset, &                                                                 !< which source is my current damage mechanism?
+   source_damage_isoDuctile_instance                                                                  !< instance of damage source mechanism
+
+ integer(pInt),                       dimension(:,:),         allocatable, target, public :: &
+   source_damage_isoDuctile_sizePostResult                                                            !< size of each post result output
+
+ character(len=64),                   dimension(:,:),         allocatable, target, public :: &
+   source_damage_isoDuctile_output                                                                    !< name of each post result output
+   
+ integer(pInt),                       dimension(:),           allocatable, target, public :: &
+   source_damage_isoDuctile_Noutput                                                                   !< number of outputs per instance of this damage 
+
+ real(pReal),                         dimension(:),     allocatable,         private :: &
+   source_damage_isoDuctile_aTol, &
+   source_damage_isoDuctile_critPlasticStrain, &
+   source_damage_isoDuctile_N
+
+ enum, bind(c) 
+   enumerator :: undefined_ID, &
+                 damage_drivingforce_ID
+ end enum                                                 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11 ToDo
+ 
+ integer(kind(undefined_ID)),         dimension(:,:),         allocatable,          private :: & 
+   source_damage_isoDuctile_outputID                                                                  !< ID of each post result output
+
+
+ public :: &
+   source_damage_isoDuctile_init, &
+   source_damage_isoDuctile_dotState, &
+   source_damage_isoDuctile_getRateAndItsTangent, &
+   source_damage_isoDuctile_postResults
+
+contains
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief module initialization
+!> @details reads in material parameters, allocates arrays, and does sanity checks
+!--------------------------------------------------------------------------------------------------
+subroutine source_damage_isoDuctile_init(fileUnit)
+ use, intrinsic :: iso_fortran_env                                                                  ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
+ use debug, only: &
+   debug_level,&
+   debug_constitutive,&
+   debug_levelBasic
+ use IO, only: &
+   IO_read, &
+   IO_lc, &
+   IO_getTag, &
+   IO_isBlank, &
+   IO_stringPos, &
+   IO_stringValue, &
+   IO_floatValue, &
+   IO_intValue, &
+   IO_warning, &
+   IO_error, &
+   IO_timeStamp, &
+   IO_EOF
+ use material, only: &
+   phase_source, &
+   phase_Nsources, &
+   phase_Noutput, &
+   SOURCE_damage_isoDuctile_label, &
+   SOURCE_damage_isoDuctile_ID, &
+   material_Nphase, &
+   material_phase, &  
+   sourceState, &
+   MATERIAL_partPhase
+ use numerics,only: &
+   analyticJaco, &
+   worldrank, &
+   numerics_integrator
+
+ implicit none
+ integer(pInt), intent(in) :: fileUnit
+
+ integer(pInt), allocatable, dimension(:) :: chunkPos
+ integer(pInt) :: maxNinstance,mySize=0_pInt,phase,instance,source,sourceOffset,o
+ integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
+ integer(pInt) :: NofMyPhase   
+ character(len=65536) :: &
+   tag  = '', &
+   line = ''
+
+ mainProcess: if (worldrank == 0) then 
+   write(6,'(/,a)')   ' <<<+-  source_'//SOURCE_damage_isoDuctile_label//' init  -+>>>'
+   write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
+#include "compilation_info.f90"
+ endif mainProcess
+
+ maxNinstance = int(count(phase_source == SOURCE_damage_isoDuctile_ID),pInt)
+ if (maxNinstance == 0_pInt) return
+ 
+ if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
+   write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
+ 
+ allocate(source_damage_isoDuctile_offset(material_Nphase), source=0_pInt)
+ allocate(source_damage_isoDuctile_instance(material_Nphase), source=0_pInt)
+ do phase = 1, material_Nphase
+   source_damage_isoDuctile_instance(phase) = count(phase_source(:,1:phase) == source_damage_isoDuctile_ID)
+   do source = 1, phase_Nsources(phase)
+     if (phase_source(source,phase) == source_damage_isoDuctile_ID) &
+       source_damage_isoDuctile_offset(phase) = source
+   enddo    
+ enddo
+   
+ allocate(source_damage_isoDuctile_sizePostResults(maxNinstance),                     source=0_pInt)
+ allocate(source_damage_isoDuctile_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
+ allocate(source_damage_isoDuctile_output(maxval(phase_Noutput),maxNinstance))
+          source_damage_isoDuctile_output = ''
+ allocate(source_damage_isoDuctile_outputID(maxval(phase_Noutput),maxNinstance),      source=undefined_ID)
+ allocate(source_damage_isoDuctile_Noutput(maxNinstance),                             source=0_pInt) 
+ allocate(source_damage_isoDuctile_critPlasticStrain(maxNinstance),                     source=0.0_pReal) 
+ allocate(source_damage_isoDuctile_N(maxNinstance),                                     source=0.0_pReal) 
+ allocate(source_damage_isoDuctile_aTol(maxNinstance),                                source=0.0_pReal) 
+
+ rewind(fileUnit)
+ phase = 0_pInt
+ do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase)         ! wind forward to <phase>
+   line = IO_read(fileUnit)
+ enddo
+ 
+ parsingFile: do while (trim(line) /= IO_EOF)                                                       ! read through sections of phase part
+   line = IO_read(fileUnit)
+   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
+   if (IO_getTag(line,'<','>') /= '') then                                                          ! stop at next part
+     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
+     exit                                                                                           
+   endif   
+   if (IO_getTag(line,'[',']') /= '') then                                                          ! next phase section
+     phase = phase + 1_pInt                                                                         ! advance phase section counter
+     cycle                                                                                          ! skip to next line
+   endif
+   if (phase > 0_pInt ) then; if (any(phase_source(:,phase) == SOURCE_damage_isoDuctile_ID)) then   ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
+     instance = source_damage_isoDuctile_instance(phase)                                            ! which instance of my damage is present phase
+     chunkPos = IO_stringPos(line)
+     tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt))                                             ! extract key
+     select case(tag)
+       case ('(output)')
+         select case(IO_lc(IO_stringValue(line,chunkPos,2_pInt)))
+           case ('isoductile_drivingforce')
+             source_damage_isoDuctile_Noutput(instance) = source_damage_isoDuctile_Noutput(instance) + 1_pInt
+             source_damage_isoDuctile_outputID(source_damage_isoDuctile_Noutput(instance),instance) = damage_drivingforce_ID
+             source_damage_isoDuctile_output(source_damage_isoDuctile_Noutput(instance),instance) = &
+                                                       IO_lc(IO_stringValue(line,chunkPos,2_pInt))
+          end select
+
+       case ('isoductile_criticalplasticstrain')
+         source_damage_isoDuctile_critPlasticStrain(instance) = IO_floatValue(line,chunkPos,2_pInt)
+
+       case ('isoductile_ratesensitivity')
+         source_damage_isoDuctile_N(instance) = IO_floatValue(line,chunkPos,2_pInt)
+
+       case ('isoductile_atol')
+         source_damage_isoDuctile_aTol(instance) = IO_floatValue(line,chunkPos,2_pInt)
+
+     end select
+   endif; endif
+ enddo parsingFile
+
+
+!--------------------------------------------------------------------------------------------------
+!  sanity checks
+ sanityChecks: do phase = 1_pInt, material_Nphase 
+   myPhase: if (any(phase_source(:,phase) == SOURCE_damage_isoDuctile_ID)) then
+     instance = source_damage_isoDuctile_instance(phase)
+     if (source_damage_isoDuctile_aTol(instance) < 0.0_pReal) &
+       source_damage_isoDuctile_aTol(instance) = 1.0e-3_pReal                                              ! default absolute tolerance 1e-3
+     if (source_damage_isoDuctile_critPlasticStrain(instance) <= 0.0_pReal) &
+       call IO_error(211_pInt,el=instance,ext_msg='critical plastic strain ('//SOURCE_damage_isoDuctile_LABEL//')')
+   endif myPhase
+ enddo sanityChecks
+
+ initializeInstances: do phase = 1_pInt, material_Nphase
+   if (any(phase_source(:,phase) == SOURCE_damage_isoDuctile_ID)) then
+     NofMyPhase=count(material_phase==phase)
+     instance = source_damage_isoDuctile_instance(phase)
+     sourceOffset = source_damage_isoDuctile_offset(phase)
+!--------------------------------------------------------------------------------------------------
+!  Determine size of postResults array
+     outputsLoop: do o = 1_pInt,source_damage_isoDuctile_Noutput(instance)
+       select case(source_damage_isoDuctile_outputID(o,instance))
+         case(damage_drivingforce_ID)
+           mySize = 1_pInt
+       end select
+ 
+       if (mySize > 0_pInt) then  ! any meaningful output found
+          source_damage_isoDuctile_sizePostResult(o,instance) = mySize
+          source_damage_isoDuctile_sizePostResults(instance)  = source_damage_isoDuctile_sizePostResults(instance) + mySize
+       endif
+     enddo outputsLoop
+! Determine size of state array
+     sizeDotState              =   1_pInt
+     sizeDeltaState            =   0_pInt
+     sizeState                 =   1_pInt
+                
+     sourceState(phase)%p(sourceOffset)%sizeState      = sizeState
+     sourceState(phase)%p(sourceOffset)%sizeDotState   = sizeDotState
+     sourceState(phase)%p(sourceOffset)%sizeDeltaState = sizeDeltaState
+     sourceState(phase)%p(sourceOffset)%sizePostResults = source_damage_isoDuctile_sizePostResults(instance)
+     allocate(sourceState(phase)%p(sourceOffset)%aTolState           (sizeState),                &
+              source=source_damage_isoDuctile_aTol(instance))
+     allocate(sourceState(phase)%p(sourceOffset)%state0              (sizeState,NofMyPhase),     source=.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%partionedState0     (sizeState,NofMyPhase),     source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%subState0           (sizeState,NofMyPhase),     source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%state               (sizeState,NofMyPhase),     source=0.0_pReal)
+
+     allocate(sourceState(phase)%p(sourceOffset)%dotState            (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%deltaState        (sizeDeltaState,NofMyPhase),  source=0.0_pReal)
+     if (.not. analyticJaco) then
+       allocate(sourceState(phase)%p(sourceOffset)%state_backup      (sizeState,NofMyPhase),     source=0.0_pReal)
+       allocate(sourceState(phase)%p(sourceOffset)%dotState_backup   (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     endif
+     if (any(numerics_integrator == 1_pInt)) then
+       allocate(sourceState(phase)%p(sourceOffset)%previousDotState  (sizeDotState,NofMyPhase),  source=0.0_pReal)
+       allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     endif
+     if (any(numerics_integrator == 4_pInt)) &
+       allocate(sourceState(phase)%p(sourceOffset)%RK4dotState       (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     if (any(numerics_integrator == 5_pInt)) &
+       allocate(sourceState(phase)%p(sourceOffset)%RKCK45dotState    (6,sizeDotState,NofMyPhase),source=0.0_pReal)
+
+   endif
+ 
+ enddo initializeInstances
+end subroutine source_damage_isoDuctile_init
+
+!--------------------------------------------------------------------------------------------------
+!> @brief calculates derived quantities from state
+!--------------------------------------------------------------------------------------------------
+subroutine source_damage_isoDuctile_dotState(ipc, ip, el)
+ use material, only: &
+   phaseAt, phasememberAt, &
+   plasticState, &
+   sourceState, &
+   material_homog, &
+   damage, &
+   damageMapping
+
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< component-ID of integration point
+   ip, &                                                                                            !< integration point
+   el                                                                                               !< element
+ integer(pInt) :: &
+   phase, constituent, instance, homog, sourceOffset, damageOffset
+
+ phase = phaseAt(ipc,ip,el)
+ constituent = phasememberAt(ipc,ip,el)
+ instance = source_damage_isoDuctile_instance(phase)
+ sourceOffset = source_damage_isoDuctile_offset(phase)
+ homog = material_homog(ip,el)
+ damageOffset = damageMapping(homog)%p(ip,el)
+
+ sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = &
+   sum(plasticState(phase)%slipRate(:,constituent))/ &
+   ((damage(homog)%p(damageOffset))**source_damage_isoDuctile_N(instance))/ & 
+   source_damage_isoDuctile_critPlasticStrain(instance) 
+
+end subroutine source_damage_isoDuctile_dotState
+ 
+!--------------------------------------------------------------------------------------------------
+!> @brief returns local part of nonlocal damage driving force
+!--------------------------------------------------------------------------------------------------
+subroutine source_damage_isoDuctile_getRateAndItsTangent(localphiDot, dLocalphiDot_dPhi, phi, ipc, ip, el)
+ use material, only: &
+   phaseAt, phasememberAt, &
+   sourceState
+
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< component-ID of integration point
+   ip, &                                                                                            !< integration point
+   el                                                                                               !< element
+ real(pReal),  intent(in) :: &
+   phi
+ real(pReal),  intent(out) :: &
+   localphiDot, &
+   dLocalphiDot_dPhi
+ integer(pInt) :: &
+   phase, constituent, sourceOffset
+
+ phase = phaseAt(ipc,ip,el)
+ constituent = phasememberAt(ipc,ip,el)
+ sourceOffset = source_damage_isoDuctile_offset(phase)
+ 
+ localphiDot = 1.0_pReal - &
+               sourceState(phase)%p(sourceOffset)%state(1,constituent)* &
+               phi
+ 
+ dLocalphiDot_dPhi = -sourceState(phase)%p(sourceOffset)%state(1,constituent)
+ 
+end subroutine source_damage_isoDuctile_getRateAndItsTangent
+ 
+!--------------------------------------------------------------------------------------------------
+!> @brief return array of local damage results
+!--------------------------------------------------------------------------------------------------
+function source_damage_isoDuctile_postResults(ipc,ip,el)
+ use material, only: &
+   phaseAt, phasememberAt, &
+   sourceState
+
+ implicit none
+ integer(pInt),              intent(in) :: &
+   ipc, &                                                                                           !< component-ID of integration point
+   ip, &                                                                                            !< integration point
+   el                                                                                               !< element
+ real(pReal), dimension(source_damage_isoDuctile_sizePostResults( &
+                          source_damage_isoDuctile_instance(phaseAt(ipc,ip,el)))) :: &
+   source_damage_isoDuctile_postResults
+
+ integer(pInt) :: &
+   instance, phase, constituent, sourceOffset, o, c
+   
+ phase = phaseAt(ipc,ip,el)
+ constituent = phasememberAt(ipc,ip,el)
+ instance = source_damage_isoDuctile_instance(phase)
+ sourceOffset = source_damage_isoDuctile_offset(phase)
+
+ c = 0_pInt
+ source_damage_isoDuctile_postResults = 0.0_pReal
+
+ do o = 1_pInt,source_damage_isoDuctile_Noutput(instance)
+    select case(source_damage_isoDuctile_outputID(o,instance))
+      case (damage_drivingforce_ID)
+        source_damage_isoDuctile_postResults(c+1_pInt) = sourceState(phase)%p(sourceOffset)%state(1,constituent)
+        c = c + 1
+
+    end select
+ enddo
+end function source_damage_isoDuctile_postResults
+
+end module source_damage_isoDuctile

code/source_thermal_dissipation.f90

@@ -0,0 +1,220 @@
+!--------------------------------------------------------------------------------------------------
+! $Id$
+!--------------------------------------------------------------------------------------------------
+!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
+!> @brief material subroutine for thermal source due to plastic dissipation
+!> @details to be done
+!--------------------------------------------------------------------------------------------------
+module source_thermal_dissipation
+ use prec, only: &
+   pReal, &
+   pInt
+
+ implicit none
+ private
+ integer(pInt),                       dimension(:),           allocatable,         public, protected :: &
+   source_thermal_dissipation_sizePostResults, &                                                        !< cumulative size of post results
+   source_thermal_dissipation_offset, &                                                                 !< which source is my current thermal dissipation mechanism?
+   source_thermal_dissipation_instance                                                                  !< instance of thermal dissipation source mechanism
+
+ integer(pInt),                       dimension(:,:),         allocatable, target, public :: &
+   source_thermal_dissipation_sizePostResult                                                            !< size of each post result output
+
+ character(len=64),                   dimension(:,:),         allocatable, target, public :: &
+   source_thermal_dissipation_output                                                                    !< name of each post result output
+   
+ integer(pInt),                       dimension(:),           allocatable, target, public :: &
+   source_thermal_dissipation_Noutput                                                                   !< number of outputs per instance of this source 
+
+ real(pReal),                         dimension(:),           allocatable,        private :: &
+   source_thermal_dissipation_coldworkCoeff
+
+ public :: &
+   source_thermal_dissipation_init, &
+   source_thermal_dissipation_getRateAndItsTangent
+
+contains
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief module initialization
+!> @details reads in material parameters, allocates arrays, and does sanity checks
+!--------------------------------------------------------------------------------------------------
+subroutine source_thermal_dissipation_init(fileUnit)
+ use, intrinsic :: iso_fortran_env                                                                  ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
+ use debug, only: &
+   debug_level,&
+   debug_constitutive,&
+   debug_levelBasic
+ use IO, only: &
+   IO_read, &
+   IO_lc, &
+   IO_getTag, &
+   IO_isBlank, &
+   IO_stringPos, &
+   IO_stringValue, &
+   IO_floatValue, &
+   IO_intValue, &
+   IO_warning, &
+   IO_error, &
+   IO_timeStamp, &
+   IO_EOF
+ use material, only: &
+   phase_source, &
+   phase_Nsources, &
+   phase_Noutput, &
+   SOURCE_thermal_dissipation_label, &
+   SOURCE_thermal_dissipation_ID, &
+   material_Nphase, &
+   material_phase, &  
+   sourceState, &
+   MATERIAL_partPhase
+ use numerics,only: &
+   analyticJaco, &
+   worldrank, &
+   numerics_integrator
+
+ implicit none
+ integer(pInt), intent(in) :: fileUnit
+
+ integer(pInt), allocatable, dimension(:) :: chunkPos
+ integer(pInt) :: maxNinstance,phase,instance,source,sourceOffset
+ integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
+ integer(pInt) :: NofMyPhase   
+ character(len=65536) :: &
+   tag  = '', &
+   line = ''
+
+ mainProcess: if (worldrank == 0) then 
+   write(6,'(/,a)')   ' <<<+-  source_'//SOURCE_thermal_dissipation_label//' init  -+>>>'
+   write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
+#include "compilation_info.f90"
+ endif mainProcess
+ 
+ maxNinstance = int(count(phase_source == SOURCE_thermal_dissipation_ID),pInt)
+ if (maxNinstance == 0_pInt) return
+ if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
+   write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
+ 
+ allocate(source_thermal_dissipation_offset(material_Nphase), source=0_pInt)
+ allocate(source_thermal_dissipation_instance(material_Nphase), source=0_pInt)
+ do phase = 1, material_Nphase
+   source_thermal_dissipation_instance(phase) = count(phase_source(:,1:phase) == SOURCE_thermal_dissipation_ID)
+   do source = 1, phase_Nsources(phase)
+     if (phase_source(source,phase) == SOURCE_thermal_dissipation_ID) &
+       source_thermal_dissipation_offset(phase) = source
+   enddo    
+ enddo
+   
+ allocate(source_thermal_dissipation_sizePostResults(maxNinstance),                     source=0_pInt)
+ allocate(source_thermal_dissipation_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
+ allocate(source_thermal_dissipation_output  (maxval(phase_Noutput),maxNinstance))
+          source_thermal_dissipation_output = ''
+ allocate(source_thermal_dissipation_Noutput(maxNinstance),                             source=0_pInt) 
+ allocate(source_thermal_dissipation_coldworkCoeff(maxNinstance),                       source=0.0_pReal) 
+
+ rewind(fileUnit)
+ phase = 0_pInt
+ do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase)         ! wind forward to <phase>
+   line = IO_read(fileUnit)
+ enddo
+ 
+ parsingFile: do while (trim(line) /= IO_EOF)                                                       ! read through sections of phase part
+   line = IO_read(fileUnit)
+   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
+   if (IO_getTag(line,'<','>') /= '') then                                                          ! stop at next part
+     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
+     exit                                                                                           
+   endif   
+   if (IO_getTag(line,'[',']') /= '') then                                                          ! next phase section
+     phase = phase + 1_pInt                                                                         ! advance phase section counter
+     cycle                                                                                          ! skip to next line
+   endif
+
+   if (phase > 0_pInt ) then; if (any(phase_source(:,phase) == SOURCE_thermal_dissipation_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
+
+     instance = source_thermal_dissipation_instance(phase)                                          ! which instance of my source is present phase
+     chunkPos = IO_stringPos(line)
+     tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt))                                             ! extract key
+     select case(tag)
+       case ('dissipation_coldworkcoeff')
+         source_thermal_dissipation_coldworkCoeff(instance) = IO_floatValue(line,chunkPos,2_pInt)
+
+     end select
+   endif; endif
+ enddo parsingFile
+ 
+ initializeInstances: do phase = 1_pInt, material_Nphase
+   if (any(phase_source(:,phase) == SOURCE_thermal_dissipation_ID)) then
+     NofMyPhase=count(material_phase==phase)
+     instance = source_thermal_dissipation_instance(phase)
+     sourceOffset = source_thermal_dissipation_offset(phase)
+
+     sizeDotState              =   0_pInt
+     sizeDeltaState            =   0_pInt
+     sizeState                 =   0_pInt
+     sourceState(phase)%p(sourceOffset)%sizeState =       sizeState
+     sourceState(phase)%p(sourceOffset)%sizeDotState =    sizeDotState
+     sourceState(phase)%p(sourceOffset)%sizeDeltaState =  sizeDeltaState
+     sourceState(phase)%p(sourceOffset)%sizePostResults = source_thermal_dissipation_sizePostResults(instance)
+     allocate(sourceState(phase)%p(sourceOffset)%aTolState           (sizeState),                source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%state0              (sizeState,NofMyPhase),     source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%partionedState0     (sizeState,NofMyPhase),     source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%subState0           (sizeState,NofMyPhase),     source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%state               (sizeState,NofMyPhase),     source=0.0_pReal)
+
+     allocate(sourceState(phase)%p(sourceOffset)%dotState            (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%deltaState        (sizeDeltaState,NofMyPhase),  source=0.0_pReal)
+     if (.not. analyticJaco) then
+       allocate(sourceState(phase)%p(sourceOffset)%state_backup      (sizeState,NofMyPhase),     source=0.0_pReal)
+       allocate(sourceState(phase)%p(sourceOffset)%dotState_backup   (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     endif
+     if (any(numerics_integrator == 1_pInt)) then
+       allocate(sourceState(phase)%p(sourceOffset)%previousDotState  (sizeDotState,NofMyPhase),  source=0.0_pReal)
+       allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     endif
+     if (any(numerics_integrator == 4_pInt)) &
+       allocate(sourceState(phase)%p(sourceOffset)%RK4dotState       (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     if (any(numerics_integrator == 5_pInt)) &
+       allocate(sourceState(phase)%p(sourceOffset)%RKCK45dotState    (6,sizeDotState,NofMyPhase),source=0.0_pReal)      
+
+   endif
+ 
+ enddo initializeInstances
+end subroutine source_thermal_dissipation_init
+
+!--------------------------------------------------------------------------------------------------
+!> @brief returns local vacancy generation rate 
+!--------------------------------------------------------------------------------------------------
+subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDOT_dT, Tstar_v, Lp, ipc, ip, el)
+ use math, only: &
+   math_Mandel6to33
+ use material, only: &
+   phaseAt, phasememberAt
+
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< grain number
+   ip, &                                                                                            !< integration point number
+   el                                                                                               !< element number
+ real(pReal),  intent(in), dimension(6) :: &
+   Tstar_v                                                                                          !< 2nd Piola Kirchhoff stress tensor (Mandel)
+ real(pReal),  intent(in), dimension(3,3) :: &
+   Lp
+ real(pReal),  intent(out) :: &
+   TDot, &
+   dTDOT_dT
+ integer(pInt) :: &
+   instance, phase, constituent 
+
+ phase = phaseAt(ipc,ip,el)
+ constituent = phasememberAt(ipc,ip,el)
+ instance = source_thermal_dissipation_instance(phase)
+ 
+ TDot = source_thermal_dissipation_coldworkCoeff(instance)* &
+        sum(abs(math_Mandel6to33(Tstar_v)*Lp))
+ dTDOT_dT = 0.0_pReal       
+ 
+end subroutine source_thermal_dissipation_getRateAndItsTangent
+ 
+end module source_thermal_dissipation

code/source_thermal_externalheat.f90

@@ -0,0 +1,277 @@
+!--------------------------------------------------------------------------------------------------
+! $Id$
+!--------------------------------------------------------------------------------------------------
+!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
+!> @brief material subroutine for thermal source due to plastic dissipation
+!> @details to be done
+!--------------------------------------------------------------------------------------------------
+module source_thermal_externalheat
+ use prec, only: &
+   pReal, &
+   pInt
+
+ implicit none
+ private
+ integer(pInt),                       dimension(:),           allocatable,         public, protected :: &
+   source_thermal_externalheat_sizePostResults, &                                                        !< cumulative size of post results
+   source_thermal_externalheat_offset, &                                                                 !< which source is my current thermal dissipation mechanism?
+   source_thermal_externalheat_instance                                                                  !< instance of thermal dissipation source mechanism
+
+ integer(pInt),                       dimension(:,:),         allocatable, target, public :: &
+   source_thermal_externalheat_sizePostResult                                                            !< size of each post result output
+
+ character(len=64),                   dimension(:,:),         allocatable, target, public :: &
+   source_thermal_externalheat_output                                                                    !< name of each post result output
+   
+ integer(pInt),                       dimension(:),           allocatable, target, public :: &
+   source_thermal_externalheat_Noutput                                                                   !< number of outputs per instance of this source 
+
+ integer(pInt),                       dimension(:),           allocatable,        private :: &
+   source_thermal_externalheat_nIntervals
+
+ real(pReal),                         dimension(:,:),         allocatable,        private :: &
+   source_thermal_externalheat_time, &
+   source_thermal_externalheat_rate
+
+ public :: &
+   source_thermal_externalheat_init, &
+   source_thermal_externalheat_dotState, &
+   source_thermal_externalheat_getRateAndItsTangent
+
+contains
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief module initialization
+!> @details reads in material parameters, allocates arrays, and does sanity checks
+!--------------------------------------------------------------------------------------------------
+subroutine source_thermal_externalheat_init(fileUnit)
+ use, intrinsic :: iso_fortran_env                                                                  ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
+ use debug, only: &
+   debug_level,&
+   debug_constitutive,&
+   debug_levelBasic
+ use IO, only: &
+   IO_read, &
+   IO_lc, &
+   IO_getTag, &
+   IO_isBlank, &
+   IO_stringPos, &
+   IO_stringValue, &
+   IO_floatValue, &
+   IO_intValue, &
+   IO_warning, &
+   IO_error, &
+   IO_timeStamp, &
+   IO_EOF
+ use material, only: &
+   phase_source, &
+   phase_Nsources, &
+   phase_Noutput, &
+   SOURCE_thermal_externalheat_label, &
+   SOURCE_thermal_externalheat_ID, &
+   material_Nphase, &
+   material_phase, &  
+   sourceState, &
+   MATERIAL_partPhase
+ use numerics,only: &
+   analyticJaco, &
+   worldrank, &
+   numerics_integrator
+
+ implicit none
+ integer(pInt), intent(in) :: fileUnit
+
+ integer(pInt), allocatable, dimension(:) :: chunkPos
+ integer(pInt) :: maxNinstance,phase,instance,source,sourceOffset
+ integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
+ integer(pInt) :: NofMyPhase,interval   
+ character(len=65536) :: &
+   tag  = '', &
+   line = ''
+ real(pReal), allocatable, dimension(:,:) :: temp_time, temp_rate  
+
+ mainProcess: if (worldrank == 0) then 
+   write(6,'(/,a)')   ' <<<+-  source_'//SOURCE_thermal_externalheat_label//' init  -+>>>'
+   write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
+#include "compilation_info.f90"
+ endif mainProcess
+ 
+ maxNinstance = int(count(phase_source == SOURCE_thermal_externalheat_ID),pInt)
+ if (maxNinstance == 0_pInt) return
+ if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
+   write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
+ 
+ allocate(source_thermal_externalheat_offset(material_Nphase), source=0_pInt)
+ allocate(source_thermal_externalheat_instance(material_Nphase), source=0_pInt)
+ do phase = 1, material_Nphase
+   source_thermal_externalheat_instance(phase) = count(phase_source(:,1:phase) == SOURCE_thermal_externalheat_ID)
+   do source = 1, phase_Nsources(phase)
+     if (phase_source(source,phase) == SOURCE_thermal_externalheat_ID) &
+       source_thermal_externalheat_offset(phase) = source
+   enddo    
+ enddo
+   
+ allocate(source_thermal_externalheat_sizePostResults(maxNinstance),                     source=0_pInt)
+ allocate(source_thermal_externalheat_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
+ allocate(source_thermal_externalheat_output  (maxval(phase_Noutput),maxNinstance))
+          source_thermal_externalheat_output = ''
+ allocate(source_thermal_externalheat_Noutput(maxNinstance),                             source=0_pInt) 
+ allocate(source_thermal_externalheat_nIntervals(maxNinstance),                          source=0_pInt) 
+
+ allocate(temp_time(maxNinstance,1000), source=0.0_pReal) 
+ allocate(temp_rate(maxNinstance,1000), source=0.0_pReal) 
+
+ rewind(fileUnit)
+ phase = 0_pInt
+ do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase)         ! wind forward to <phase>
+   line = IO_read(fileUnit)
+ enddo
+ 
+ parsingFile: do while (trim(line) /= IO_EOF)                                                       ! read through sections of phase part
+   line = IO_read(fileUnit)
+   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
+   if (IO_getTag(line,'<','>') /= '') then                                                          ! stop at next part
+     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
+     exit                                                                                           
+   endif   
+   if (IO_getTag(line,'[',']') /= '') then                                                          ! next phase section
+     phase = phase + 1_pInt                                                                         ! advance phase section counter
+     cycle                                                                                          ! skip to next line
+   endif
+
+   if (phase > 0_pInt ) then; if (any(phase_source(:,phase) == SOURCE_thermal_externalheat_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
+
+     instance = source_thermal_externalheat_instance(phase)                                          ! which instance of my source is present phase
+     chunkPos = IO_stringPos(line)
+     tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt))                                             ! extract key
+     select case(tag)
+       case ('externalheat_time')
+         if (chunkPos(1) <= 2_pInt) &
+           call IO_error(150_pInt,ext_msg=trim(tag)//' ('//SOURCE_thermal_externalheat_label//')')
+         source_thermal_externalheat_nIntervals(instance) = chunkPos(1) - 2_pInt
+         do interval = 1, source_thermal_externalheat_nIntervals(instance) + 1_pInt
+           temp_time(instance, interval) = IO_floatValue(line,chunkPos,1_pInt + interval)
+         enddo  
+
+       case ('externalheat_rate')
+         do interval = 1, source_thermal_externalheat_nIntervals(instance) + 1_pInt
+           temp_rate(instance, interval) = IO_floatValue(line,chunkPos,1_pInt + interval)
+         enddo  
+
+     end select
+   endif; endif
+ enddo parsingFile
+ 
+ allocate(source_thermal_externalheat_time(maxNinstance,maxval(source_thermal_externalheat_nIntervals)+1_pInt), source=0.0_pReal) 
+ allocate(source_thermal_externalheat_rate(maxNinstance,maxval(source_thermal_externalheat_nIntervals)+1_pInt), source=0.0_pReal) 
+
+ initializeInstances: do phase = 1_pInt, material_Nphase
+   if (any(phase_source(:,phase) == SOURCE_thermal_externalheat_ID)) then
+     NofMyPhase=count(material_phase==phase)
+     instance = source_thermal_externalheat_instance(phase)
+     sourceOffset = source_thermal_externalheat_offset(phase)
+     source_thermal_externalheat_time(instance,1:source_thermal_externalheat_nIntervals(instance)+1_pInt) = &
+       temp_time(instance,1:source_thermal_externalheat_nIntervals(instance)+1_pInt)
+     source_thermal_externalheat_rate(instance,1:source_thermal_externalheat_nIntervals(instance)+1_pInt) = &
+       temp_rate(instance,1:source_thermal_externalheat_nIntervals(instance)+1_pInt)
+
+     sizeDotState              =   1_pInt
+     sizeDeltaState            =   0_pInt
+     sizeState                 =   1_pInt
+     sourceState(phase)%p(sourceOffset)%sizeState =       sizeState
+     sourceState(phase)%p(sourceOffset)%sizeDotState =    sizeDotState
+     sourceState(phase)%p(sourceOffset)%sizeDeltaState =  sizeDeltaState
+     sourceState(phase)%p(sourceOffset)%sizePostResults = source_thermal_externalheat_sizePostResults(instance)
+     allocate(sourceState(phase)%p(sourceOffset)%aTolState           (sizeState),                source=0.00001_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%state0              (sizeState,NofMyPhase),     source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%partionedState0     (sizeState,NofMyPhase),     source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%subState0           (sizeState,NofMyPhase),     source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%state               (sizeState,NofMyPhase),     source=0.0_pReal)
+
+     allocate(sourceState(phase)%p(sourceOffset)%dotState            (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%deltaState        (sizeDeltaState,NofMyPhase),  source=0.0_pReal)
+     if (.not. analyticJaco) then
+       allocate(sourceState(phase)%p(sourceOffset)%state_backup      (sizeState,NofMyPhase),     source=0.0_pReal)
+       allocate(sourceState(phase)%p(sourceOffset)%dotState_backup   (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     endif
+     if (any(numerics_integrator == 1_pInt)) then
+       allocate(sourceState(phase)%p(sourceOffset)%previousDotState  (sizeDotState,NofMyPhase),  source=0.0_pReal)
+       allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     endif
+     if (any(numerics_integrator == 4_pInt)) &
+       allocate(sourceState(phase)%p(sourceOffset)%RK4dotState       (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     if (any(numerics_integrator == 5_pInt)) &
+       allocate(sourceState(phase)%p(sourceOffset)%RKCK45dotState    (6,sizeDotState,NofMyPhase),source=0.0_pReal)      
+
+   endif
+ 
+ enddo initializeInstances
+end subroutine source_thermal_externalheat_init
+
+!--------------------------------------------------------------------------------------------------
+!> @brief calculates derived quantities from state
+!--------------------------------------------------------------------------------------------------
+subroutine source_thermal_externalheat_dotState(ipc, ip, el)
+ use material, only: &
+   phaseAt, phasememberAt, &
+   sourceState
+
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< component-ID of integration point
+   ip, &                                                                                            !< integration point
+   el                                                                                               !< element
+ integer(pInt) :: &
+   phase, &
+   constituent, &
+   sourceOffset
+
+ phase = phaseAt(ipc,ip,el)
+ constituent = phasememberAt(ipc,ip,el)
+ sourceOffset = source_thermal_externalheat_offset(phase)
+ 
+ sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = 1.0_pReal
+
+end subroutine source_thermal_externalheat_dotState
+
+!--------------------------------------------------------------------------------------------------
+!> @brief returns local vacancy generation rate 
+!--------------------------------------------------------------------------------------------------
+subroutine source_thermal_externalheat_getRateAndItsTangent(TDot, dTDot_dT, ipc, ip, el)
+ use material, only: &
+   phaseAt, phasememberAt, &
+   sourceState
+
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< grain number
+   ip, &                                                                                            !< integration point number
+   el                                                                                               !< element number
+ real(pReal),  intent(out) :: &
+   TDot, &
+   dTDot_dT
+ integer(pInt) :: &
+   instance, phase, constituent, sourceOffset, interval
+ real(pReal) :: &
+   norm_time   
+
+ phase = phaseAt(ipc,ip,el)
+ constituent = phasememberAt(ipc,ip,el)
+ instance = source_thermal_externalheat_instance(phase)
+ sourceOffset = source_thermal_externalheat_offset(phase)
+
+ do interval = 1, source_thermal_externalheat_nIntervals(instance)
+   norm_time = (sourceState(phase)%p(sourceOffset)%state(1,constituent) - &
+                source_thermal_externalheat_time(instance,interval)) / &
+               (source_thermal_externalheat_time(instance,interval+1) - &
+                source_thermal_externalheat_time(instance,interval))
+   if (norm_time >= 0.0_pReal .and. norm_time < 1.0_pReal) &
+     TDot = source_thermal_externalheat_rate(instance,interval  ) * (1.0_pReal - norm_time) + &
+            source_thermal_externalheat_rate(instance,interval+1) * norm_time
+ enddo
+ dTDot_dT = 0.0
+ 
+end subroutine source_thermal_externalheat_getRateAndItsTangent
+ 
+end module source_thermal_externalheat

code/source_vacancy_irradiation.f90

@@ -0,0 +1,253 @@
+!--------------------------------------------------------------------------------------------------
+! $Id$
+!--------------------------------------------------------------------------------------------------
+!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
+!> @brief material subroutine for vacancy generation due to irradiation
+!> @details to be done
+!--------------------------------------------------------------------------------------------------
+module source_vacancy_irradiation
+ use prec, only: &
+   pReal, &
+   pInt
+
+ implicit none
+ private
+ integer(pInt),                       dimension(:),           allocatable,         public, protected :: &
+   source_vacancy_irradiation_sizePostResults, &                                                        !< cumulative size of post results
+   source_vacancy_irradiation_offset, &                                                                 !< which source is my current damage mechanism?
+   source_vacancy_irradiation_instance                                                                  !< instance of damage source mechanism
+
+ integer(pInt),                       dimension(:,:),         allocatable, target, public :: &
+   source_vacancy_irradiation_sizePostResult                                                            !< size of each post result output
+
+ character(len=64),                   dimension(:,:),         allocatable, target, public :: &
+   source_vacancy_irradiation_output                                                                    !< name of each post result output
+   
+ integer(pInt),                       dimension(:),           allocatable, target, public :: &
+   source_vacancy_irradiation_Noutput                                                                   !< number of outputs per instance of this damage 
+
+ real(pReal),                         dimension(:),           allocatable,        private :: &
+   source_vacancy_irradiation_cascadeProb, &
+   source_vacancy_irradiation_cascadeVolume
+
+ public :: &
+   source_vacancy_irradiation_init, &
+   source_vacancy_irradiation_deltaState, &
+   source_vacancy_irradiation_getRateAndItsTangent
+
+contains
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief module initialization
+!> @details reads in material parameters, allocates arrays, and does sanity checks
+!--------------------------------------------------------------------------------------------------
+subroutine source_vacancy_irradiation_init(fileUnit)
+ use, intrinsic :: iso_fortran_env                                                                  ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
+ use debug, only: &
+   debug_level,&
+   debug_constitutive,&
+   debug_levelBasic
+ use IO, only: &
+   IO_read, &
+   IO_lc, &
+   IO_getTag, &
+   IO_isBlank, &
+   IO_stringPos, &
+   IO_stringValue, &
+   IO_floatValue, &
+   IO_intValue, &
+   IO_warning, &
+   IO_error, &
+   IO_timeStamp, &
+   IO_EOF
+ use material, only: &
+   phase_source, &
+   phase_Nsources, &
+   phase_Noutput, &
+   SOURCE_vacancy_irradiation_label, &
+   SOURCE_vacancy_irradiation_ID, &
+   material_Nphase, &
+   material_phase, &  
+   sourceState, &
+   MATERIAL_partPhase
+ use numerics,only: &
+   analyticJaco, &
+   worldrank, &
+   numerics_integrator
+
+ implicit none
+ integer(pInt), intent(in) :: fileUnit
+
+ integer(pInt), allocatable, dimension(:) :: chunkPos
+ integer(pInt) :: maxNinstance,phase,instance,source,sourceOffset
+ integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
+ integer(pInt) :: NofMyPhase   
+ character(len=65536) :: &
+   tag  = '', &
+   line = ''
+
+ mainProcess: if (worldrank == 0) then 
+   write(6,'(/,a)')   ' <<<+-  source_'//SOURCE_vacancy_irradiation_label//' init  -+>>>'
+   write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
+#include "compilation_info.f90"
+ endif mainProcess
+ 
+ maxNinstance = int(count(phase_source == SOURCE_vacancy_irradiation_ID),pInt)
+ if (maxNinstance == 0_pInt) return
+ if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
+   write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
+ 
+ allocate(source_vacancy_irradiation_offset(material_Nphase), source=0_pInt)
+ allocate(source_vacancy_irradiation_instance(material_Nphase), source=0_pInt)
+ do phase = 1, material_Nphase
+   source_vacancy_irradiation_instance(phase) = count(phase_source(:,1:phase) == source_vacancy_irradiation_ID)
+   do source = 1, phase_Nsources(phase)
+     if (phase_source(source,phase) == source_vacancy_irradiation_ID) &
+       source_vacancy_irradiation_offset(phase) = source
+   enddo    
+ enddo
+   
+ allocate(source_vacancy_irradiation_sizePostResults(maxNinstance),                     source=0_pInt)
+ allocate(source_vacancy_irradiation_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
+ allocate(source_vacancy_irradiation_output(maxval(phase_Noutput),maxNinstance))
+          source_vacancy_irradiation_output = ''
+ allocate(source_vacancy_irradiation_Noutput(maxNinstance),                             source=0_pInt) 
+ allocate(source_vacancy_irradiation_cascadeProb(maxNinstance),                         source=0.0_pReal) 
+ allocate(source_vacancy_irradiation_cascadeVolume(maxNinstance),                       source=0.0_pReal) 
+
+ rewind(fileUnit)
+ phase = 0_pInt
+ do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase)         ! wind forward to <phase>
+   line = IO_read(fileUnit)
+ enddo
+ 
+ parsingFile: do while (trim(line) /= IO_EOF)                                                       ! read through sections of phase part
+   line = IO_read(fileUnit)
+   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
+   if (IO_getTag(line,'<','>') /= '') then                                                          ! stop at next part
+     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
+     exit                                                                                           
+   endif   
+   if (IO_getTag(line,'[',']') /= '') then                                                          ! next phase section
+     phase = phase + 1_pInt                                                                         ! advance phase section counter
+     cycle                                                                                          ! skip to next line
+   endif
+
+   if (phase > 0_pInt ) then; if (any(phase_source(:,phase) == SOURCE_vacancy_irradiation_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
+
+     instance = source_vacancy_irradiation_instance(phase)                                          ! which instance of my vacancy is present phase
+     chunkPos = IO_stringPos(line)
+     tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt))                                             ! extract key
+     select case(tag)
+       case ('irradiation_cascadeprobability')
+         source_vacancy_irradiation_cascadeProb(instance) = IO_floatValue(line,chunkPos,2_pInt)
+
+       case ('irradiation_cascadevolume')
+         source_vacancy_irradiation_cascadeVolume(instance) = IO_floatValue(line,chunkPos,2_pInt)
+
+     end select
+   endif; endif
+ enddo parsingFile
+ 
+ initializeInstances: do phase = 1_pInt, material_Nphase
+   if (any(phase_source(:,phase) == SOURCE_vacancy_irradiation_ID)) then
+     NofMyPhase=count(material_phase==phase)
+     instance = source_vacancy_irradiation_instance(phase)
+     sourceOffset = source_vacancy_irradiation_offset(phase)
+
+     sizeDotState              =   2_pInt
+     sizeDeltaState            =   2_pInt
+     sizeState                 =   2_pInt
+     sourceState(phase)%p(sourceOffset)%sizeState =       sizeState
+     sourceState(phase)%p(sourceOffset)%sizeDotState =    sizeDotState
+     sourceState(phase)%p(sourceOffset)%sizeDeltaState =  sizeDeltaState
+     sourceState(phase)%p(sourceOffset)%sizePostResults = source_vacancy_irradiation_sizePostResults(instance)
+     allocate(sourceState(phase)%p(sourceOffset)%aTolState           (sizeState),                source=0.1_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%state0              (sizeState,NofMyPhase),     source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%partionedState0     (sizeState,NofMyPhase),     source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%subState0           (sizeState,NofMyPhase),     source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%state               (sizeState,NofMyPhase),     source=0.0_pReal)
+
+     allocate(sourceState(phase)%p(sourceOffset)%dotState            (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%deltaState        (sizeDeltaState,NofMyPhase),  source=0.0_pReal)
+     if (.not. analyticJaco) then
+       allocate(sourceState(phase)%p(sourceOffset)%state_backup      (sizeState,NofMyPhase),     source=0.0_pReal)
+       allocate(sourceState(phase)%p(sourceOffset)%dotState_backup   (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     endif
+     if (any(numerics_integrator == 1_pInt)) then
+       allocate(sourceState(phase)%p(sourceOffset)%previousDotState  (sizeDotState,NofMyPhase),  source=0.0_pReal)
+       allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     endif
+     if (any(numerics_integrator == 4_pInt)) &
+       allocate(sourceState(phase)%p(sourceOffset)%RK4dotState       (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     if (any(numerics_integrator == 5_pInt)) &
+       allocate(sourceState(phase)%p(sourceOffset)%RKCK45dotState    (6,sizeDotState,NofMyPhase),source=0.0_pReal)      
+
+   endif
+ 
+ enddo initializeInstances
+end subroutine source_vacancy_irradiation_init
+
+!--------------------------------------------------------------------------------------------------
+!> @brief calculates derived quantities from state
+!--------------------------------------------------------------------------------------------------
+subroutine source_vacancy_irradiation_deltaState(ipc, ip, el)
+ use material, only: &
+   phaseAt, phasememberAt, &
+   sourceState
+
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< component-ID of integration point
+   ip, &                                                                                            !< integration point
+   el                                                                                               !< element
+ integer(pInt) :: &
+   phase, constituent, sourceOffset
+ real(pReal) :: &
+   randNo   
+
+ phase       = phaseAt(ipc,ip,el)
+ constituent = phasememberAt(ipc,ip,el)
+ sourceOffset = source_vacancy_irradiation_offset(phase)
+
+ call random_number(randNo)
+ sourceState(phase)%p(sourceOffset)%deltaState(1,constituent) = &
+   randNo - sourceState(phase)%p(sourceOffset)%state(1,constituent)
+ call random_number(randNo)
+ sourceState(phase)%p(sourceOffset)%deltaState(2,constituent) = &
+   randNo - sourceState(phase)%p(sourceOffset)%state(2,constituent)
+
+end subroutine source_vacancy_irradiation_deltaState
+
+!--------------------------------------------------------------------------------------------------
+!> @brief returns local vacancy generation rate 
+!--------------------------------------------------------------------------------------------------
+subroutine source_vacancy_irradiation_getRateAndItsTangent(CvDot, dCvDot_dCv, ipc, ip, el)
+ use material, only: &
+   phaseAt, phasememberAt, &
+   sourceState
+
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< grain number
+   ip, &                                                                                            !< integration point number
+   el                                                                                               !< element number
+ real(pReal), intent(out) :: &
+   CvDot, dCvDot_dCv
+ integer(pInt) :: &
+   instance, phase, constituent, sourceOffset 
+
+ phase = phaseAt(ipc,ip,el)
+ constituent = phasememberAt(ipc,ip,el)
+ instance = source_vacancy_irradiation_instance(phase)
+ sourceOffset = source_vacancy_irradiation_offset(phase)
+ 
+ CvDot = 0.0_pReal
+ dCvDot_dCv = 0.0_pReal
+ if (sourceState(phase)%p(sourceOffset)%state0(1,constituent) < source_vacancy_irradiation_cascadeProb(instance)) &
+   CvDot = sourceState(phase)%p(sourceOffset)%state0(2,constituent)*source_vacancy_irradiation_cascadeVolume(instance)
+ 
+end subroutine source_vacancy_irradiation_getRateAndItsTangent
+ 
+end module source_vacancy_irradiation

code/source_vacancy_phenoplasticity.f90

@@ -0,0 +1,215 @@
+!--------------------------------------------------------------------------------------------------
+! $Id$
+!--------------------------------------------------------------------------------------------------
+!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
+!> @brief material subroutine for vacancy generation due to plasticity
+!> @details to be done
+!--------------------------------------------------------------------------------------------------
+module source_vacancy_phenoplasticity
+ use prec, only: &
+   pReal, &
+   pInt
+
+ implicit none
+ private
+ integer(pInt),                       dimension(:),           allocatable,         public, protected :: &
+   source_vacancy_phenoplasticity_sizePostResults, &                                                        !< cumulative size of post results
+   source_vacancy_phenoplasticity_offset, &                                                                 !< which source is my current damage mechanism?
+   source_vacancy_phenoplasticity_instance                                                                  !< instance of damage source mechanism
+
+ integer(pInt),                       dimension(:,:),         allocatable, target, public :: &
+   source_vacancy_phenoplasticity_sizePostResult                                                            !< size of each post result output
+
+ character(len=64),                   dimension(:,:),         allocatable, target, public :: &
+   source_vacancy_phenoplasticity_output                                                                    !< name of each post result output
+   
+ integer(pInt),                       dimension(:),           allocatable, target, public :: &
+   source_vacancy_phenoplasticity_Noutput                                                                   !< number of outputs per instance of this damage 
+
+ real(pReal),                         dimension(:),           allocatable,        private :: &
+   source_vacancy_phenoplasticity_rateCoeff
+
+ public :: &
+   source_vacancy_phenoplasticity_init, &
+   source_vacancy_phenoplasticity_getRateAndItsTangent
+
+contains
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief module initialization
+!> @details reads in material parameters, allocates arrays, and does sanity checks
+!--------------------------------------------------------------------------------------------------
+subroutine source_vacancy_phenoplasticity_init(fileUnit)
+ use, intrinsic :: iso_fortran_env                                                                  ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
+ use debug, only: &
+   debug_level,&
+   debug_constitutive,&
+   debug_levelBasic
+ use IO, only: &
+   IO_read, &
+   IO_lc, &
+   IO_getTag, &
+   IO_isBlank, &
+   IO_stringPos, &
+   IO_stringValue, &
+   IO_floatValue, &
+   IO_intValue, &
+   IO_warning, &
+   IO_error, &
+   IO_timeStamp, &
+   IO_EOF
+ use material, only: &
+   phase_source, &
+   phase_Nsources, &
+   phase_Noutput, &
+   SOURCE_vacancy_phenoplasticity_label, &
+   SOURCE_vacancy_phenoplasticity_ID, &
+   material_Nphase, &
+   material_phase, &  
+   sourceState, &
+   MATERIAL_partPhase
+ use numerics,only: &
+   analyticJaco, &
+   worldrank, &
+   numerics_integrator
+
+ implicit none
+ integer(pInt), intent(in) :: fileUnit
+
+ integer(pInt), allocatable, dimension(:) :: chunkPos
+ integer(pInt) :: maxNinstance,phase,instance,source,sourceOffset
+ integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
+ integer(pInt) :: NofMyPhase   
+ character(len=65536) :: &
+   tag  = '', &
+   line = ''
+
+ mainProcess: if (worldrank == 0) then 
+   write(6,'(/,a)')   ' <<<+-  source_'//SOURCE_vacancy_phenoplasticity_label//' init  -+>>>'
+   write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
+#include "compilation_info.f90"
+ endif mainProcess
+ 
+ maxNinstance = int(count(phase_source == SOURCE_vacancy_phenoplasticity_ID),pInt)
+ if (maxNinstance == 0_pInt) return
+ if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
+   write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
+ 
+ allocate(source_vacancy_phenoplasticity_offset(material_Nphase), source=0_pInt)
+ allocate(source_vacancy_phenoplasticity_instance(material_Nphase), source=0_pInt)
+ do phase = 1, material_Nphase
+   source_vacancy_phenoplasticity_instance(phase) = count(phase_source(:,1:phase) == source_vacancy_phenoplasticity_ID)
+   do source = 1, phase_Nsources(phase)
+     if (phase_source(source,phase) == source_vacancy_phenoplasticity_ID) &
+       source_vacancy_phenoplasticity_offset(phase) = source
+   enddo    
+ enddo
+   
+ allocate(source_vacancy_phenoplasticity_sizePostResults(maxNinstance),                     source=0_pInt)
+ allocate(source_vacancy_phenoplasticity_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
+ allocate(source_vacancy_phenoplasticity_output(maxval(phase_Noutput),maxNinstance))
+          source_vacancy_phenoplasticity_output = ''
+ allocate(source_vacancy_phenoplasticity_Noutput(maxNinstance),                             source=0_pInt) 
+ allocate(source_vacancy_phenoplasticity_rateCoeff(maxNinstance),                           source=0.0_pReal) 
+
+ rewind(fileUnit)
+ phase = 0_pInt
+ do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase)         ! wind forward to <phase>
+   line = IO_read(fileUnit)
+ enddo
+ 
+ parsingFile: do while (trim(line) /= IO_EOF)                                                       ! read through sections of phase part
+   line = IO_read(fileUnit)
+   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
+   if (IO_getTag(line,'<','>') /= '') then                                                          ! stop at next part
+     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
+     exit                                                                                           
+   endif   
+   if (IO_getTag(line,'[',']') /= '') then                                                          ! next phase section
+     phase = phase + 1_pInt                                                                         ! advance phase section counter
+     cycle                                                                                          ! skip to next line
+   endif
+
+   if (phase > 0_pInt ) then; if (any(phase_source(:,phase) == SOURCE_vacancy_phenoplasticity_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
+
+     instance = source_vacancy_phenoplasticity_instance(phase)                                          ! which instance of my vacancy is present phase
+     chunkPos = IO_stringPos(line)
+     tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt))                                             ! extract key
+     select case(tag)
+       case ('phenoplasticity_ratecoeff')
+         source_vacancy_phenoplasticity_rateCoeff(instance) = IO_floatValue(line,chunkPos,2_pInt)
+
+     end select
+   endif; endif
+ enddo parsingFile
+ 
+ initializeInstances: do phase = 1_pInt, material_Nphase
+   if (any(phase_source(:,phase) == SOURCE_vacancy_phenoplasticity_ID)) then
+     NofMyPhase=count(material_phase==phase)
+     instance = source_vacancy_phenoplasticity_instance(phase)
+     sourceOffset = source_vacancy_phenoplasticity_offset(phase)
+
+     sizeDotState              =   0_pInt
+     sizeDeltaState            =   0_pInt
+     sizeState                 =   0_pInt
+     sourceState(phase)%p(sourceOffset)%sizeState =       sizeState
+     sourceState(phase)%p(sourceOffset)%sizeDotState =    sizeDotState
+     sourceState(phase)%p(sourceOffset)%sizeDeltaState =  sizeDeltaState
+     sourceState(phase)%p(sourceOffset)%sizePostResults = source_vacancy_phenoplasticity_sizePostResults(instance)
+     allocate(sourceState(phase)%p(sourceOffset)%aTolState           (sizeState),                source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%state0              (sizeState,NofMyPhase),     source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%partionedState0     (sizeState,NofMyPhase),     source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%subState0           (sizeState,NofMyPhase),     source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%state               (sizeState,NofMyPhase),     source=0.0_pReal)
+
+     allocate(sourceState(phase)%p(sourceOffset)%dotState            (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%deltaState        (sizeDeltaState,NofMyPhase),  source=0.0_pReal)
+     if (.not. analyticJaco) then
+       allocate(sourceState(phase)%p(sourceOffset)%state_backup      (sizeState,NofMyPhase),     source=0.0_pReal)
+       allocate(sourceState(phase)%p(sourceOffset)%dotState_backup   (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     endif
+     if (any(numerics_integrator == 1_pInt)) then
+       allocate(sourceState(phase)%p(sourceOffset)%previousDotState  (sizeDotState,NofMyPhase),  source=0.0_pReal)
+       allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     endif
+     if (any(numerics_integrator == 4_pInt)) &
+       allocate(sourceState(phase)%p(sourceOffset)%RK4dotState       (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     if (any(numerics_integrator == 5_pInt)) &
+       allocate(sourceState(phase)%p(sourceOffset)%RKCK45dotState    (6,sizeDotState,NofMyPhase),source=0.0_pReal)      
+
+   endif
+ 
+ enddo initializeInstances
+end subroutine source_vacancy_phenoplasticity_init
+
+!--------------------------------------------------------------------------------------------------
+!> @brief returns local vacancy generation rate 
+!--------------------------------------------------------------------------------------------------
+subroutine source_vacancy_phenoplasticity_getRateAndItsTangent(CvDot, dCvDot_dCv, ipc, ip, el)
+ use material, only: &
+   phaseAt, phasememberAt, &
+   plasticState
+
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< grain number
+   ip, &                                                                                            !< integration point number
+   el                                                                                               !< element number
+ real(pReal), intent(out) :: &
+   CvDot, dCvDot_dCv
+ integer(pInt) :: &
+   instance, phase, constituent 
+
+ phase = phaseAt(ipc,ip,el)
+ constituent = phasememberAt(ipc,ip,el)
+ instance = source_vacancy_phenoplasticity_instance(phase)
+ 
+ CvDot = &
+   source_vacancy_phenoplasticity_rateCoeff(instance)* &
+   sum(plasticState(phase)%slipRate(:,constituent))
+ dCvDot_dCv = 0.0_pReal  
+ 
+end subroutine source_vacancy_phenoplasticity_getRateAndItsTangent
+ 
+end module source_vacancy_phenoplasticity

code/source_vacancy_thermalfluc.f90

@@ -0,0 +1,255 @@
+!--------------------------------------------------------------------------------------------------
+! $Id$
+!--------------------------------------------------------------------------------------------------
+!> @author Pratheek Shanthraj, Max-Planck-Institut für Eisenforschung GmbH
+!> @brief material subroutine for vacancy generation due to thermal fluctuations
+!> @details to be done
+!--------------------------------------------------------------------------------------------------
+module source_vacancy_thermalfluc
+ use prec, only: &
+   pReal, &
+   pInt
+
+ implicit none
+ private
+ integer(pInt),                       dimension(:),           allocatable,         public, protected :: &
+   source_vacancy_thermalfluc_sizePostResults, &                                                        !< cumulative size of post results
+   source_vacancy_thermalfluc_offset, &                                                                 !< which source is my current damage mechanism?
+   source_vacancy_thermalfluc_instance                                                                  !< instance of damage source mechanism
+
+ integer(pInt),                       dimension(:,:),         allocatable, target, public :: &
+   source_vacancy_thermalfluc_sizePostResult                                                            !< size of each post result output
+
+ character(len=64),                   dimension(:,:),         allocatable, target, public :: &
+   source_vacancy_thermalfluc_output                                                                    !< name of each post result output
+   
+ integer(pInt),                       dimension(:),           allocatable, target, public :: &
+   source_vacancy_thermalfluc_Noutput                                                                   !< number of outputs per instance of this damage 
+
+ real(pReal),                         dimension(:),           allocatable,        private :: &
+   source_vacancy_thermalfluc_amplitude, &
+   source_vacancy_thermalfluc_normVacancyEnergy
+
+ public :: &
+   source_vacancy_thermalfluc_init, &
+   source_vacancy_thermalfluc_deltaState, &
+   source_vacancy_thermalfluc_getRateAndItsTangent
+
+contains
+
+
+!--------------------------------------------------------------------------------------------------
+!> @brief module initialization
+!> @details reads in material parameters, allocates arrays, and does sanity checks
+!--------------------------------------------------------------------------------------------------
+subroutine source_vacancy_thermalfluc_init(fileUnit)
+ use, intrinsic :: iso_fortran_env                                                                  ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
+ use debug, only: &
+   debug_level,&
+   debug_constitutive,&
+   debug_levelBasic
+ use IO, only: &
+   IO_read, &
+   IO_lc, &
+   IO_getTag, &
+   IO_isBlank, &
+   IO_stringPos, &
+   IO_stringValue, &
+   IO_floatValue, &
+   IO_intValue, &
+   IO_warning, &
+   IO_error, &
+   IO_timeStamp, &
+   IO_EOF
+ use lattice, only: &
+   lattice_vacancyFormationEnergy
+ use material, only: &
+   phase_source, &
+   phase_Nsources, &
+   phase_Noutput, &
+   SOURCE_vacancy_thermalfluc_label, &
+   SOURCE_vacancy_thermalfluc_ID, &
+   material_Nphase, &
+   material_phase, &  
+   sourceState, &
+   MATERIAL_partPhase
+ use numerics,only: &
+   analyticJaco, &
+   worldrank, &
+   numerics_integrator
+
+ implicit none
+ integer(pInt), intent(in) :: fileUnit
+
+ integer(pInt), allocatable, dimension(:) :: chunkPos
+ integer(pInt) :: maxNinstance,phase,instance,source,sourceOffset
+ integer(pInt) :: sizeState, sizeDotState, sizeDeltaState
+ integer(pInt) :: NofMyPhase   
+ character(len=65536) :: &
+   tag  = '', &
+   line = ''
+
+ mainProcess: if (worldrank == 0) then 
+   write(6,'(/,a)')   ' <<<+-  source_'//SOURCE_vacancy_thermalfluc_label//' init  -+>>>'
+   write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
+#include "compilation_info.f90"
+ endif mainProcess
+ 
+ maxNinstance = int(count(phase_source == SOURCE_vacancy_thermalfluc_ID),pInt)
+ if (maxNinstance == 0_pInt) return
+ if (iand(debug_level(debug_constitutive),debug_levelBasic) /= 0_pInt) &
+   write(6,'(a16,1x,i5,/)') '# instances:',maxNinstance
+ 
+ allocate(source_vacancy_thermalfluc_offset(material_Nphase), source=0_pInt)
+ allocate(source_vacancy_thermalfluc_instance(material_Nphase), source=0_pInt)
+ do phase = 1, material_Nphase
+   source_vacancy_thermalfluc_instance(phase) = count(phase_source(:,1:phase) == source_vacancy_thermalfluc_ID)
+   do source = 1, phase_Nsources(phase)
+     if (phase_source(source,phase) == source_vacancy_thermalfluc_ID) &
+       source_vacancy_thermalfluc_offset(phase) = source
+   enddo    
+ enddo
+   
+ allocate(source_vacancy_thermalfluc_sizePostResults(maxNinstance),                     source=0_pInt)
+ allocate(source_vacancy_thermalfluc_sizePostResult(maxval(phase_Noutput),maxNinstance),source=0_pInt)
+ allocate(source_vacancy_thermalfluc_output(maxval(phase_Noutput),maxNinstance))
+          source_vacancy_thermalfluc_output = ''
+ allocate(source_vacancy_thermalfluc_Noutput(maxNinstance),                             source=0_pInt) 
+ allocate(source_vacancy_thermalfluc_amplitude(maxNinstance),                           source=0.0_pReal) 
+ allocate(source_vacancy_thermalfluc_normVacancyEnergy(maxNinstance),                   source=0.0_pReal) 
+
+ rewind(fileUnit)
+ phase = 0_pInt
+ do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>')) /= MATERIAL_partPhase)         ! wind forward to <phase>
+   line = IO_read(fileUnit)
+ enddo
+ 
+ parsingFile: do while (trim(line) /= IO_EOF)                                                       ! read through sections of phase part
+   line = IO_read(fileUnit)
+   if (IO_isBlank(line)) cycle                                                                      ! skip empty lines
+   if (IO_getTag(line,'<','>') /= '') then                                                          ! stop at next part
+     line = IO_read(fileUnit, .true.)                                                               ! reset IO_read
+     exit                                                                                           
+   endif   
+   if (IO_getTag(line,'[',']') /= '') then                                                          ! next phase section
+     phase = phase + 1_pInt                                                                         ! advance phase section counter
+     cycle                                                                                          ! skip to next line
+   endif
+
+   if (phase > 0_pInt ) then; if (any(phase_source(:,phase) == SOURCE_vacancy_thermalfluc_ID)) then ! do not short-circuit here (.and. with next if statemen). It's not safe in Fortran
+
+     instance = source_vacancy_thermalfluc_instance(phase)                                          ! which instance of my vacancy is present phase
+     chunkPos = IO_stringPos(line)
+     tag = IO_lc(IO_stringValue(line,chunkPos,1_pInt))                                             ! extract key
+     select case(tag)
+       case ('thermalfluctuation_amplitude')
+         source_vacancy_thermalfluc_amplitude(instance) = IO_floatValue(line,chunkPos,2_pInt)
+
+     end select
+   endif; endif
+ enddo parsingFile
+ 
+ initializeInstances: do phase = 1_pInt, material_Nphase
+   if (any(phase_source(:,phase) == SOURCE_vacancy_thermalfluc_ID)) then
+     NofMyPhase=count(material_phase==phase)
+     instance = source_vacancy_thermalfluc_instance(phase)
+     source_vacancy_thermalfluc_normVacancyEnergy(instance) = &
+       lattice_vacancyFormationEnergy(phase)/1.3806488e-23_pReal
+     sourceOffset = source_vacancy_thermalfluc_offset(phase)
+
+     sizeDotState              =   1_pInt
+     sizeDeltaState            =   1_pInt
+     sizeState                 =   1_pInt
+     sourceState(phase)%p(sourceOffset)%sizeState =       sizeState
+     sourceState(phase)%p(sourceOffset)%sizeDotState =    sizeDotState
+     sourceState(phase)%p(sourceOffset)%sizeDeltaState =  sizeDeltaState
+     sourceState(phase)%p(sourceOffset)%sizePostResults = source_vacancy_thermalfluc_sizePostResults(instance)
+     allocate(sourceState(phase)%p(sourceOffset)%aTolState           (sizeState),                source=0.1_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%state0              (sizeState,NofMyPhase),     source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%partionedState0     (sizeState,NofMyPhase),     source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%subState0           (sizeState,NofMyPhase),     source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%state               (sizeState,NofMyPhase),     source=0.0_pReal)
+
+     allocate(sourceState(phase)%p(sourceOffset)%dotState            (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     allocate(sourceState(phase)%p(sourceOffset)%deltaState        (sizeDeltaState,NofMyPhase),  source=0.0_pReal)
+     if (.not. analyticJaco) then
+       allocate(sourceState(phase)%p(sourceOffset)%state_backup      (sizeState,NofMyPhase),     source=0.0_pReal)
+       allocate(sourceState(phase)%p(sourceOffset)%dotState_backup   (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     endif
+     if (any(numerics_integrator == 1_pInt)) then
+       allocate(sourceState(phase)%p(sourceOffset)%previousDotState  (sizeDotState,NofMyPhase),  source=0.0_pReal)
+       allocate(sourceState(phase)%p(sourceOffset)%previousDotState2 (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     endif
+     if (any(numerics_integrator == 4_pInt)) &
+       allocate(sourceState(phase)%p(sourceOffset)%RK4dotState       (sizeDotState,NofMyPhase),  source=0.0_pReal)
+     if (any(numerics_integrator == 5_pInt)) &
+       allocate(sourceState(phase)%p(sourceOffset)%RKCK45dotState    (6,sizeDotState,NofMyPhase),source=0.0_pReal)      
+
+   endif
+ 
+ enddo initializeInstances
+end subroutine source_vacancy_thermalfluc_init
+
+!--------------------------------------------------------------------------------------------------
+!> @brief calculates derived quantities from state
+!--------------------------------------------------------------------------------------------------
+subroutine source_vacancy_thermalfluc_deltaState(ipc, ip, el)
+ use material, only: &
+   phaseAt, phasememberAt, &
+   sourceState
+
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< component-ID of integration point
+   ip, &                                                                                            !< integration point
+   el                                                                                               !< element
+ integer(pInt) :: &
+   phase, constituent, sourceOffset 
+ real(pReal) :: &
+   randNo   
+
+ phase       = phaseAt(ipc,ip,el)
+ constituent = phasememberAt(ipc,ip,el)
+ sourceOffset = source_vacancy_thermalfluc_offset(phase)
+
+ call random_number(randNo)
+ sourceState(phase)%p(sourceOffset)%deltaState(1,constituent) = &
+   randNo - 0.5_pReal - sourceState(phase)%p(sourceOffset)%state(1,constituent)
+
+end subroutine source_vacancy_thermalfluc_deltaState
+
+!--------------------------------------------------------------------------------------------------
+!> @brief returns local vacancy generation rate 
+!--------------------------------------------------------------------------------------------------
+subroutine source_vacancy_thermalfluc_getRateAndItsTangent(CvDot, dCvDot_dCv, ipc, ip, el)
+ use material, only: &
+   phaseAt, phasememberAt, &
+   material_homog, &
+   temperature, &
+   thermalMapping, &
+   sourceState
+
+ implicit none
+ integer(pInt), intent(in) :: &
+   ipc, &                                                                                           !< grain number
+   ip, &                                                                                            !< integration point number
+   el                                                                                               !< element number
+ real(pReal), intent(out) :: &
+   CvDot, dCvDot_dCv
+ integer(pInt) :: &
+   instance, phase, constituent, sourceOffset 
+
+ phase = phaseAt(ipc,ip,el)
+ constituent = phasememberAt(ipc,ip,el)
+ instance = source_vacancy_thermalfluc_instance(phase)
+ sourceOffset = source_vacancy_thermalfluc_offset(phase)
+ 
+ CvDot = source_vacancy_thermalfluc_amplitude(instance)* &
+         sourceState(phase)%p(sourceOffset)%state0(2,constituent)* &
+         exp(-source_vacancy_thermalfluc_normVacancyEnergy(instance)/ &
+              temperature(material_homog(ip,el))%p(thermalMapping(material_homog(ip,el))%p(ip,el)))
+ dCvDot_dCv = 0.0_pReal
+ 
+end subroutine source_vacancy_thermalfluc_getRateAndItsTangent
+ 
+end module source_vacancy_thermalfluc

compile.log

@@ -0,0 +1,3 @@
+make DAMASK_spectral.exe -C code
+make[1]: Entering directory `/mnt/research/CMM/DAMASK/zhangc43_git/code'
+make[1]: Leaving directory `/mnt/research/CMM/DAMASK/zhangc43_git/code'