DAMASK_EICMD/src/constitutive.f90

!--------------------------------------------------------------------------------------------------
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief elasticity, plasticity, damage & thermal internal microstructure state
!--------------------------------------------------------------------------------------------------
module constitutive
  use prec
  use math
  use rotations
  use IO
  use config
  use material
  use results
  use lattice
  use discretization
  use parallelization
  use HDF5_utilities
  use DAMASK_interface
  use FEsolving
  use results

  implicit none
  private
  enum, bind(c); enumerator :: &
    PLASTICITY_UNDEFINED_ID, &
    PLASTICITY_NONE_ID, &
    PLASTICITY_ISOTROPIC_ID, &
    PLASTICITY_PHENOPOWERLAW_ID, &
    PLASTICITY_KINEHARDENING_ID, &
    PLASTICITY_DISLOTWIN_ID, &
    PLASTICITY_DISLOTUNGSTEN_ID, &
    PLASTICITY_NONLOCAL_ID, &
    SOURCE_UNDEFINED_ID ,&
    SOURCE_THERMAL_DISSIPATION_ID, &
    SOURCE_THERMAL_EXTERNALHEAT_ID, &
    SOURCE_DAMAGE_ISOBRITTLE_ID, &
    SOURCE_DAMAGE_ISODUCTILE_ID, &
    SOURCE_DAMAGE_ANISOBRITTLE_ID, &
    SOURCE_DAMAGE_ANISODUCTILE_ID, &
    KINEMATICS_UNDEFINED_ID ,&
    KINEMATICS_CLEAVAGE_OPENING_ID, &
    KINEMATICS_SLIPPLANE_OPENING_ID, &
    KINEMATICS_THERMAL_EXPANSION_ID
  end enum
  real(pReal),               dimension(:,:,:),        allocatable, public :: &
    crystallite_dt                                                                                  !< requested time increment of each grain
  real(pReal),               dimension(:,:,:),        allocatable :: &
    crystallite_subdt                                                                               !< substepped time increment of each grain
  type(rotation),            dimension(:,:,:),        allocatable :: &
    crystallite_orientation                                                                         !< current orientation
  real(pReal),               dimension(:,:,:,:,:),    allocatable :: &
    crystallite_F0, &                                                                               !< def grad at start of FE inc
    crystallite_subF,  &                                                                            !< def grad to be reached at end of crystallite inc
    crystallite_subF0, &                                                                            !< def grad at start of crystallite inc
    crystallite_Fe, &                                                                               !< current "elastic" def grad (end of converged time step)
    crystallite_subFp0,&                                                                            !< plastic def grad at start of crystallite inc
    crystallite_subFi0,&                                                                            !< intermediate def grad at start of crystallite inc
    crystallite_Lp0, &                                                                              !< plastic velocitiy grad at start of FE inc
    crystallite_partitionedLp0, &                                                                   !< plastic velocity grad at start of homog inc
    crystallite_S0, &                                                                               !< 2nd Piola-Kirchhoff stress vector at start of FE inc
    crystallite_partitionedS0                                                                       !< 2nd Piola-Kirchhoff stress vector at start of homog inc
  real(pReal),               dimension(:,:,:,:,:),    allocatable, public :: &
    crystallite_P, &                                                                                !< 1st Piola-Kirchhoff stress per grain
    crystallite_Lp, &                                                                               !< current plastic velocitiy grad (end of converged time step)
    crystallite_S, &                                                                                !< current 2nd Piola-Kirchhoff stress vector (end of converged time step)
    crystallite_partitionedF0                                                                       !< def grad at start of homog inc
  real(pReal),               dimension(:,:,:,:,:),    allocatable, public :: &
    crystallite_partitionedF                                                                        !< def grad to be reached at end of homog inc

  logical,                    dimension(:,:,:),       allocatable :: &
    crystallite_converged                                                                           !< convergence flag


  type :: tTensorContainer
    real(pReal), dimension(:,:,:), allocatable :: data
  end type

  type(tTensorContainer), dimension(:), allocatable :: &
    constitutive_mech_Fi, &
    constitutive_mech_Fi0, &
    constitutive_mech_partionedFi0, &
    constitutive_mech_Li, &
    constitutive_mech_Li0, &
    constitutive_mech_partionedLi0, &
    constitutive_mech_Fp, &
    constitutive_mech_Fp0, &
    constitutive_mech_partionedFp0


  type :: tNumerics
    integer :: &
      iJacoLpresiduum, &                                                                            !< frequency of Jacobian update of residuum in Lp
      nState, &                                                                                     !< state loop limit
      nStress                                                                                       !< stress loop limit
    real(pReal) :: &
      subStepMinCryst, &                                                                            !< minimum (relative) size of sub-step allowed during cutback
      subStepSizeCryst, &                                                                           !< size of first substep when cutback
      subStepSizeLp, &                                                                              !< size of first substep when cutback in Lp calculation
      subStepSizeLi, &                                                                              !< size of first substep when cutback in Li calculation
      stepIncreaseCryst, &                                                                          !< increase of next substep size when previous substep converged
      rtol_crystalliteState, &                                                                      !< relative tolerance in state loop
      rtol_crystalliteStress, &                                                                     !< relative tolerance in stress loop
      atol_crystalliteStress                                                                        !< absolute tolerance in stress loop
  end type tNumerics

  type(tNumerics) :: num                                                                            ! numerics parameters. Better name?

  type :: tDebugOptions
    logical :: &
      basic, &
      extensive, &
      selective
    integer :: &
      element, &
      ip, &
      grain
  end type tDebugOptions

  type(tDebugOptions) :: debugCrystallite

  procedure(integrateStateFPI), pointer :: integrateState

  integer(kind(PLASTICITY_undefined_ID)), dimension(:),   allocatable :: &
    phase_plasticity                                                                                !< plasticity of each phase

  integer(kind(SOURCE_undefined_ID)),     dimension(:,:), allocatable :: &
    phase_source, &                                                                                 !< active sources mechanisms of each phase
    phase_kinematics                                                                                !< active kinematic mechanisms of each phase

  integer, dimension(:), allocatable, public :: &                                                   !< ToDo: should be protected (bug in Intel compiler)
    phase_Nsources, &                                                                               !< number of source mechanisms active in each phase
    phase_Nkinematics, &                                                                            !< number of kinematic mechanisms active in each phase
    phase_NstiffnessDegradations, &                                                                 !< number of stiffness degradation mechanisms active in each phase
    phase_plasticityInstance, &                                                                     !< instance of particular plasticity of each phase
    phase_elasticityInstance                                                                        !< instance of particular elasticity of each phase

  logical, dimension(:), allocatable, public :: &                                                   ! ToDo: should be protected (bug in Intel Compiler)
    phase_localPlasticity                                                                           !< flags phases with local constitutive law

  type(tPlasticState), allocatable, dimension(:), public :: &
    plasticState
  type(tSourceState),  allocatable, dimension(:), public :: &
    sourceState


  integer, public, protected :: &
    constitutive_plasticity_maxSizeDotState, &
    constitutive_source_maxSizeDotState

  interface

! == cleaned:begin =================================================================================
    module subroutine mech_init
    end subroutine mech_init

    module subroutine damage_init
    end subroutine damage_init

    module subroutine thermal_init
    end subroutine thermal_init


    module subroutine mech_results(group,ph)
      character(len=*), intent(in) :: group
      integer,          intent(in) :: ph
    end subroutine mech_results

    module subroutine damage_results(group,ph)
      character(len=*), intent(in) :: group
      integer,          intent(in) :: ph
    end subroutine damage_results


    module subroutine mech_restart_read(fileHandle)
      integer(HID_T), intent(in) :: fileHandle
    end subroutine mech_restart_read

    module subroutine mech_initializeRestorationPoints(ph,me)
      integer, intent(in) :: ph, me
    end subroutine mech_initializeRestorationPoints

    module subroutine constitutive_mech_windForward(ph,me)
      integer, intent(in) :: ph, me
    end subroutine constitutive_mech_windForward

    module subroutine constitutive_mech_forward
    end subroutine constitutive_mech_forward

