should become mech only

This commit is contained in:
Martin Diehl 2020-12-24 10:22:41 +01:00
parent 935b531d27
commit acc998d242
2 changed files with 127 additions and 119 deletions

View File

@ -117,7 +117,7 @@ module constitutive
type(tDebugOptions) :: debugCrystallite
procedure(integrateStateFPI), pointer :: integrateState
integer(kind(PLASTICITY_undefined_ID)), dimension(:), allocatable, public :: &
phase_plasticity !< plasticity of each phase
@ -187,6 +187,12 @@ module constitutive
! == cleaned:end ===================================================================================
module function crystallite_stress(dt,co,ip,el)
real(pReal), intent(in) :: dt
integer, intent(in) :: co, ip, el
logical :: crystallite_stress
end function crystallite_stress
module function constitutive_homogenizedC(co,ip,el) result(C)
integer, intent(in) :: co, ip, el
real(pReal), dimension(6,6) :: C
@ -362,10 +368,6 @@ module constitutive
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
@ -392,6 +394,7 @@ module constitutive
crystallite_orientations, &
crystallite_push33ToRef, &
crystallite_restartWrite, &
integrateSourceState, &
crystallite_restartRead, &
constitutive_initializeRestorationPoints, &
constitutive_windForward, &
@ -983,120 +986,7 @@ subroutine crystallite_init
end subroutine crystallite_init
!--------------------------------------------------------------------------------------------------
!> @brief calculate stress (P)
!--------------------------------------------------------------------------------------------------
function crystallite_stress(dt,co,ip,el)
real(pReal), intent(in) :: dt
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 * dt
crystallite_converged(co,ip,el) = .false.
call integrateState(co,ip,el)
call integrateSourceState(co,ip,el)
endif
enddo cutbackLooping
! return whether converged or not
crystallite_stress = crystallite_converged(co,ip,el)
end function crystallite_stress
!--------------------------------------------------------------------------------------------------

View File

@ -290,6 +290,8 @@ submodule(constitutive) constitutive_mech
end type tOutput
type(tOutput), allocatable, dimension(:) :: output_constituent
procedure(integrateStateFPI), pointer :: integrateState
contains
@ -959,7 +961,7 @@ end function integrateStress
!> @brief integrate stress, state with adaptive 1st order explicit Euler method
!> using Fixed Point Iteration to adapt the stepsize
!--------------------------------------------------------------------------------------------------
module subroutine integrateStateFPI(co,ip,el)
subroutine integrateStateFPI(co,ip,el)
integer, intent(in) :: &
el, & !< element index in element loop
@ -1475,5 +1477,121 @@ module function constitutive_homogenizedC(co,ip,el) result(C)
end function constitutive_homogenizedC
!--------------------------------------------------------------------------------------------------
!> @brief calculate stress (P)
!--------------------------------------------------------------------------------------------------
module function crystallite_stress(dt,co,ip,el)
real(pReal), intent(in) :: dt
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 * dt
crystallite_converged(co,ip,el) = .false.
call integrateState(co,ip,el)
call integrateSourceState(co,ip,el)
endif
enddo cutbackLooping
! return whether converged or not
crystallite_stress = crystallite_converged(co,ip,el)
end function crystallite_stress
end submodule constitutive_mech