constitutive.f90

Go to the documentation of this file.

# 1 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/constitutive.f90"
# 1 "<built-in>"
# 1 "<command-line>"
# 1 "/home/damask_user/GitLabCI_Pipeline_4301/DAMASK/src/constitutive.f90"
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
module constitutive
  use prec
  use math
  use rotations
  use debug
  use numerics
  use io
  use config
  use material
  use results
  use lattice
  use discretization
  use geometry_plastic_nonlocal
  use source_thermal_dissipation
  use source_thermal_externalheat
  use source_damage_isobrittle
  use source_damage_isoductile
  use source_damage_anisobrittle
  use source_damage_anisoductile
  use kinematics_cleavage_opening
  use kinematics_slipplane_opening
  use kinematics_thermal_expansion
 
  implicit none
  private
 
  integer, public, protected :: &
    constitutive_plasticity_maxsizedotstate, &
    constitutive_source_maxsizedotstate
 
  interface
 
    module subroutine plastic_none_init
    end subroutine plastic_none_init
 
    module subroutine plastic_isotropic_init
    end subroutine plastic_isotropic_init
 
    module subroutine plastic_phenopowerlaw_init
    end subroutine plastic_phenopowerlaw_init
 
    module subroutine plastic_kinehardening_init
    end subroutine plastic_kinehardening_init
 
    module subroutine plastic_dislotwin_init
    end subroutine plastic_dislotwin_init
 
    module subroutine plastic_disloucla_init
    end subroutine plastic_disloucla_init
 
    module subroutine plastic_nonlocal_init
    end subroutine plastic_nonlocal_init
 
 
    module subroutine plastic_isotropic_lpanditstangent(lp,dlp_dmp,mp,instance,of)
      real(preal), dimension(3,3),     intent(out) :: &
        lp
      real(preal), dimension(3,3,3,3), intent(out) :: &
        dlp_dmp
 
      real(preal), dimension(3,3),     intent(in) :: &
        mp
      integer,                         intent(in) :: &
        instance, &
        of
    end subroutine plastic_isotropic_lpanditstangent
 
    pure module subroutine plastic_phenopowerlaw_lpanditstangent(lp,dlp_dmp,mp,instance,of)
      real(preal), dimension(3,3),     intent(out) :: &
        lp
      real(preal), dimension(3,3,3,3), intent(out) :: &
        dlp_dmp
 
      real(preal), dimension(3,3),     intent(in) :: &
        mp
      integer,                         intent(in) :: &
        instance, &
        of
    end subroutine plastic_phenopowerlaw_lpanditstangent
 
    pure module subroutine plastic_kinehardening_lpanditstangent(lp,dlp_dmp,mp,instance,of)
      real(preal), dimension(3,3),     intent(out) :: &
        lp
      real(preal), dimension(3,3,3,3), intent(out) :: &
        dlp_dmp
 
      real(preal), dimension(3,3),     intent(in) :: &
        mp
      integer,                         intent(in) :: &
        instance, &
        of
    end subroutine plastic_kinehardening_lpanditstangent
 
    module subroutine plastic_dislotwin_lpanditstangent(lp,dlp_dmp,mp,t,instance,of)
      real(preal), dimension(3,3),     intent(out) :: &
        lp
      real(preal), dimension(3,3,3,3), intent(out) :: &
        dlp_dmp
 
      real(preal), dimension(3,3),     intent(in) :: &
        mp
      real(preal),                     intent(in) :: &
        t
      integer,                         intent(in) :: &
        instance, &
        of
    end subroutine plastic_dislotwin_lpanditstangent
 
    pure module subroutine plastic_disloucla_lpanditstangent(lp,dlp_dmp,mp,t,instance,of)
      real(preal), dimension(3,3),     intent(out) :: &
        lp
      real(preal), dimension(3,3,3,3), intent(out) :: &
        dlp_dmp
 
      real(preal), dimension(3,3),     intent(in) :: &
        mp
      real(preal),                     intent(in) :: &
        t
      integer,                         intent(in) :: &
        instance, &
        of
    end subroutine plastic_disloucla_lpanditstangent
 
    module subroutine plastic_nonlocal_lpanditstangent(lp,dlp_dmp, &
                                                       mp,temperature,instance,of,ip,el)
      real(preal), dimension(3,3),     intent(out) :: &
        lp
      real(preal), dimension(3,3,3,3), intent(out) :: &
        dlp_dmp
 
      real(preal), dimension(3,3),     intent(in) :: &
        mp
      real(preal),                     intent(in) :: &
        temperature
      integer,                         intent(in) :: &
        instance, &
        of, &
        ip, &                                                                                       !< current integration point
        el
    end subroutine plastic_nonlocal_lpanditstangent
 
 
    module subroutine plastic_isotropic_lianditstangent(li,dli_dmi,mi,instance,of)
      real(preal), dimension(3,3),     intent(out) :: &
        li
      real(preal), dimension(3,3,3,3), intent(out)  :: &
        dli_dmi
 
      real(preal), dimension(3,3),     intent(in) :: &
        mi
      integer,                         intent(in) :: &
        instance, &
        of
    end subroutine plastic_isotropic_lianditstangent
 
 
    module subroutine plastic_isotropic_dotstate(mp,instance,of)
      real(preal), dimension(3,3),  intent(in) :: &
        mp
      integer,                      intent(in) :: &
        instance, &
        of
    end subroutine plastic_isotropic_dotstate
 
    module subroutine plastic_phenopowerlaw_dotstate(mp,instance,of)
      real(preal), dimension(3,3),  intent(in) :: &
        mp
      integer,                      intent(in) :: &
        instance, &
        of
    end subroutine plastic_phenopowerlaw_dotstate
 
    module subroutine plastic_kinehardening_dotstate(mp,instance,of)
      real(preal), dimension(3,3),  intent(in) :: &
        mp
      integer,                      intent(in) :: &
        instance, &
        of
    end subroutine plastic_kinehardening_dotstate
 
    module subroutine plastic_dislotwin_dotstate(mp,t,instance,of)
      real(preal), dimension(3,3),  intent(in) :: &
        mp
      real(preal),                  intent(in) :: &
        t
      integer,                      intent(in) :: &
        instance, &
        of
    end subroutine plastic_dislotwin_dotstate
 
    module subroutine plastic_disloucla_dotstate(mp,t,instance,of)
      real(preal), dimension(3,3),  intent(in) :: &
        mp
      real(preal),                  intent(in) :: &
        t
      integer,                      intent(in) :: &
        instance, &
        of
    end subroutine plastic_disloucla_dotstate
 
    module subroutine plastic_nonlocal_dotstate(mp, f, fp, temperature,timestep, &
                                                instance,of,ip,el)
      real(preal), dimension(3,3), intent(in) ::&
        mp
      real(preal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: &
        f, &                                                                                        !< deformation gradient
        fp
      real(preal), intent(in) :: &
        temperature, &                                                                              !< temperature
        timestep
      integer, intent(in) :: &
        instance, &
        of, &
        ip, &                                                                                       !< current integration point
        el
    end subroutine plastic_nonlocal_dotstate
 
 
    module subroutine plastic_dislotwin_dependentstate(t,instance,of)
      integer,       intent(in) :: &
        instance, &
        of
      real(preal),   intent(in) :: &
        t
    end subroutine plastic_dislotwin_dependentstate
 
    module subroutine plastic_disloucla_dependentstate(instance,of)
      integer,       intent(in) :: &
        instance, &
        of
    end subroutine plastic_disloucla_dependentstate
 
    module subroutine plastic_nonlocal_dependentstate(f, fp, instance, of, ip, el)
      real(preal), dimension(3,3), intent(in) :: &
        f, &
        fp
      integer, intent(in) :: &
        instance, &
        of, &
        ip, &
        el
    end subroutine plastic_nonlocal_dependentstate
 
 
    module subroutine plastic_kinehardening_deltastate(mp,instance,of)
      real(preal), dimension(3,3),  intent(in) :: &
        mp
      integer,                      intent(in) :: &
        instance, &
        of
    end subroutine plastic_kinehardening_deltastate
 
    module subroutine plastic_nonlocal_deltastate(mp,instance,of,ip,el)
      real(preal), dimension(3,3), intent(in) :: &
        mp
      integer, intent(in) :: &
        instance, &
        of, &
        ip, &
        el
    end subroutine plastic_nonlocal_deltastate
 
 
    module function plastic_dislotwin_homogenizedc(ipc,ip,el) result(homogenizedc)
      real(preal), dimension(6,6) :: &
        homogenizedc
      integer,     intent(in) :: &
        ipc, &                                                                                      !< component-ID of integration point
        ip, &                                                                                       !< integration point
        el
    end function plastic_dislotwin_homogenizedc
 
    module subroutine plastic_nonlocal_updatecompatibility(orientation,instance,i,e)
      integer, intent(in) :: &
        instance, &
        i, &
        e
      type(rotation), dimension(1,discretization_nIP,discretization_nElem), intent(in) :: &
        orientation
    end subroutine plastic_nonlocal_updatecompatibility
 
 
    module subroutine plastic_isotropic_results(instance,group)
      integer,          intent(in) :: instance
      character(len=*), intent(in) :: group
    end subroutine plastic_isotropic_results
 
    module subroutine plastic_phenopowerlaw_results(instance,group)
      integer,          intent(in) :: instance
      character(len=*), intent(in) :: group
    end subroutine plastic_phenopowerlaw_results
 
    module subroutine plastic_kinehardening_results(instance,group)
      integer,          intent(in) :: instance
      character(len=*), intent(in) :: group
    end subroutine plastic_kinehardening_results
 
    module subroutine plastic_dislotwin_results(instance,group)
      integer,          intent(in) :: instance
      character(len=*), intent(in) :: group
    end subroutine plastic_dislotwin_results
 
    module subroutine plastic_disloucla_results(instance,group)
      integer,          intent(in) :: instance
      character(len=*), intent(in) :: group
    end subroutine plastic_disloucla_results
 
    module subroutine plastic_nonlocal_results(instance,group)
      integer,          intent(in) :: instance
      character(len=*), intent(in) :: group
    end subroutine plastic_nonlocal_results
 
  end interface
 
  public :: &
    plastic_nonlocal_updatecompatibility, &
    constitutive_init, &
    constitutive_homogenizedc, &
    constitutive_dependentstate, &
    constitutive_lpanditstangents, &
    constitutive_lianditstangents, &
    constitutive_initialfi, &
    constitutive_sanditstangents, &
    constitutive_collectdotstate, &
    constitutive_collectdeltastate, &
    constitutive_results
 
contains
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
subroutine constitutive_init
 
  integer :: &
    ph, &                                                                                           !< counter in phase loop
    s
 
!--------------------------------------------------------------------------------------------------
! initialized plasticity
  if (any(phase_plasticity == plasticity_none_id))          call plastic_none_init
  if (any(phase_plasticity == plasticity_isotropic_id))     call plastic_isotropic_init
  if (any(phase_plasticity == plasticity_phenopowerlaw_id)) call plastic_phenopowerlaw_init
  if (any(phase_plasticity == plasticity_kinehardening_id)) call plastic_kinehardening_init
  if (any(phase_plasticity == plasticity_dislotwin_id))     call plastic_dislotwin_init
  if (any(phase_plasticity == plasticity_disloucla_id))     call plastic_disloucla_init
  if (any(phase_plasticity == plasticity_nonlocal_id)) then
    call plastic_nonlocal_init
  else
    call geometry_plastic_nonlocal_disable
  endif
!--------------------------------------------------------------------------------------------------
! initialize source mechanisms
  if (any(phase_source == source_thermal_dissipation_id))     call source_thermal_dissipation_init
  if (any(phase_source == source_thermal_externalheat_id))    call source_thermal_externalheat_init
  if (any(phase_source == source_damage_isobrittle_id))       call source_damage_isobrittle_init
  if (any(phase_source == source_damage_isoductile_id))       call source_damage_isoductile_init
  if (any(phase_source == source_damage_anisobrittle_id))     call source_damage_anisobrittle_init
  if (any(phase_source == source_damage_anisoductile_id))     call source_damage_anisoductile_init
 
!--------------------------------------------------------------------------------------------------
! initialize kinematic mechanisms
  if (any(phase_kinematics == kinematics_cleavage_opening_id))  call kinematics_cleavage_opening_init
  if (any(phase_kinematics == kinematics_slipplane_opening_id)) call kinematics_slipplane_opening_init
  if (any(phase_kinematics == kinematics_thermal_expansion_id)) call kinematics_thermal_expansion_init
 
  write(6,'(/,a)')   ' <<<+-  constitutive init  -+>>>'; flush(6)
 
  constitutive_source_maxsizedotstate = 0
  phaseloop2:do ph = 1,material_nphase
!--------------------------------------------------------------------------------------------------
! partition and inititalize state
    plasticstate(ph)%partionedState0 = plasticstate(ph)%state0
    plasticstate(ph)%state           = plasticstate(ph)%partionedState0
    forall(s = 1:phase_nsources(ph))
      sourcestate(ph)%p(s)%partionedState0 = sourcestate(ph)%p(s)%state0
      sourcestate(ph)%p(s)%state           = sourcestate(ph)%p(s)%partionedState0
    end forall
!--------------------------------------------------------------------------------------------------
! determine max size of source state
    constitutive_source_maxsizedotstate   = max(constitutive_source_maxsizedotstate, &
                                                maxval(sourcestate(ph)%p%sizeDotState))
  enddo phaseloop2
  constitutive_plasticity_maxsizedotstate = maxval(plasticstate%sizeDotState)
 
end subroutine constitutive_init
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
function constitutive_homogenizedc(ipc,ip,el)
 
  real(preal), dimension(6,6) :: constitutive_homogenizedc
  integer, intent(in) :: &
    ipc, &                                                                                          !< component-ID of integration point
    ip, &                                                                                           !< integration point
    el
 
  plasticitytype: select case (phase_plasticity(material_phaseat(ipc,el)))
    case (plasticity_dislotwin_id) plasticitytype
      constitutive_homogenizedc = plastic_dislotwin_homogenizedc(ipc,ip,el)
    case default plasticitytype
      constitutive_homogenizedc = lattice_c66(1:6,1:6,material_phaseat(ipc,el))
  end select plasticitytype
 
end function constitutive_homogenizedc
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
subroutine constitutive_dependentstate(F, Fp, ipc, ip, el)
 
  integer, intent(in) :: &
    ipc, &                                                                                          !< component-ID of integration point
    ip, &                                                                                           !< integration point
    el
  real(preal),   intent(in), dimension(3,3) :: &
    f, &                                                                                           !< elastic deformation gradient
    fp
  integer :: &
    ho, &                                                                                           !< homogenization
    tme, &                                                                                          !< thermal member position
    instance, of
 
  ho  = material_homogenizationat(el)
  tme = thermalmapping(ho)%p(ip,el)
  of  = material_phasememberat(ipc,ip,el)
  instance = phase_plasticityinstance(material_phaseat(ipc,el))
 
  plasticitytype: select case (phase_plasticity(material_phaseat(ipc,el)))
    case (plasticity_dislotwin_id) plasticitytype
      call plastic_dislotwin_dependentstate(temperature(ho)%p(tme),instance,of)
    case (plasticity_disloucla_id) plasticitytype
      call plastic_disloucla_dependentstate(instance,of)
    case (plasticity_nonlocal_id) plasticitytype
      call plastic_nonlocal_dependentstate (f,fp,instance,of,ip,el)
  end select plasticitytype
 
end subroutine constitutive_dependentstate
 
 
!--------------------------------------------------------------------------------------------------
! ToDo: Discuss wheter it makes sense if crystallite handles the configuration conversion, i.e.
! Mp in, dLp_dMp out
!--------------------------------------------------------------------------------------------------
subroutine constitutive_lpanditstangents(Lp, dLp_dS, dLp_dFi, &
                                         S, Fi, ipc, ip, el)
  integer, intent(in) :: &
    ipc, &                                                                                          !< component-ID of integration point
    ip, &                                                                                           !< integration point
    el
  real(preal),   intent(in),  dimension(3,3) :: &
    s, &                                                                                            !< 2nd Piola-Kirchhoff stress
    fi
  real(preal),   intent(out), dimension(3,3) :: &
    lp
  real(preal),   intent(out), dimension(3,3,3,3) :: &
    dlp_ds, &
    dlp_dfi
  real(preal), dimension(3,3,3,3) :: &
    dlp_dmp
  real(preal), dimension(3,3) :: &
    mp
  integer :: &
    ho, &                                                                                           !< homogenization
    tme
  integer :: &
    i, j, instance, of
 
  ho = material_homogenizationat(el)
  tme = thermalmapping(ho)%p(ip,el)
 
  mp = matmul(matmul(transpose(fi),fi),s)
  of = material_phasememberat(ipc,ip,el)
  instance = phase_plasticityinstance(material_phaseat(ipc,el))
 
  plasticitytype: select case (phase_plasticity(material_phaseat(ipc,el)))
 
    case (plasticity_none_id) plasticitytype
      lp = 0.0_preal
      dlp_dmp = 0.0_preal
 
    case (plasticity_isotropic_id) plasticitytype
      call plastic_isotropic_lpanditstangent   (lp,dlp_dmp,mp,instance,of)
 
    case (plasticity_phenopowerlaw_id) plasticitytype
      call plastic_phenopowerlaw_lpanditstangent(lp,dlp_dmp,mp,instance,of)
 
    case (plasticity_kinehardening_id) plasticitytype
      call plastic_kinehardening_lpanditstangent(lp,dlp_dmp,mp,instance,of)
 
    case (plasticity_nonlocal_id) plasticitytype
      call plastic_nonlocal_lpanditstangent     (lp,dlp_dmp,mp, temperature(ho)%p(tme),instance,of,ip,el)
 
    case (plasticity_dislotwin_id) plasticitytype
      call plastic_dislotwin_lpanditstangent    (lp,dlp_dmp,mp,temperature(ho)%p(tme),instance,of)
 
    case (plasticity_disloucla_id) plasticitytype
      call plastic_disloucla_lpanditstangent    (lp,dlp_dmp,mp,temperature(ho)%p(tme),instance,of)
 
  end select plasticitytype
 
  do i=1,3; do j=1,3
    dlp_dfi(i,j,1:3,1:3) = matmul(matmul(fi,s),transpose(dlp_dmp(i,j,1:3,1:3))) + &
                           matmul(matmul(fi,dlp_dmp(i,j,1:3,1:3)),s)
    dlp_ds(i,j,1:3,1:3)  = matmul(matmul(transpose(fi),fi),dlp_dmp(i,j,1:3,1:3))                     ! ToDo: @PS: why not:   dLp_dMp:(FiT Fi)
  enddo; enddo
 
end subroutine constitutive_lpanditstangents
 
 
!--------------------------------------------------------------------------------------------------
! ToDo: MD: S is Mi?
!--------------------------------------------------------------------------------------------------
subroutine constitutive_lianditstangents(Li, dLi_dS, dLi_dFi, &
                                         S, Fi, ipc, ip, el)
 
  integer, intent(in) :: &
    ipc, &                                                                                          !< component-ID of integration point
    ip, &                                                                                           !< integration point
    el
  real(preal),   intent(in),  dimension(3,3) :: &
    s
  real(preal),   intent(in),  dimension(3,3) :: &
    fi
  real(preal),   intent(out), dimension(3,3) :: &
    li
  real(preal),   intent(out), dimension(3,3,3,3) :: &
    dli_ds, &                                                                                       !< derivative of Li with respect to S
    dli_dfi
 
  real(preal), dimension(3,3) :: &
    my_li, &                                                                                        !< intermediate velocity gradient
    fiinv, &
    temp_33
  real(preal), dimension(3,3,3,3) :: &
    my_dli_ds
  real(preal) :: &
    detfi
  integer :: &
    k, i, j, &
    instance, of
 
  li = 0.0_preal
  dli_ds  = 0.0_preal
  dli_dfi = 0.0_preal
 
  plasticitytype: select case (phase_plasticity(material_phaseat(ipc,el)))
    case (plasticity_isotropic_id) plasticitytype
      of = material_phasememberat(ipc,ip,el)
      instance = phase_plasticityinstance(material_phaseat(ipc,el))
      call plastic_isotropic_lianditstangent(my_li, my_dli_ds, s ,instance,of)
    case default plasticitytype
      my_li = 0.0_preal
      my_dli_ds = 0.0_preal
  end select plasticitytype
 
  li = li + my_li
  dli_ds = dli_ds + my_dli_ds
 
  kinematicsloop: do k = 1, phase_nkinematics(material_phaseat(ipc,el))
    kinematicstype: select case (phase_kinematics(k,material_phaseat(ipc,el)))
      case (kinematics_cleavage_opening_id) kinematicstype
        call kinematics_cleavage_opening_lianditstangent(my_li, my_dli_ds, s, ipc, ip, el)
      case (kinematics_slipplane_opening_id) kinematicstype
        call kinematics_slipplane_opening_lianditstangent(my_li, my_dli_ds, s, ipc, ip, el)
      case (kinematics_thermal_expansion_id) kinematicstype
        call kinematics_thermal_expansion_lianditstangent(my_li, my_dli_ds, ipc, ip, el)
      case default kinematicstype
        my_li = 0.0_preal
        my_dli_ds = 0.0_preal
    end select kinematicstype
    li = li + my_li
    dli_ds = dli_ds + my_dli_ds
  enddo kinematicsloop
 
  fiinv = math_inv33(fi)
  detfi = math_det33(fi)
  li = matmul(matmul(fi,li),fiinv)*detfi                                                            
  temp_33 = matmul(fiinv,li)
 
  do i = 1,3; do j = 1,3
    dli_ds(1:3,1:3,i,j)  = matmul(matmul(fi,dli_ds(1:3,1:3,i,j)),fiinv)*detfi
    dli_dfi(1:3,1:3,i,j) = dli_dfi(1:3,1:3,i,j) + li*fiinv(j,i)
    dli_dfi(1:3,i,1:3,j) = dli_dfi(1:3,i,1:3,j) + math_i3*temp_33(j,i) + li*fiinv(j,i)
  enddo; enddo
 
end subroutine constitutive_lianditstangents
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
pure function constitutive_initialfi(ipc, ip, el)
 
  integer, intent(in) :: &
    ipc, &                                                                                          !< component-ID of integration point
    ip, &                                                                                           !< integration point
    el
  real(preal), dimension(3,3) :: &
    constitutive_initialfi
  integer :: &
    k
  integer :: &
    phase, &
    homog, offset
 
  constitutive_initialfi = math_i3
  phase = material_phaseat(ipc,el)
 
  kinematicsloop: do k = 1, phase_nkinematics(phase)                                                
    kinematicstype: select case (phase_kinematics(k,phase))
      case (kinematics_thermal_expansion_id) kinematicstype
        homog = material_homogenizationat(el)
        offset = thermalmapping(homog)%p(ip,el)
        constitutive_initialfi = &
          constitutive_initialfi + kinematics_thermal_expansion_initialstrain(homog,phase,offset)
    end select kinematicstype
  enddo kinematicsloop
 
end function constitutive_initialfi
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
subroutine constitutive_sanditstangents(S, dS_dFe, dS_dFi, Fe, Fi, ipc, ip, el)
 
  integer, intent(in) :: &
    ipc, &                                                                                          !< component-ID of integration point
    ip, &                                                                                           !< integration point
    el
  real(preal),   intent(in),  dimension(3,3) :: &
    fe, &                                                                                           !< elastic deformation gradient
    fi
  real(preal),   intent(out), dimension(3,3) :: &
    s
  real(preal),   intent(out), dimension(3,3,3,3) :: &
    ds_dfe, &                                                                                       !< derivative of 2nd P-K stress with respect to elastic deformation gradient
    ds_dfi
 
  call constitutive_hooke_sanditstangents(s, ds_dfe, ds_dfi, fe, fi, ipc, ip, el)
 
 
end subroutine constitutive_sanditstangents
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
subroutine constitutive_hooke_sanditstangents(S, dS_dFe, dS_dFi, &
                                              Fe, Fi, ipc, ip, el)
 
  integer, intent(in) :: &
    ipc, &                                                                                          !< component-ID of integration point
    ip, &                                                                                           !< integration point
    el
  real(pReal),   intent(in),  dimension(3,3) :: &
    Fe, &                                                                                           !< elastic deformation gradient
    Fi
  real(pReal),   intent(out), dimension(3,3) :: &
    S
  real(pReal),   intent(out), dimension(3,3,3,3) :: &
    dS_dFe, &                                                                                       !< derivative of 2nd P-K stress with respect to elastic deformation gradient
    dS_dFi
  real(pReal), dimension(3,3) :: E
  real(pReal), dimension(3,3,3,3) :: C
  integer :: &
    ho, &                                                                                           !< homogenization
    d
  integer :: &
    i, j
 
  ho = material_homogenizationat(el)
  c = math_66tosym3333(constitutive_homogenizedc(ipc,ip,el))
 
  degradationloop: do d = 1, phase_nstiffnessdegradations(material_phaseat(ipc,el))
    degradationtype: select case(phase_stiffnessdegradation(d,material_phaseat(ipc,el)))
      case (stiffness_degradation_damage_id) degradationtype
        c = c * damage(ho)%p(damagemapping(ho)%p(ip,el))**2
    end select degradationtype
  enddo degradationloop
 
  e = 0.5_preal*(matmul(transpose(fe),fe)-math_i3)                                                  
  s = math_mul3333xx33(c,matmul(matmul(transpose(fi),e),fi))                                        
 
  do i =1, 3;do j=1,3
    ds_dfe(i,j,1:3,1:3) = matmul(fe,matmul(matmul(fi,c(i,j,1:3,1:3)),transpose(fi)))                
    ds_dfi(i,j,1:3,1:3) = 2.0_preal*matmul(matmul(e,fi),c(i,j,1:3,1:3))                             
  enddo; enddo
 
end subroutine constitutive_hooke_sanditstangents
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
subroutine constitutive_collectdotstate(S, FArray, Fi, FpArray, subdt, ipc, ip, el)
 
  integer, intent(in) :: &
    ipc, &                                                                                          !< component-ID of integration point
    ip, &                                                                                           !< integration point
    el
  real(preal),  intent(in) :: &
    subdt
  real(preal),  intent(in), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: &
    farray, &                                                                                       !< elastic deformation gradient
    fparray
  real(preal),  intent(in), dimension(3,3) :: &
    fi
  real(preal),  intent(in), dimension(3,3) :: &
    s
  real(preal),              dimension(3,3) :: &
    mp
  integer :: &
    ho, &                                                                                           !< homogenization
    tme, &                                                                                          !< thermal member position
    i, &                                                                                            !< counter in source loop
    instance, of
 
  ho = material_homogenizationat(el)
  tme = thermalmapping(ho)%p(ip,el)
  of = material_phasememberat(ipc,ip,el)
  instance = phase_plasticityinstance(material_phaseat(ipc,el))
 
  mp = matmul(matmul(transpose(fi),fi),s)
 
  plasticitytype: select case (phase_plasticity(material_phaseat(ipc,el)))
 
    case (plasticity_isotropic_id) plasticitytype
      call plastic_isotropic_dotstate    (mp,instance,of)
 
    case (plasticity_phenopowerlaw_id) plasticitytype
      call plastic_phenopowerlaw_dotstate(mp,instance,of)
 
    case (plasticity_kinehardening_id) plasticitytype
      call plastic_kinehardening_dotstate(mp,instance,of)
 
    case (plasticity_dislotwin_id) plasticitytype
      call plastic_dislotwin_dotstate    (mp,temperature(ho)%p(tme),instance,of)
 
    case (plasticity_disloucla_id) plasticitytype
      call plastic_disloucla_dotstate    (mp,temperature(ho)%p(tme),instance,of)
 
    case (plasticity_nonlocal_id) plasticitytype
      call plastic_nonlocal_dotstate     (mp,farray,fparray,temperature(ho)%p(tme),subdt, &
                                          instance,of,ip,el)
  end select plasticitytype
 
  sourceloop: do i = 1, phase_nsources(material_phaseat(ipc,el))
 
    sourcetype: select case (phase_source(i,material_phaseat(ipc,el)))
 
      case (source_damage_anisobrittle_id) sourcetype
        call source_damage_anisobrittle_dotstate (s, ipc, ip, el) 
 
      case (source_damage_isoductile_id) sourcetype
        call source_damage_isoductile_dotstate   (   ipc, ip, el)
 
      case (source_damage_anisoductile_id) sourcetype
        call source_damage_anisoductile_dotstate (   ipc, ip, el)
 
      case (source_thermal_externalheat_id) sourcetype
        call source_thermal_externalheat_dotstate(material_phaseat(ipc,el),of)
 
    end select sourcetype
 
  enddo sourceloop
 
end subroutine constitutive_collectdotstate
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
subroutine constitutive_collectdeltastate(S, Fe, Fi, ipc, ip, el)
 
  integer, intent(in) :: &
    ipc, &                                                                                          !< component-ID of integration point
    ip, &                                                                                           !< integration point
    el
  real(preal),   intent(in), dimension(3,3) :: &
    s, &                                                                                            !< 2nd Piola Kirchhoff stress
    fe, &                                                                                           !< elastic deformation gradient
    fi
  real(preal),               dimension(3,3) :: &
    mp
  integer :: &
    i, &
    instance, of
 
  mp  = matmul(matmul(transpose(fi),fi),s)
  of = material_phasememberat(ipc,ip,el)
  instance = phase_plasticityinstance(material_phaseat(ipc,el))
 
  plasticitytype: select case (phase_plasticity(material_phaseat(ipc,el)))
 
    case (plasticity_kinehardening_id) plasticitytype
      call plastic_kinehardening_deltastate(mp,instance,of)
 
    case (plasticity_nonlocal_id) plasticitytype
      call plastic_nonlocal_deltastate(mp,instance,of,ip,el)
 
  end select plasticitytype
 
  sourceloop: do i = 1, phase_nsources(material_phaseat(ipc,el))
 
     sourcetype: select case (phase_source(i,material_phaseat(ipc,el)))
 
      case (source_damage_isobrittle_id) sourcetype
        call source_damage_isobrittle_deltastate  (constitutive_homogenizedc(ipc,ip,el), fe, &
                                                   ipc, ip, el)
 
    end select sourcetype
 
  enddo sourceloop
 
end subroutine constitutive_collectdeltastate
 
 
!--------------------------------------------------------------------------------------------------
!--------------------------------------------------------------------------------------------------
subroutine constitutive_results
 
  integer :: p
  character(len=pStringLen) :: group
  do p=1,size(config_name_phase)
    group = trim('current/constituent')//'/'//trim(config_name_phase(p))
    call results_closegroup(results_addgroup(group))
 
    group = trim(group)//'/plastic'
 
    call results_closegroup(results_addgroup(group))
    select case(phase_plasticity(p))
 
      case(plasticity_isotropic_id)
        call plastic_isotropic_results(phase_plasticityinstance(p),group)
 
      case(plasticity_phenopowerlaw_id)
        call plastic_phenopowerlaw_results(phase_plasticityinstance(p),group)
 
      case(plasticity_kinehardening_id)
        call plastic_kinehardening_results(phase_plasticityinstance(p),group)
 
      case(plasticity_dislotwin_id)
        call plastic_dislotwin_results(phase_plasticityinstance(p),group)
 
      case(plasticity_disloucla_id)
        call plastic_disloucla_results(phase_plasticityinstance(p),group)
 
      case(plasticity_nonlocal_id)
        call plastic_nonlocal_results(phase_plasticityinstance(p),group)
    end select
 
  enddo
 
end subroutine constitutive_results
 
end module constitutive