! == cleaned:end ===================================================================================

    module subroutine source_damage_anisoBrittle_dotState(S, co, ip, el)
      integer, intent(in) :: &
        co, &                                                                                      !< component-ID of integration point
        ip, &                                                                                       !< integration point
        el                                                                                          !< element
      real(pReal),  intent(in), dimension(3,3) :: &
        S
    end subroutine source_damage_anisoBrittle_dotState

    module subroutine source_damage_anisoDuctile_dotState(co, ip, el)
      integer, intent(in) :: &
        co, &                                                                                      !< component-ID of integration point
        ip, &                                                                                       !< integration point
        el                                                                                          !< element
    end subroutine source_damage_anisoDuctile_dotState

    module subroutine source_damage_isoDuctile_dotState(co, ip, el)
      integer, intent(in) :: &
        co, &                                                                                      !< component-ID of integration point
        ip, &                                                                                       !< integration point
        el                                                                                          !< element
    end subroutine source_damage_isoDuctile_dotState

    module subroutine source_thermal_externalheat_dotState(phase, of)
      integer, intent(in) :: &
        phase, &
        of
    end subroutine source_thermal_externalheat_dotState

    module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el)
      integer, intent(in) :: &
        ip, &                                                                                       !< integration point number
        el                                                                                          !< element number
      real(pReal), intent(in) :: &
        phi                                                                                         !< damage parameter
      real(pReal), intent(inout) :: &
        phiDot, &
        dPhiDot_dPhi
    end subroutine constitutive_damage_getRateAndItsTangents

    module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, S, Lp, ip, el)
      integer, intent(in) :: &
        ip, &                                                                                       !< integration point number
        el                                                                                          !< element number
      real(pReal), intent(in) :: &
        T
      real(pReal),  intent(in), dimension(:,:,:,:,:) :: &
        S, &                                                                                        !< current 2nd Piola Kitchoff stress vector
        Lp                                                                                          !< plastic velocity gradient
      real(pReal), intent(inout) :: &
        TDot, &
        dTDot_dT
    end subroutine constitutive_thermal_getRateAndItsTangents

    module function plastic_dislotwin_homogenizedC(co,ip,el) result(homogenizedC)
      real(pReal), dimension(6,6) :: &
        homogenizedC
      integer,     intent(in) :: &
        co, &                                                                                      !< component-ID of integration point
        ip, &                                                                                       !< integration point
        el                                                                                          !< element
    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_nIPs,discretization_Nelems), intent(in) :: &
        orientation                                                                                 !< crystal orientation
    end subroutine plastic_nonlocal_updateCompatibility

    module subroutine plastic_isotropic_LiAndItsTangent(Li,dLi_dMi,Mi,instance,of)
      real(pReal), dimension(3,3),     intent(out) :: &
        Li                                                                                          !< inleastic velocity gradient
      real(pReal), dimension(3,3,3,3), intent(out)  :: &
        dLi_dMi                                                                                     !< derivative of Li with respect to Mandel stress
      real(pReal), dimension(3,3),     intent(in) :: &
        Mi                                                                                          !< Mandel stress
      integer,                         intent(in) :: &
        instance, &
        of
    end subroutine plastic_isotropic_LiAndItsTangent

    module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, co, ip, el)
      integer, intent(in) :: &
        co, &                                                                                      !< grain number
        ip, &                                                                                       !< integration point number
        el                                                                                          !< element number
      real(pReal),   intent(in),  dimension(3,3) :: &
        S
      real(pReal),   intent(out), dimension(3,3) :: &
        Ld                                                                                          !< damage velocity gradient
      real(pReal),   intent(out), dimension(3,3,3,3) :: &
        dLd_dTstar                                                                                  !< derivative of Ld with respect to Tstar (4th-order tensor)
    end subroutine kinematics_cleavage_opening_LiAndItsTangent

    module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, co, ip, el)
      integer, intent(in) :: &
        co, &                                                                                      !< grain number
        ip, &                                                                                       !< integration point number
        el                                                                                          !< element number
      real(pReal),   intent(in),  dimension(3,3) :: &
        S
      real(pReal),   intent(out), dimension(3,3) :: &
        Ld                                                                                          !< damage velocity gradient
      real(pReal),   intent(out), dimension(3,3,3,3) :: &
        dLd_dTstar                                                                                  !< derivative of Ld with respect to Tstar (4th-order tensor)
    end subroutine kinematics_slipplane_opening_LiAndItsTangent

    module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, co, ip, el)
      integer, intent(in) :: &
        co, &                                                                                      !< 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_dTstar                                                                                  !< derivative of Li with respect to Tstar (4th-order tensor defined to be zero)
    end subroutine kinematics_thermal_expansion_LiAndItsTangent


    module subroutine source_damage_isoBrittle_deltaState(C, Fe, co, ip, el)
      integer, intent(in) :: &
        co, &                                                                                      !< 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
    end subroutine source_damage_isoBrittle_deltaState


    module subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
                                         S, Fi, co, ip, el)
      integer, intent(in) :: &
        co, &                                                                                      !< component-ID of integration point
        ip, &                                                                                       !< integration point
        el                                                                                          !< element
      real(pReal),   intent(in),  dimension(3,3) :: &
        S, &                                                                                        !< 2nd Piola-Kirchhoff stress
        Fi                                                                                          !< intermediate deformation gradient
      real(pReal),   intent(out), dimension(3,3) :: &
        Lp                                                                                          !< plastic velocity gradient
      real(pReal),   intent(out), dimension(3,3,3,3) :: &
        dLp_dS, &
        dLp_dFi                                                                                     !< derivative of Lp with respect to Fi
    end subroutine constitutive_plastic_LpAndItsTangents


    module subroutine constitutive_plastic_dependentState(F, co, ip, el)
      integer, intent(in) :: &
        co, &                                                                                      !< component-ID of integration point
        ip, &                                                                                       !< integration point
        el                                                                                          !< element
      real(pReal),   intent(in), dimension(3,3) :: &
        F                                                                                           !< elastic deformation gradient
    end subroutine constitutive_plastic_dependentState



    module subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, Fe, Fi, co, ip, el)
      integer, intent(in) :: &
        co, &                                                                                      !< component-ID of integration point
        ip, &                                                                                       !< integration point
        el                                                                                          !< element
      real(pReal),   intent(in),  dimension(3,3) :: &
        Fe, &                                                                                       !< elastic deformation gradient
        Fi                                                                                          !< intermediate deformation gradient
      real(pReal),   intent(out), dimension(3,3) :: &
        S                                                                                           !< 2nd Piola-Kirchhoff stress tensor
      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                                                                                      !< derivative of 2nd P-K stress with respect to intermediate deformation gradient
    end subroutine constitutive_hooke_SandItsTangents

    module subroutine integrateStateFPI(co,ip,el)
      integer, intent(in) :: co, ip, el
    end subroutine integrateStateFPI

  end interface



  type(tDebugOptions) :: debugConstitutive

  public :: &
    constitutive_init, &
    constitutive_homogenizedC, &
    constitutive_LiAndItsTangents, &
    constitutive_damage_getRateAndItsTangents, &
    constitutive_thermal_getRateAndItsTangents, &
    constitutive_results, &
    constitutive_allocateState, &
    constitutive_forward, &
    constitutive_restore, &
    plastic_nonlocal_updateCompatibility, &
    source_active, &
    kinematics_active, &
    converged, &
    crystallite_init, &
    crystallite_stress, &
    crystallite_stressTangent, &
    crystallite_orientations, &
    crystallite_push33ToRef, &
    crystallite_restartWrite, &
    crystallite_restartRead, &
    constitutive_initializeRestorationPoints, &
    constitutive_windForward, &
    crystallite_restore

contains


!--------------------------------------------------------------------------------------------------
!> @brief Initialze constitutive models for individual physics
!--------------------------------------------------------------------------------------------------
subroutine constitutive_init

  integer :: &
    p, &                                                                                            !< counter in phase loop
    s                                                                                               !< counter in source loop
  class (tNode), pointer :: &
    debug_constitutive, &
    phases

  debug_constitutive => config_debug%get('constitutive', defaultVal=emptyList)
  debugConstitutive%basic      =  debug_constitutive%contains('basic')
  debugConstitutive%extensive  =  debug_constitutive%contains('extensive')
  debugConstitutive%selective  =  debug_constitutive%contains('selective')
  debugConstitutive%element    =  config_debug%get_asInt('element',defaultVal = 1)
  debugConstitutive%ip         =  config_debug%get_asInt('integrationpoint',defaultVal = 1)
  debugConstitutive%grain      =  config_debug%get_asInt('grain',defaultVal = 1)

