diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 582222dae..e7f5628ed 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -709,12 +709,14 @@ end subroutine constitutive_hooke_SandItsTangents !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the rate of change of microstructure !-------------------------------------------------------------------------------------------------- -function constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el) result(broken) +function constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el,phase,of) result(broken) integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point - el !< element + el, & !< element + phase, & + of real(pReal), intent(in) :: & subdt !< timestep real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & @@ -727,17 +729,14 @@ function constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el real(pReal), dimension(3,3) :: & Mp integer :: & - phase, & ho, & !< homogenization tme, & !< thermal member position i, & !< counter in source loop - instance, of + instance logical :: broken ho = material_homogenizationAt(el) tme = thermalMapping(ho)%p(ip,el) - of = material_phasememberAt(ipc,ip,el) - phase = material_phaseAt(ipc,el) instance = phase_plasticityInstance(phase) Mp = matmul(matmul(transpose(Fi),Fi),S) @@ -794,12 +793,14 @@ end function constitutive_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_deltaState(S, Fe, Fi, ipc, ip, el) result(broken) +function constitutive_deltaState(S, Fe, Fi, ipc, ip, el, phase, of) result(broken) integer, intent(in) :: & ipc, & !< component-ID of integration point ip, & !< integration point - el !< element + el, & !< element + phase, & + of real(pReal), intent(in), dimension(3,3) :: & S, & !< 2nd Piola Kirchhoff stress Fe, & !< elastic deformation gradient @@ -808,27 +809,28 @@ function constitutive_deltaState(S, Fe, Fi, ipc, ip, el) result(broken) Mp integer :: & i, & - instance, of, & - phase + instance logical :: & broken - Mp = matmul(matmul(transpose(Fi),Fi),S) - of = material_phasememberAt(ipc,ip,el) - phase = material_phaseAt(ipc,el) - instance = phase_plasticityInstance(material_phaseAt(ipc,el)) + instance = phase_plasticityInstance(phase) plasticityType: select case (phase_plasticity(phase)) case (PLASTICITY_KINEHARDENING_ID) plasticityType call plastic_kinehardening_deltaState(Mp,instance,of) + broken = any(IEEE_is_NaN(plasticState(phase)%deltaState(:,of))) case (PLASTICITY_NONLOCAL_ID) plasticityType call plastic_nonlocal_deltaState(Mp,instance,of,ip,el) + broken = any(IEEE_is_NaN(plasticState(phase)%deltaState(:,of))) + + case default + broken = .false. end select plasticityType - broken = any(IEEE_is_NaN(plasticState(phase)%deltaState(:,of))) + sourceLoop: do i = 1, phase_Nsources(phase) @@ -837,11 +839,10 @@ function constitutive_deltaState(S, Fe, Fi, ipc, ip, el) result(broken) case (SOURCE_damage_isoBrittle_ID) sourceType call source_damage_isoBrittle_deltaState (constitutive_homogenizedC(ipc,ip,el), Fe, & ipc, ip, el) + broken = broken .or. any(IEEE_is_NaN(sourceState(phase)%p(i)%deltaState(:,of))) end select sourceType - broken = broken .or. any(IEEE_is_NaN(sourceState(phase)%p(i)%deltaState(:,of))) - enddo SourceLoop end function constitutive_deltaState diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 611a7551d..054bc2105 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -1038,7 +1038,7 @@ subroutine integrateStateFPI(todo) crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e) + crystallite_subdt(g,i,e), g,i,e,p,c) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle @@ -1073,7 +1073,7 @@ subroutine integrateStateFPI(todo) crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e) + crystallite_subdt(g,i,e), g,i,e,p,c) if(broken) exit iteration sizeDotState = plasticState(p)%sizeDotState @@ -1107,7 +1107,7 @@ subroutine integrateStateFPI(todo) enddo if(crystallite_converged(g,i,e)) then - broken = stateJump(g,i,e) + broken = stateJump(g,i,e,p,c) exit iteration endif @@ -1177,7 +1177,7 @@ subroutine integrateStateEuler(todo) crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e) + crystallite_subdt(g,i,e), g,i,e,p,c) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle @@ -1193,7 +1193,7 @@ subroutine integrateStateEuler(todo) * crystallite_subdt(g,i,e) enddo - broken = stateJump(g,i,e) + broken = stateJump(g,i,e,p,c) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle @@ -1246,7 +1246,7 @@ subroutine integrateStateAdaptiveEuler(todo) crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e) + crystallite_subdt(g,i,e), g,i,e,p,c) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle @@ -1265,7 +1265,7 @@ subroutine integrateStateAdaptiveEuler(todo) + sourceState(p)%p(s)%dotstate(1:sizeDotState,c) * crystallite_subdt(g,i,e) enddo - broken = stateJump(g,i,e) + broken = stateJump(g,i,e,p,c) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle @@ -1277,7 +1277,7 @@ subroutine integrateStateAdaptiveEuler(todo) crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e) + crystallite_subdt(g,i,e), g,i,e,p,c) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle @@ -1357,7 +1357,7 @@ subroutine integrateStateRK4(todo) crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e) + crystallite_subdt(g,i,e), g,i,e,p,c) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle @@ -1401,7 +1401,7 @@ subroutine integrateStateRK4(todo) crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & - crystallite_subdt(g,i,e)*CC(stage), g,i,e) + crystallite_subdt(g,i,e)*CC(stage), g,i,e,p,c) if(broken) exit enddo @@ -1428,7 +1428,7 @@ subroutine integrateStateRK4(todo) * crystallite_subdt(g,i,e) enddo - broken = stateJump(g,i,e) + broken = stateJump(g,i,e,p,c) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle @@ -1502,7 +1502,7 @@ subroutine integrateStateRKCK45(todo) crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & - crystallite_subdt(g,i,e), g,i,e) + crystallite_subdt(g,i,e), g,i,e,p,c) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle @@ -1546,7 +1546,7 @@ subroutine integrateStateRKCK45(todo) crystallite_partionedF0, & crystallite_Fi(1:3,1:3,g,i,e), & crystallite_partionedFp0, & - crystallite_subdt(g,i,e)*CC(stage), g,i,e) + crystallite_subdt(g,i,e)*CC(stage), g,i,e,p,c) if(broken) exit enddo @@ -1582,7 +1582,7 @@ subroutine integrateStateRKCK45(todo) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle - broken = stateJump(g,i,e) + broken = stateJump(g,i,e,p,c) if(broken .and. plasticState(p)%nonlocal) nonlocalBroken = .true. if(broken) cycle @@ -1631,7 +1631,7 @@ end function converged !> @brief calculates a jump in the state according to the current state and the current stress !> returns true, if state jump was successfull or not needed. false indicates NaN in delta state !-------------------------------------------------------------------------------------------------- -function stateJump(ipc,ip,el) result(broken) +function stateJump(ipc,ip,el,p,c) result(broken) integer, intent(in):: & el, & ! element index @@ -1646,13 +1646,10 @@ function stateJump(ipc,ip,el) result(broken) mySize logical :: broken - c = material_phaseMemberAt(ipc,ip,el) - p = material_phaseAt(ipc,el) - broken = constitutive_deltaState(crystallite_S(1:3,1:3,ipc,ip,el), & crystallite_Fe(1:3,1:3,ipc,ip,el), & crystallite_Fi(1:3,1:3,ipc,ip,el), & - ipc,ip,el) + ipc,ip,el,p,c) if(broken) return myOffset = plasticState(p)%offsetDeltaState