!--------------------------------------------------------------------------------------------------
! initialize constitutive laws
  call mech_init
  call damage_init
  call thermal_init

  print'(/,a)', ' <<<+-  constitutive init  -+>>>'; flush(IO_STDOUT)

  phases => config_material%get('phase')
  constitutive_source_maxSizeDotState = 0
  PhaseLoop2:do p = 1,phases%length
!--------------------------------------------------------------------------------------------------
! partition and initialize state
    plasticState(p)%partitionedState0 = plasticState(p)%state0
    plasticState(p)%state             = plasticState(p)%partitionedState0
    forall(s = 1:phase_Nsources(p))
      sourceState(p)%p(s)%partitionedState0 = sourceState(p)%p(s)%state0
      sourceState(p)%p(s)%state             = sourceState(p)%p(s)%partitionedState0
    end forall

    constitutive_source_maxSizeDotState   = max(constitutive_source_maxSizeDotState, &
                                                maxval(sourceState(p)%p%sizeDotState))
  enddo PhaseLoop2
  constitutive_plasticity_maxSizeDotState = maxval(plasticState%sizeDotState)

end subroutine constitutive_init


!--------------------------------------------------------------------------------------------------
!> @brief checks if a source mechanism is active or not
!--------------------------------------------------------------------------------------------------
function source_active(source_label,src_length)  result(active_source)

  character(len=*), intent(in)         :: source_label                                              !< name of source mechanism
  integer,          intent(in)         :: src_length                                                !< max. number of sources in system
  logical, dimension(:,:), allocatable :: active_source

  class(tNode), pointer :: &
    phases, &
    phase, &
    sources, &
    src
  integer :: p,s

  phases => config_material%get('phase')
  allocate(active_source(src_length,phases%length), source = .false. )
  do p = 1, phases%length
    phase => phases%get(p)
    sources => phase%get('source',defaultVal=emptyList)
    do s = 1, sources%length
      src => sources%get(s)
      if(src%get_asString('type') == source_label) active_source(s,p) = .true.
    enddo
  enddo


end function source_active


!--------------------------------------------------------------------------------------------------
!> @brief checks if a kinematic mechanism is active or not
!--------------------------------------------------------------------------------------------------
function kinematics_active(kinematics_label,kinematics_length)  result(active_kinematics)

  character(len=*), intent(in)         :: kinematics_label                                          !< name of kinematic mechanism
  integer,          intent(in)         :: kinematics_length                                         !< max. number of kinematics in system
  logical, dimension(:,:), allocatable :: active_kinematics

  class(tNode), pointer :: &
    phases, &
    phase, &
    kinematics, &
    kinematics_type
  integer :: p,k

  phases => config_material%get('phase')
  allocate(active_kinematics(kinematics_length,phases%length), source = .false. )
  do p = 1, phases%length
    phase => phases%get(p)
    kinematics => phase%get('kinematics',defaultVal=emptyList)
    do k = 1, kinematics%length
      kinematics_type => kinematics%get(k)
      if(kinematics_type%get_asString('type') == kinematics_label) active_kinematics(k,p) = .true.
    enddo
  enddo


end function kinematics_active


!--------------------------------------------------------------------------------------------------
!> @brief returns the homogenize elasticity matrix
!> ToDo: homogenizedC66 would be more consistent
!--------------------------------------------------------------------------------------------------
function constitutive_homogenizedC(co,ip,el)

  real(pReal), dimension(6,6) :: &
    constitutive_homogenizedC
  integer,      intent(in)     :: &
    co, &                                                                                          !< component-ID of integration point
    ip, &                                                                                           !< integration point
    el                                                                                              !< element

  plasticityType: select case (phase_plasticity(material_phaseAt(co,el)))
    case (PLASTICITY_DISLOTWIN_ID) plasticityType
     constitutive_homogenizedC = plastic_dislotwin_homogenizedC(co,ip,el)
    case default plasticityType
     constitutive_homogenizedC = lattice_C66(1:6,1:6,material_phaseAt(co,el))
  end select plasticityType

end function constitutive_homogenizedC


!--------------------------------------------------------------------------------------------------
!> @brief  contains the constitutive equation for calculating the velocity gradient
! ToDo: MD: S is Mi?
!--------------------------------------------------------------------------------------------------
subroutine constitutive_LiAndItsTangents(Li, dLi_dS, dLi_dFi, &
                                         S, Fi, co, ip, el)

  integer, intent(in) :: &
    co, &                                                                                          !< component-ID of integration point
    ip, &                                                                                           !< integration point
    el                                                                                              !< element
  real(pReal),   intent(in),  dimension(3,3) :: &
    S                                                                                               !< 2nd Piola-Kirchhoff stress
  real(pReal),   intent(in),  dimension(3,3) :: &
    Fi                                                                                              !< intermediate deformation gradient
  real(pReal),   intent(out), dimension(3,3) :: &
    Li                                                                                              !< intermediate velocity gradient
  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(co,el)))
    case (PLASTICITY_isotropic_ID) plasticityType
      of = material_phasememberAt(co,ip,el)
      instance = phase_plasticityInstance(material_phaseAt(co,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(co,el))
    kinematicsType: select case (phase_kinematics(k,material_phaseAt(co,el)))
      case (KINEMATICS_cleavage_opening_ID) kinematicsType
        call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el)
      case (KINEMATICS_slipplane_opening_ID) kinematicsType
        call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, co, ip, el)
      case (KINEMATICS_thermal_expansion_ID) kinematicsType
        call kinematics_thermal_expansion_LiAndItsTangent(my_Li, my_dLi_dS, co, 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                                                            !< push forward to intermediate configuration
  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


!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
function constitutive_damage_collectDotState(S, co, ip, el,phase,of) result(broken)

  integer, intent(in) :: &
    co, &                                                                                          !< component-ID of integration point
    ip, &                                                                                           !< integration point
    el, &                                                                                           !< element
    phase, &
    of
  real(pReal),  intent(in), dimension(3,3) :: &
    S                                                                                               !< 2nd Piola Kirchhoff stress (vector notation)
  integer :: &
    i                                                                                               !< counter in source loop
  logical :: broken


  broken = .false.

  SourceLoop: do i = 1, phase_Nsources(phase)

    sourceType: select case (phase_source(i,phase))

      case (SOURCE_damage_anisoBrittle_ID) sourceType
        call source_damage_anisoBrittle_dotState(S, co, ip, el) ! correct stress?

      case (SOURCE_damage_isoDuctile_ID) sourceType
        call source_damage_isoDuctile_dotState(co, ip, el)

      case (SOURCE_damage_anisoDuctile_ID) sourceType
        call source_damage_anisoDuctile_dotState(co, ip, el)

    end select sourceType

    broken = broken .or. any(IEEE_is_NaN(sourceState(phase)%p(i)%dotState(:,of)))

  enddo SourceLoop

end function constitutive_damage_collectDotState


!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
function constitutive_thermal_collectDotState(ph,me) result(broken)

  integer, intent(in) :: ph, me
  logical :: broken

  integer :: i


  broken = .false.

  SourceLoop: do i = 1, phase_Nsources(ph)

    if (phase_source(i,ph) == SOURCE_thermal_externalheat_ID) &
      call source_thermal_externalheat_dotState(ph,me)

    broken = broken .or. any(IEEE_is_NaN(sourceState(ph)%p(i)%dotState(:,me)))

  enddo SourceLoop

end function constitutive_thermal_collectDotState


!--------------------------------------------------------------------------------------------------
!> @brief for constitutive models having an instantaneous change of state
!> will return false if delta state is not needed/supported by the constitutive model
!--------------------------------------------------------------------------------------------------
function constitutive_damage_deltaState(Fe, co, ip, el, ph, of) result(broken)

  integer, intent(in) :: &
    co, &                                                                                          !< component-ID of integration point
    ip, &                                                                                           !< integration point
    el, &                                                                                           !< element
    ph, &
    of
  real(pReal),   intent(in), dimension(3,3) :: &
    Fe                                                                                              !< elastic deformation gradient
  integer :: &
    so, &
    myOffset, &
    mySize
  logical :: &
    broken


  broken = .false.

  sourceLoop: do so = 1, phase_Nsources(ph)

     sourceType: select case (phase_source(so,ph))

      case (SOURCE_damage_isoBrittle_ID) sourceType
        call source_damage_isoBrittle_deltaState  (constitutive_homogenizedC(co,ip,el), Fe, &
                                                   co, ip, el)
        broken = any(IEEE_is_NaN(sourceState(ph)%p(so)%deltaState(:,of)))
        if(.not. broken) then
          myOffset = sourceState(ph)%p(so)%offsetDeltaState
          mySize   = sourceState(ph)%p(so)%sizeDeltaState
          sourceState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,of) = &
          sourceState(ph)%p(so)%state(myOffset + 1: myOffset + mySize,of) + sourceState(ph)%p(so)%deltaState(1:mySize,of)
        endif

    end select sourceType

  enddo SourceLoop

end function constitutive_damage_deltaState


!--------------------------------------------------------------------------------------------------
!> @brief Allocate the components of the state structure for a given phase
!--------------------------------------------------------------------------------------------------
subroutine constitutive_allocateState(state, &
                                  Nconstituents,sizeState,sizeDotState,sizeDeltaState)

  class(tState), intent(out) :: &
    state
  integer, intent(in) :: &
    Nconstituents, &
    sizeState, &
    sizeDotState, &
    sizeDeltaState

  state%sizeState        = sizeState
  state%sizeDotState     = sizeDotState
  state%sizeDeltaState   = sizeDeltaState
  state%offsetDeltaState = sizeState-sizeDeltaState                                                 ! deltaState occupies latter part of state by definition

  allocate(state%atol           (sizeState),             source=0.0_pReal)
  allocate(state%state0         (sizeState,Nconstituents), source=0.0_pReal)
  allocate(state%partitionedState0(sizeState,Nconstituents), source=0.0_pReal)
  allocate(state%subState0      (sizeState,Nconstituents), source=0.0_pReal)
  allocate(state%state          (sizeState,Nconstituents), source=0.0_pReal)

  allocate(state%dotState    (sizeDotState,Nconstituents), source=0.0_pReal)

  allocate(state%deltaState(sizeDeltaState,Nconstituents), source=0.0_pReal)


end subroutine constitutive_allocateState


!--------------------------------------------------------------------------------------------------
!> @brief Restore data after homog cutback.
!--------------------------------------------------------------------------------------------------
subroutine constitutive_restore(ip,el)

  integer, intent(in) :: &
    ip, &                                                                                            !< integration point number
    el                                                                                               !< element number
  integer :: &
    co, &                                                                                            !< constituent number
    s

  do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
    do s = 1, phase_Nsources(material_phaseAt(co,el))
      sourceState(material_phaseAt(co,el))%p(s)%state(          :,material_phasememberAt(co,ip,el)) = &
      sourceState(material_phaseAt(co,el))%p(s)%partitionedState0(:,material_phasememberAt(co,ip,el))
    enddo
  enddo

end subroutine constitutive_restore


!--------------------------------------------------------------------------------------------------
!> @brief Forward data after successful increment.
! ToDo: Any guessing for the current states possible?
!--------------------------------------------------------------------------------------------------
subroutine constitutive_forward

  integer :: i, j

  crystallite_F0  = crystallite_partitionedF
  crystallite_Lp0 = crystallite_Lp
  crystallite_S0  = crystallite_S

  call constitutive_mech_forward()

  do i = 1, size(sourceState)
    do j = 1,phase_Nsources(i)
      sourceState(i)%p(j)%state0 = sourceState(i)%p(j)%state
  enddo; enddo

end subroutine constitutive_forward


!--------------------------------------------------------------------------------------------------
!> @brief writes constitutive results to HDF5 output file
!--------------------------------------------------------------------------------------------------
subroutine constitutive_results

  integer :: ph
  character(len=:), allocatable :: group


  call results_closeGroup(results_addGroup('/current/phase/'))

  do ph = 1, size(material_name_phase)

    group = '/current/phase/'//trim(material_name_phase(ph))//'/'
    call results_closeGroup(results_addGroup(group))

    call mech_results(group,ph)
    call damage_results(group,ph)

  enddo

end subroutine constitutive_results


!--------------------------------------------------------------------------------------------------
!> @brief allocates and initialize per grain variables
!--------------------------------------------------------------------------------------------------
subroutine crystallite_init

  integer :: &
    Nconstituents, &
    ph, &
    me, &
    co, &                                                                                            !< counter in integration point component loop
    ip, &                                                                                            !< counter in integration point loop
    el, &                                                                                            !< counter in element loop
    cMax, &                                                                                         !< maximum number of  integration point components
    iMax, &                                                                                         !< maximum number of integration points
    eMax                                                                                            !< maximum number of elements

  class(tNode), pointer :: &
    num_crystallite, &
    debug_crystallite, &                                                                            ! pointer to debug options for crystallite
    phases, &
    phase, &
    mech


  print'(/,a)', ' <<<+-  crystallite init  -+>>>'

  debug_crystallite => config_debug%get('crystallite', defaultVal=emptyList)
  debugCrystallite%extensive = debug_crystallite%contains('extensive')

  cMax = homogenization_maxNconstituents
  iMax = discretization_nIPs
  eMax = discretization_Nelems

  allocate(crystallite_partitionedF(3,3,cMax,iMax,eMax),source=0.0_pReal)

  allocate(crystallite_S0, &
           crystallite_F0,crystallite_Lp0, &
           crystallite_partitionedS0, &
           crystallite_partitionedF0,&
                                   crystallite_partitionedLp0, &
           crystallite_S,crystallite_P, &
           crystallite_Fe,crystallite_Lp, &
           crystallite_subF,crystallite_subF0, &
           crystallite_subFp0,crystallite_subFi0, &
           source = crystallite_partitionedF)

  allocate(crystallite_dt(cMax,iMax,eMax),source=0.0_pReal)
  allocate(crystallite_subdt, &
           source = crystallite_dt)

  allocate(crystallite_orientation(cMax,iMax,eMax))

  allocate(crystallite_converged(cMax,iMax,eMax),             source=.true.)

  num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict)

  num%subStepMinCryst        = num_crystallite%get_asFloat ('subStepMin',       defaultVal=1.0e-3_pReal)
  num%subStepSizeCryst       = num_crystallite%get_asFloat ('subStepSize',      defaultVal=0.25_pReal)
  num%stepIncreaseCryst      = num_crystallite%get_asFloat ('stepIncrease',     defaultVal=1.5_pReal)
  num%subStepSizeLp          = num_crystallite%get_asFloat ('subStepSizeLp',    defaultVal=0.5_pReal)
  num%subStepSizeLi          = num_crystallite%get_asFloat ('subStepSizeLi',    defaultVal=0.5_pReal)
  num%rtol_crystalliteState  = num_crystallite%get_asFloat ('rtol_State',       defaultVal=1.0e-6_pReal)
  num%rtol_crystalliteStress = num_crystallite%get_asFloat ('rtol_Stress',      defaultVal=1.0e-6_pReal)
  num%atol_crystalliteStress = num_crystallite%get_asFloat ('atol_Stress',      defaultVal=1.0e-8_pReal)
  num%iJacoLpresiduum        = num_crystallite%get_asInt   ('iJacoLpresiduum',  defaultVal=1)
  num%nState                 = num_crystallite%get_asInt   ('nState',           defaultVal=20)
  num%nStress                = num_crystallite%get_asInt   ('nStress',          defaultVal=40)

  if(num%subStepMinCryst   <= 0.0_pReal)      call IO_error(301,ext_msg='subStepMinCryst')
  if(num%subStepSizeCryst  <= 0.0_pReal)      call IO_error(301,ext_msg='subStepSizeCryst')
  if(num%stepIncreaseCryst <= 0.0_pReal)      call IO_error(301,ext_msg='stepIncreaseCryst')

  if(num%subStepSizeLp <= 0.0_pReal)          call IO_error(301,ext_msg='subStepSizeLp')
  if(num%subStepSizeLi <= 0.0_pReal)          call IO_error(301,ext_msg='subStepSizeLi')

  if(num%rtol_crystalliteState  <= 0.0_pReal) call IO_error(301,ext_msg='rtol_crystalliteState')
  if(num%rtol_crystalliteStress <= 0.0_pReal) call IO_error(301,ext_msg='rtol_crystalliteStress')
  if(num%atol_crystalliteStress <= 0.0_pReal) call IO_error(301,ext_msg='atol_crystalliteStress')

  if(num%iJacoLpresiduum < 1)                 call IO_error(301,ext_msg='iJacoLpresiduum')

  if(num%nState < 1)                          call IO_error(301,ext_msg='nState')
  if(num%nStress< 1)                          call IO_error(301,ext_msg='nStress')


  phases => config_material%get('phase')

  allocate(constitutive_mech_Fi(phases%length))
  allocate(constitutive_mech_Fi0(phases%length))
  allocate(constitutive_mech_partionedFi0(phases%length))
  allocate(constitutive_mech_Fp(phases%length))
  allocate(constitutive_mech_Fp0(phases%length))
  allocate(constitutive_mech_partionedFp0(phases%length))
  allocate(constitutive_mech_Li(phases%length))
  allocate(constitutive_mech_Li0(phases%length))
  allocate(constitutive_mech_partionedLi0(phases%length))
  do ph = 1, phases%length
    Nconstituents = count(material_phaseAt == ph) * discretization_nIPs

    allocate(constitutive_mech_Fi(ph)%data(3,3,Nconstituents))
    allocate(constitutive_mech_Fi0(ph)%data(3,3,Nconstituents))
    allocate(constitutive_mech_partionedFi0(ph)%data(3,3,Nconstituents))
    allocate(constitutive_mech_Fp(ph)%data(3,3,Nconstituents))
    allocate(constitutive_mech_Fp0(ph)%data(3,3,Nconstituents))
    allocate(constitutive_mech_partionedFp0(ph)%data(3,3,Nconstituents))
    allocate(constitutive_mech_Li(ph)%data(3,3,Nconstituents))
    allocate(constitutive_mech_Li0(ph)%data(3,3,Nconstituents))
    allocate(constitutive_mech_partionedLi0(ph)%data(3,3,Nconstituents))
  enddo

  print'(a42,1x,i10)', '    # of elements:                       ', eMax
  print'(a42,1x,i10)', '    # of integration points/element:     ', iMax
  print'(a42,1x,i10)', 'max # of constituents/integration point: ', cMax
  flush(IO_STDOUT)

 !$OMP PARALLEL DO PRIVATE(ph,me)
  do el = FEsolving_execElem(1),FEsolving_execElem(2)
    do ip = FEsolving_execIP(1), FEsolving_execIP(2); do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))

      ph = material_phaseAt(co,el)
      me = material_phaseMemberAt(co,ip,el)
      constitutive_mech_Fp0(ph)%data(1:3,1:3,me) = material_orientation0(co,ip,el)%asMatrix()                      ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005)
      constitutive_mech_Fp0(ph)%data(1:3,1:3,me) = constitutive_mech_Fp0(ph)%data(1:3,1:3,me) &
                                     / math_det33(constitutive_mech_Fp0(ph)%data(1:3,1:3,me))**(1.0_pReal/3.0_pReal)
      constitutive_mech_Fi0(ph)%data(1:3,1:3,me) = math_I3

      crystallite_F0(1:3,1:3,co,ip,el)  = math_I3

      crystallite_Fe(1:3,1:3,co,ip,el)  = math_inv33(matmul(constitutive_mech_Fi0(ph)%data(1:3,1:3,me), &
                                                         constitutive_mech_Fp0(ph)%data(1:3,1:3,me)))           ! assuming that euler angles are given in internal strain free configuration
      constitutive_mech_Fp(ph)%data(1:3,1:3,me)  = constitutive_mech_Fp0(ph)%data(1:3,1:3,me)
      constitutive_mech_Fi(ph)%data(1:3,1:3,me) = constitutive_mech_Fi0(ph)%data(1:3,1:3,me)

      constitutive_mech_partionedFi0(ph)%data(1:3,1:3,me) = constitutive_mech_Fi0(ph)%data(1:3,1:3,me)
      constitutive_mech_partionedFp0(ph)%data(1:3,1:3,me) = constitutive_mech_Fp0(ph)%data(1:3,1:3,me)

    enddo; enddo
  enddo
  !$OMP END PARALLEL DO

  crystallite_partitionedF0  = crystallite_F0
  crystallite_partitionedF   = crystallite_F0

  call crystallite_orientations()

  !$OMP PARALLEL DO PRIVATE(ph,me)
  do el = FEsolving_execElem(1),FEsolving_execElem(2)
    do ip = FEsolving_execIP(1),FEsolving_execIP(2)
      do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
        ph = material_phaseAt(co,el)
        me = material_phaseMemberAt(co,ip,el)
        call constitutive_plastic_dependentState(crystallite_partitionedF0(1:3,1:3,co,ip,el), &
                                         co,ip,el)                                                     ! update dependent state variables to be consistent with basic states
     enddo
    enddo
  enddo
  !$OMP END PARALLEL DO


end subroutine crystallite_init


!--------------------------------------------------------------------------------------------------
!> @brief calculate stress (P)
!--------------------------------------------------------------------------------------------------
function crystallite_stress(co,ip,el)

  integer, intent(in) :: &
    co, &
    ip, &
    el

  logical :: crystallite_stress

  real(pReal) :: &
    formerSubStep
  integer :: &
    NiterationCrystallite, &                                                                        ! number of iterations in crystallite loop
    s, ph, me
  logical :: todo
  real(pReal) :: subFrac,subStep
  real(pReal), dimension(3,3) :: &
    subLp0, &                                                                                       !< plastic velocity grad at start of crystallite inc
    subLi0                                                                                          !< intermediate velocity grad at start of crystallite inc


  ph = material_phaseAt(co,el)
  me = material_phaseMemberAt(co,ip,el)
  subLi0 = constitutive_mech_partionedLi0(ph)%data(1:3,1:3,me)
  subLp0 = crystallite_partitionedLp0(1:3,1:3,co,ip,el)
  plasticState    (material_phaseAt(co,el))%subState0(      :,material_phaseMemberAt(co,ip,el)) = &
  plasticState    (material_phaseAt(co,el))%partitionedState0(:,material_phaseMemberAt(co,ip,el))

  do s = 1, phase_Nsources(material_phaseAt(co,el))
    sourceState(material_phaseAt(co,el))%p(s)%subState0(      :,material_phaseMemberAt(co,ip,el)) = &
    sourceState(material_phaseAt(co,el))%p(s)%partitionedState0(:,material_phaseMemberAt(co,ip,el))
  enddo
  crystallite_subFp0(1:3,1:3,co,ip,el) = constitutive_mech_partionedFp0(ph)%data(1:3,1:3,me)
  crystallite_subFi0(1:3,1:3,co,ip,el) = constitutive_mech_partionedFi0(ph)%data(1:3,1:3,me)
  crystallite_subF0(1:3,1:3,co,ip,el)  = crystallite_partitionedF0(1:3,1:3,co,ip,el)
  subFrac = 0.0_pReal
  subStep = 1.0_pReal/num%subStepSizeCryst
  todo = .true.
  crystallite_converged(co,ip,el) = .false.                                                      ! pretend failed step of 1/subStepSizeCryst

  todo = .true.
  NiterationCrystallite = 0
  cutbackLooping: do while (todo)
    NiterationCrystallite = NiterationCrystallite + 1

!--------------------------------------------------------------------------------------------------
!  wind forward
   if (crystallite_converged(co,ip,el)) then
     formerSubStep = subStep
     subFrac = subFrac + subStep
     subStep = min(1.0_pReal - subFrac, num%stepIncreaseCryst * subStep)

     todo = subStep > 0.0_pReal                        ! still time left to integrate on?

     if (todo) then
       crystallite_subF0 (1:3,1:3,co,ip,el) = crystallite_subF(1:3,1:3,co,ip,el)
       subLp0 = crystallite_Lp  (1:3,1:3,co,ip,el)
       subLi0 = constitutive_mech_Li(ph)%data(1:3,1:3,me)
       crystallite_subFp0(1:3,1:3,co,ip,el) = constitutive_mech_Fp(ph)%data(1:3,1:3,me)
       crystallite_subFi0(1:3,1:3,co,ip,el) = constitutive_mech_Fi(ph)%data(1:3,1:3,me)
       plasticState(    material_phaseAt(co,el))%subState0(:,material_phaseMemberAt(co,ip,el)) &
         = plasticState(material_phaseAt(co,el))%state(    :,material_phaseMemberAt(co,ip,el))
       do s = 1, phase_Nsources(material_phaseAt(co,el))
         sourceState(    material_phaseAt(co,el))%p(s)%subState0(:,material_phaseMemberAt(co,ip,el)) &
           = sourceState(material_phaseAt(co,el))%p(s)%state(    :,material_phaseMemberAt(co,ip,el))
       enddo
     endif

!--------------------------------------------------------------------------------------------------
!  cut back (reduced time and restore)
   else
     subStep       = num%subStepSizeCryst * subStep
     constitutive_mech_Fp(ph)%data(1:3,1:3,me) =            crystallite_subFp0(1:3,1:3,co,ip,el)
     constitutive_mech_Fi(ph)%data(1:3,1:3,me) =            crystallite_subFi0(1:3,1:3,co,ip,el)
     crystallite_S    (1:3,1:3,co,ip,el) =            crystallite_S0    (1:3,1:3,co,ip,el)
     if (subStep < 1.0_pReal) then                                        ! actual (not initial) cutback
       crystallite_Lp (1:3,1:3,co,ip,el) =            subLp0
       constitutive_mech_Li(ph)%data(1:3,1:3,me) =            subLi0
     endif
     plasticState    (material_phaseAt(co,el))%state(    :,material_phaseMemberAt(co,ip,el)) &
       = plasticState(material_phaseAt(co,el))%subState0(:,material_phaseMemberAt(co,ip,el))
     do s = 1, phase_Nsources(material_phaseAt(co,el))
       sourceState(    material_phaseAt(co,el))%p(s)%state(    :,material_phaseMemberAt(co,ip,el)) &
         = sourceState(material_phaseAt(co,el))%p(s)%subState0(:,material_phaseMemberAt(co,ip,el))
     enddo

                                                                                             ! cant restore dotState here, since not yet calculated in first cutback after initialization
     todo = subStep > num%subStepMinCryst                          ! still on track or already done (beyond repair)
   endif

!--------------------------------------------------------------------------------------------------
!  prepare for integration
   if (todo) then
     crystallite_subF(1:3,1:3,co,ip,el) = crystallite_subF0(1:3,1:3,co,ip,el) &
                                        + subStep *( crystallite_partitionedF (1:3,1:3,co,ip,el) &
                                                    -crystallite_partitionedF0(1:3,1:3,co,ip,el))
     crystallite_Fe(1:3,1:3,co,ip,el) = matmul(crystallite_subF(1:3,1:3,co,ip,el), &
                                            math_inv33(matmul(constitutive_mech_Fi(ph)%data(1:3,1:3,me), &
                                                              constitutive_mech_Fp(ph)%data(1:3,1:3,me))))
     crystallite_subdt(co,ip,el) = subStep * crystallite_dt(co,ip,el)
     crystallite_converged(co,ip,el) = .false.
     call integrateState(co,ip,el)
     call integrateSourceState(co,ip,el)
   endif

   if (.not. crystallite_converged(co,ip,el) .and. subStep > num%subStepMinCryst) &            ! do not try non-converged but fully cutbacked any further
      todo = .true.
  enddo cutbackLooping

! return whether converged or not
  crystallite_stress = crystallite_converged(co,ip,el)

end function crystallite_stress


!--------------------------------------------------------------------------------------------------
!> @brief Backup data for homog cutback.
!--------------------------------------------------------------------------------------------------
subroutine constitutive_initializeRestorationPoints(ip,el)

  integer, intent(in) :: &
    ip, &                                                                                            !< integration point number
    el                                                                                               !< element number
  integer :: &
    co, &                                                                                            !< constituent number
    s,ph, me

  do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
    ph = material_phaseAt(co,el)
    me = material_phaseMemberAt(co,ip,el)
    crystallite_partitionedLp0(1:3,1:3,co,ip,el) = crystallite_Lp0(1:3,1:3,co,ip,el)
    crystallite_partitionedF0(1:3,1:3,co,ip,el)  = crystallite_F0(1:3,1:3,co,ip,el)
    crystallite_partitionedS0(1:3,1:3,co,ip,el)  = crystallite_S0(1:3,1:3,co,ip,el)

    call mech_initializeRestorationPoints(ph,me)

    do s = 1, phase_Nsources(material_phaseAt(co,el))
      sourceState(material_phaseAt(co,el))%p(s)%partitionedState0(:,material_phasememberAt(co,ip,el)) = &
      sourceState(material_phaseAt(co,el))%p(s)%state0(           :,material_phasememberAt(co,ip,el))
    enddo
  enddo

end subroutine constitutive_initializeRestorationPoints


!--------------------------------------------------------------------------------------------------
!> @brief Wind homog inc forward.
!--------------------------------------------------------------------------------------------------
subroutine constitutive_windForward(ip,el)

  integer, intent(in) :: &
    ip, &                                                                                            !< integration point number
    el                                                                                               !< element number
  integer :: &
    co, &                                                                                            !< constituent number
    s, ph, me
  do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
    ph = material_phaseAt(co,el)
    me = material_phaseMemberAt(co,ip,el)
    crystallite_partitionedF0 (1:3,1:3,co,ip,el) = crystallite_partitionedF(1:3,1:3,co,ip,el)
    crystallite_partitionedLp0(1:3,1:3,co,ip,el) = crystallite_Lp        (1:3,1:3,co,ip,el)
    crystallite_partitionedS0 (1:3,1:3,co,ip,el) = crystallite_S         (1:3,1:3,co,ip,el)

    call constitutive_mech_windForward(ph,me)
    do s = 1, phase_Nsources(material_phaseAt(co,el))
      sourceState(ph)%p(s)%partitionedState0(:,me) = sourceState(ph)%p(s)%state(:,me)
    enddo
  enddo

end subroutine constitutive_windForward


!--------------------------------------------------------------------------------------------------
!> @brief Restore data after homog cutback.
!--------------------------------------------------------------------------------------------------
subroutine crystallite_restore(ip,el,includeL)

  integer, intent(in) :: &
    ip, &                                                                                            !< integration point number
    el                                                                                               !< element number
  logical, intent(in) :: &
    includeL                                                                                        !< protect agains fake cutback
  integer :: &
    co, p, m                                                                                            !< constituent number

  do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
  p = material_phaseAt(co,el)
  m = material_phaseMemberAt(co,ip,el)
    if (includeL) then
      crystallite_Lp(1:3,1:3,co,ip,el) = crystallite_partitionedLp0(1:3,1:3,co,ip,el)
      constitutive_mech_Li(p)%data(1:3,1:3,m) = constitutive_mech_partionedLi0(p)%data(1:3,1:3,m)
    endif                                                                                           ! maybe protecting everything from overwriting makes more sense

    constitutive_mech_Fp(p)%data(1:3,1:3,m)   = constitutive_mech_partionedFp0(p)%data(1:3,1:3,m)
    constitutive_mech_Fi(p)%data(1:3,1:3,m)   = constitutive_mech_partionedFi0(p)%data(1:3,1:3,m)
    crystallite_S (1:3,1:3,co,ip,el)   = crystallite_partitionedS0 (1:3,1:3,co,ip,el)

    plasticState    (material_phaseAt(co,el))%state(          :,material_phasememberAt(co,ip,el)) = &
    plasticState    (material_phaseAt(co,el))%partitionedState0(:,material_phasememberAt(co,ip,el))
  enddo

end subroutine crystallite_restore


!--------------------------------------------------------------------------------------------------
!> @brief Calculate tangent (dPdF).
!--------------------------------------------------------------------------------------------------
function crystallite_stressTangent(co,ip,el) result(dPdF)

  real(pReal), dimension(3,3,3,3) :: dPdF
  integer, intent(in) :: &
    co, &                                                                                            !< counter in constituent loop
    ip, &                                                                                            !< counter in integration point loop
    el                                                                                               !< counter in element loop
  integer :: &
    o, &
    p, pp, m

  real(pReal), dimension(3,3)     ::   devNull, &
                                       invSubFp0,invSubFi0,invFp,invFi, &
                                       temp_33_1, temp_33_2, temp_33_3, temp_33_4
  real(pReal), dimension(3,3,3,3) ::   dSdFe, &
                                       dSdF, &
                                       dSdFi, &
                                       dLidS, &                                                     ! tangent in lattice configuration
                                       dLidFi, &
                                       dLpdS, &
                                       dLpdFi, &
                                       dFidS, &
                                       dFpinvdF, &
                                       rhs_3333, &
                                       lhs_3333, &
                                       temp_3333
  real(pReal), dimension(9,9)::        temp_99
  logical :: error

  pp = material_phaseAt(co,el)
  m = material_phaseMemberAt(co,ip,el)

  call constitutive_hooke_SandItsTangents(devNull,dSdFe,dSdFi, &
                                   crystallite_Fe(1:3,1:3,co,ip,el), &
                                   constitutive_mech_Fi(pp)%data(1:3,1:3,m),co,ip,el)
  call constitutive_LiAndItsTangents(devNull,dLidS,dLidFi, &
                                     crystallite_S (1:3,1:3,co,ip,el), &
                                     constitutive_mech_Fi(pp)%data(1:3,1:3,m), &
                                     co,ip,el)

  invFp = math_inv33(constitutive_mech_Fp(pp)%data(1:3,1:3,m))
  invFi = math_inv33(constitutive_mech_Fi(pp)%data(1:3,1:3,m))
  invSubFp0 = math_inv33(crystallite_subFp0(1:3,1:3,co,ip,el))
  invSubFi0 = math_inv33(crystallite_subFi0(1:3,1:3,co,ip,el))

  if (sum(abs(dLidS)) < tol_math_check) then
    dFidS = 0.0_pReal
  else
    lhs_3333 = 0.0_pReal; rhs_3333 = 0.0_pReal
    do o=1,3; do p=1,3
      lhs_3333(1:3,1:3,o,p) = lhs_3333(1:3,1:3,o,p) &
                            + crystallite_subdt(co,ip,el)*matmul(invSubFi0,dLidFi(1:3,1:3,o,p))
      lhs_3333(1:3,o,1:3,p) = lhs_3333(1:3,o,1:3,p) &
                            + invFi*invFi(p,o)
      rhs_3333(1:3,1:3,o,p) = rhs_3333(1:3,1:3,o,p) &
                            - crystallite_subdt(co,ip,el)*matmul(invSubFi0,dLidS(1:3,1:3,o,p))
    enddo; enddo
    call math_invert(temp_99,error,math_3333to99(lhs_3333))
    if (error) then
      call IO_warning(warning_ID=600,el=el,ip=ip,g=co, &
                      ext_msg='inversion error in analytic tangent calculation')
      dFidS = 0.0_pReal
    else
      dFidS = math_mul3333xx3333(math_99to3333(temp_99),rhs_3333)
    endif
    dLidS = math_mul3333xx3333(dLidFi,dFidS) + dLidS
  endif

  call constitutive_plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, &
                                     crystallite_S (1:3,1:3,co,ip,el), &
                                     constitutive_mech_Fi(pp)%data(1:3,1:3,m),co,ip,el)
  dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS

!--------------------------------------------------------------------------------------------------
! calculate dSdF
  temp_33_1 = transpose(matmul(invFp,invFi))
  temp_33_2 = matmul(crystallite_subF(1:3,1:3,co,ip,el),invSubFp0)
  temp_33_3 = matmul(matmul(crystallite_subF(1:3,1:3,co,ip,el),invFp), invSubFi0)

  do o=1,3; do p=1,3
    rhs_3333(p,o,1:3,1:3)  = matmul(dSdFe(p,o,1:3,1:3),temp_33_1)
    temp_3333(1:3,1:3,p,o) = matmul(matmul(temp_33_2,dLpdS(1:3,1:3,p,o)), invFi) &
                           + matmul(temp_33_3,dLidS(1:3,1:3,p,o))
  enddo; enddo
  lhs_3333 = crystallite_subdt(co,ip,el)*math_mul3333xx3333(dSdFe,temp_3333) &
           + math_mul3333xx3333(dSdFi,dFidS)

  call math_invert(temp_99,error,math_eye(9)+math_3333to99(lhs_3333))
  if (error) then
    call IO_warning(warning_ID=600,el=el,ip=ip,g=co, &
                    ext_msg='inversion error in analytic tangent calculation')
    dSdF = rhs_3333
  else
    dSdF = math_mul3333xx3333(math_99to3333(temp_99),rhs_3333)
  endif

!--------------------------------------------------------------------------------------------------
! calculate dFpinvdF
  temp_3333 = math_mul3333xx3333(dLpdS,dSdF)
  do o=1,3; do p=1,3
    dFpinvdF(1:3,1:3,p,o) = -crystallite_subdt(co,ip,el) &
                          * matmul(invSubFp0, matmul(temp_3333(1:3,1:3,p,o),invFi))
  enddo; enddo

!--------------------------------------------------------------------------------------------------
! assemble dPdF
  temp_33_1 = matmul(crystallite_S(1:3,1:3,co,ip,el),transpose(invFp))
  temp_33_2 = matmul(invFp,temp_33_1)
  temp_33_3 = matmul(crystallite_subF(1:3,1:3,co,ip,el),invFp)
  temp_33_4 = matmul(temp_33_3,crystallite_S(1:3,1:3,co,ip,el))

  dPdF = 0.0_pReal
  do p=1,3
    dPdF(p,1:3,p,1:3) = transpose(temp_33_2)
  enddo
  do o=1,3; do p=1,3
    dPdF(1:3,1:3,p,o) = dPdF(1:3,1:3,p,o) &
                      + matmul(matmul(crystallite_subF(1:3,1:3,co,ip,el), &
                               dFpinvdF(1:3,1:3,p,o)),temp_33_1) &
                      + matmul(matmul(temp_33_3,dSdF(1:3,1:3,p,o)), &
                               transpose(invFp)) &
                      + matmul(temp_33_4,transpose(dFpinvdF(1:3,1:3,p,o)))
  enddo; enddo

end function crystallite_stressTangent


!--------------------------------------------------------------------------------------------------
!> @brief calculates orientations
!--------------------------------------------------------------------------------------------------
subroutine crystallite_orientations

  integer &
    co, &                                                                                            !< counter in integration point component loop
    ip, &                                                                                            !< counter in integration point loop
    el                                                                                               !< counter in element loop


  !$OMP PARALLEL DO
  do el = FEsolving_execElem(1),FEsolving_execElem(2)
    do ip = FEsolving_execIP(1),FEsolving_execIP(2)
      do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
        call crystallite_orientation(co,ip,el)%fromMatrix(transpose(math_rotationalPart(crystallite_Fe(1:3,1:3,co,ip,el))))
  enddo; enddo; enddo
  !$OMP END PARALLEL DO

  nonlocalPresent: if (any(plasticState%nonlocal)) then
    !$OMP PARALLEL DO
    do el = FEsolving_execElem(1),FEsolving_execElem(2)
      if (plasticState(material_phaseAt(1,el))%nonlocal) then
        do ip = FEsolving_execIP(1),FEsolving_execIP(2)
          call plastic_nonlocal_updateCompatibility(crystallite_orientation, &
                                                    phase_plasticityInstance(material_phaseAt(1,el)),ip,el)
        enddo
      endif
    enddo
    !$OMP END PARALLEL DO
  endif nonlocalPresent

end subroutine crystallite_orientations


!--------------------------------------------------------------------------------------------------
!> @brief Map 2nd order tensor to reference config
!--------------------------------------------------------------------------------------------------
function crystallite_push33ToRef(co,ip,el, tensor33)

  real(pReal), dimension(3,3), intent(in) :: tensor33
  real(pReal), dimension(3,3)             :: T
  integer, intent(in):: &
    el, &
    ip, &
    co

  real(pReal), dimension(3,3) :: crystallite_push33ToRef


  T = matmul(material_orientation0(co,ip,el)%asMatrix(), &                                         ! ToDo: initial orientation correct?
             transpose(math_inv33(crystallite_subF(1:3,1:3,co,ip,el))))
  crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T))

end function crystallite_push33ToRef


!--------------------------------------------------------------------------------------------------
!> @brief integrate stress, state with adaptive 1st order explicit Euler method
!> using Fixed Point Iteration to adapt the stepsize
!--------------------------------------------------------------------------------------------------
subroutine integrateSourceState(co,ip,el)

  integer, intent(in) :: &
    el, &                                                                                            !< element index in element loop
    ip, &                                                                                            !< integration point index in ip loop
    co                                                                                               !< grain index in grain loop

  integer :: &
    NiterationState, &                                                                              !< number of iterations in state loop
    ph, &
    me, &
    so
  integer, dimension(maxval(phase_Nsources)) :: &
    size_so
  real(pReal) :: &
    zeta
  real(pReal), dimension(constitutive_source_maxSizeDotState) :: &
    r                                                                                               ! state residuum
  real(pReal), dimension(constitutive_source_maxSizeDotState,2,maxval(phase_Nsources)) :: source_dotState
  logical :: &
    broken


  ph = material_phaseAt(co,el)
  me = material_phaseMemberAt(co,ip,el)

  broken = constitutive_thermal_collectDotState(ph,me)
  broken = broken .or. constitutive_damage_collectDotState(crystallite_S(1:3,1:3,co,ip,el), co,ip,el,ph,me)
  if(broken) return

  do so = 1, phase_Nsources(ph)
    size_so(so) = sourceState(ph)%p(so)%sizeDotState
    sourceState(ph)%p(so)%state(1:size_so(so),me) = sourceState(ph)%p(so)%subState0(1:size_so(so),me) &
                                              + sourceState(ph)%p(so)%dotState (1:size_so(so),me) &
                                              * crystallite_subdt(co,ip,el)
    source_dotState(1:size_so(so),2,so) = 0.0_pReal
  enddo

  iteration: do NiterationState = 1, num%nState

    do so = 1, phase_Nsources(ph)
      if(nIterationState > 1) source_dotState(1:size_so(so),2,so) = source_dotState(1:size_so(so),1,so)
      source_dotState(1:size_so(so),1,so) = sourceState(ph)%p(so)%dotState(:,me)
    enddo

    broken = constitutive_thermal_collectDotState(ph,me)
    broken = broken .or. constitutive_damage_collectDotState(crystallite_S(1:3,1:3,co,ip,el), co,ip,el,ph,me)
    if(broken) exit iteration

    do so = 1, phase_Nsources(ph)
      zeta = damper(sourceState(ph)%p(so)%dotState(:,me), &
                    source_dotState(1:size_so(so),1,so),&
                    source_dotState(1:size_so(so),2,so))
      sourceState(ph)%p(so)%dotState(:,me) = sourceState(ph)%p(so)%dotState(:,me) * zeta &
                                        + source_dotState(1:size_so(so),1,so)* (1.0_pReal - zeta)
      r(1:size_so(so)) = sourceState(ph)%p(so)%state    (1:size_so(so),me)  &
                      - sourceState(ph)%p(so)%subState0(1:size_so(so),me)  &
                      - sourceState(ph)%p(so)%dotState (1:size_so(so),me) * crystallite_subdt(co,ip,el)
      sourceState(ph)%p(so)%state(1:size_so(so),me) = sourceState(ph)%p(so)%state(1:size_so(so),me) &
                                                - r(1:size_so(so))
      crystallite_converged(co,ip,el) = &
      crystallite_converged(co,ip,el) .and. converged(r(1:size_so(so)), &
                                                   sourceState(ph)%p(so)%state(1:size_so(so),me), &
                                                   sourceState(ph)%p(so)%atol(1:size_so(so)))
    enddo

    if(crystallite_converged(co,ip,el)) then
      broken = constitutive_damage_deltaState(crystallite_Fe(1:3,1:3,co,ip,el),co,ip,el,ph,me)
      exit iteration
    endif

  enddo iteration


  contains

  !--------------------------------------------------------------------------------------------------
  !> @brief calculate the damping for correction of state and dot state
  !--------------------------------------------------------------------------------------------------
  real(pReal) pure function damper(current,previous,previous2)

  real(pReal), dimension(:), intent(in) ::&
    current, previous, previous2

  real(pReal) :: dot_prod12, dot_prod22

  dot_prod12 = dot_product(current  - previous,  previous - previous2)
  dot_prod22 = dot_product(previous - previous2, previous - previous2)
  if ((dot_product(current,previous) < 0.0_pReal .or. dot_prod12 < 0.0_pReal) .and. dot_prod22 > 0.0_pReal) then
    damper = 0.75_pReal + 0.25_pReal * tanh(2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22)
  else
    damper = 1.0_pReal
  endif

  end function damper

end subroutine integrateSourceState


!--------------------------------------------------------------------------------------------------
!> @brief determines whether a point is converged
!--------------------------------------------------------------------------------------------------
logical pure function converged(residuum,state,atol)

  real(pReal), intent(in), dimension(:) ::&
    residuum, state, atol
  real(pReal) :: &
    rTol

  rTol = num%rTol_crystalliteState

  converged = all(abs(residuum) <= max(atol, rtol*abs(state)))

end function converged


!--------------------------------------------------------------------------------------------------
!> @brief Write current  restart information (Field and constitutive data) to file.
! ToDo: Merge data into one file for MPI, move state to constitutive and homogenization, respectively
!--------------------------------------------------------------------------------------------------
subroutine crystallite_restartWrite

  integer :: i
  integer(HID_T) :: fileHandle, groupHandle
  character(len=pStringLen) :: fileName, datasetName

  print*, ' writing field and constitutive data required for restart to file';flush(IO_STDOUT)

  write(fileName,'(a,i0,a)') trim(getSolverJobName())//'_',worldrank,'.hdf5'
  fileHandle = HDF5_openFile(fileName,'a')

  call HDF5_write(fileHandle,crystallite_partitionedF,'F')
  call HDF5_write(fileHandle,crystallite_Lp,        'L_p')
  call HDF5_write(fileHandle,crystallite_S,           'S')

  groupHandle = HDF5_addGroup(fileHandle,'phase')
  do i = 1,size(material_name_phase)
    write(datasetName,'(i0,a)') i,'_omega'
    call HDF5_write(groupHandle,plasticState(i)%state,datasetName)
    write(datasetName,'(i0,a)') i,'_F_i'
    call HDF5_write(groupHandle,constitutive_mech_Fi(i)%data,datasetName)
    write(datasetName,'(i0,a)') i,'_L_i'
    call HDF5_write(groupHandle,constitutive_mech_Li(i)%data,datasetName)
    write(datasetName,'(i0,a)') i,'_F_p'
    call HDF5_write(groupHandle,constitutive_mech_Fp(i)%data,datasetName)
  enddo
  call HDF5_closeGroup(groupHandle)

  groupHandle = HDF5_addGroup(fileHandle,'homogenization')
  do i = 1, size(material_name_homogenization)
    write(datasetName,'(i0,a)') i,'_omega'
    call HDF5_write(groupHandle,homogState(i)%state,datasetName)
  enddo
  call HDF5_closeGroup(groupHandle)

  call HDF5_closeFile(fileHandle)

end subroutine crystallite_restartWrite


!--------------------------------------------------------------------------------------------------
!> @brief Read data for restart
! ToDo: Merge data into one file for MPI, move state to constitutive and homogenization, respectively
!--------------------------------------------------------------------------------------------------
subroutine crystallite_restartRead

  integer :: i
  integer(HID_T) :: fileHandle, groupHandle
  character(len=pStringLen) :: fileName, datasetName

  print'(/,a,i0,a)', ' reading restart information of increment from file'

  write(fileName,'(a,i0,a)') trim(getSolverJobName())//'_',worldrank,'.hdf5'
  fileHandle = HDF5_openFile(fileName)

  call HDF5_read(fileHandle,crystallite_F0, 'F')
  call HDF5_read(fileHandle,crystallite_Lp0,'L_p')
  call HDF5_read(fileHandle,crystallite_S0, 'S')

  groupHandle = HDF5_openGroup(fileHandle,'phase')
  do i = 1,size(material_name_phase)
    write(datasetName,'(i0,a)') i,'_omega'
    call HDF5_read(groupHandle,plasticState(i)%state0,datasetName)
    write(datasetName,'(i0,a)') i,'_F_i'
    call HDF5_read(groupHandle,constitutive_mech_Fi0(i)%data,datasetName)
    write(datasetName,'(i0,a)') i,'_L_i'
    call HDF5_read(groupHandle,constitutive_mech_Li0(i)%data,datasetName)
    write(datasetName,'(i0,a)') i,'_F_p'
    call HDF5_read(groupHandle,constitutive_mech_Fp0(i)%data,datasetName)
  enddo
  call HDF5_closeGroup(groupHandle)

  groupHandle = HDF5_openGroup(fileHandle,'homogenization')
  do i = 1,size(material_name_homogenization)
    write(datasetName,'(i0,a)') i,'_omega'
    call HDF5_read(groupHandle,homogState(i)%state0,datasetName)
  enddo
  call HDF5_closeGroup(groupHandle)

  call HDF5_closeFile(fileHandle)

end subroutine crystallite_restartRead


end module constitutive