788 lines
41 KiB
Fortran
788 lines
41 KiB
Fortran
|
|
!***************************************
|
|
!* Module: CRYSTALLITE *
|
|
!***************************************
|
|
!* contains: *
|
|
!* - _init *
|
|
!* - materialpoint_stressAndItsTangent *
|
|
!* - _partitionDeformation *
|
|
!* - _updateState *
|
|
!* - _averageStressAndItsTangent *
|
|
!* - _postResults *
|
|
!***************************************
|
|
|
|
MODULE crystallite
|
|
|
|
use prec, only: pReal,pInt
|
|
implicit none
|
|
!
|
|
! ****************************************************************
|
|
! *** General variables for the crystallite calculation ***
|
|
! ****************************************************************
|
|
integer(pInt), parameter :: crystallite_Nresults = 5_pInt ! phaseID, volfrac within this phase, Euler angles
|
|
|
|
real(pReal), dimension (:,:,:,:,:), allocatable :: crystallite_Fe, & ! current "elastic" def grad (end of converged time step)
|
|
crystallite_Fp, & ! current plastic def grad (end of converged time step)
|
|
crystallite_Lp, & ! current plastic velocitiy grad (end of converged time step)
|
|
crystallite_F0, & ! def grad at start of FE inc
|
|
crystallite_Fp0, & ! plastic def grad at start of FE inc
|
|
crystallite_Lp0, & ! plastic velocitiy grad at start of FE inc
|
|
crystallite_partionedF, & ! def grad to be reached at end of homog inc
|
|
crystallite_partionedF0, & ! def grad at start of homog inc
|
|
crystallite_partionedFp0,& ! plastic def grad at start of homog inc
|
|
crystallite_partionedLp0,& ! plastic velocity grad at start of homog inc
|
|
crystallite_subF, & ! def grad to be reached at end of crystallite inc
|
|
crystallite_subF0, & ! def grad at start of crystallite inc
|
|
crystallite_subFp0,& ! plastic def grad at start of crystallite inc
|
|
crystallite_subLp0,& ! plastic velocity grad at start of crystallite inc
|
|
crystallite_P ! 1st Piola-Kirchhoff stress per grain
|
|
real(pReal), dimension (:,:,:,:), allocatable :: crystallite_Tstar_v ! 2nd Piola-Kirchhoff stress (vector) per grain
|
|
real(pReal), dimension (:,:,:,:,:,:,:),allocatable :: crystallite_dPdF, & ! individual dPdF per grain
|
|
crystallite_fallbackdPdF ! dPdF fallback for non-converged grains (elastic prediction)
|
|
real(pReal), dimension (:,:,:), allocatable :: crystallite_dt, & ! requested time increment of each grain
|
|
crystallite_subdt, & ! substepped time increment of each grain
|
|
crystallite_subFrac, & ! already calculated fraction of increment
|
|
crystallite_subStep, & ! size of next integration step
|
|
crystallite_Temperature ! Temp of each grain
|
|
|
|
logical, dimension (:,:,:), allocatable :: crystallite_localConstitution, & ! indicates this grain to have purely local constitutive law
|
|
crystallite_requested, & ! flag to request crystallite calculation
|
|
crystallite_onTrack, & ! flag to indicate ongoing calculation
|
|
crystallite_converged ! convergence flag
|
|
|
|
|
|
CONTAINS
|
|
|
|
!********************************************************************
|
|
! allocate and initialize per grain variables
|
|
!********************************************************************
|
|
subroutine crystallite_init()
|
|
|
|
use prec, only: pInt,pReal
|
|
use debug, only: debug_info,debug_reset
|
|
use math, only: math_I3,math_EulerToR
|
|
use FEsolving, only: FEsolving_execElem,FEsolving_execIP
|
|
use mesh, only: mesh_element,mesh_NcpElems,mesh_maxNips
|
|
use material, only: homogenization_Ngrains,homogenization_maxNgrains,&
|
|
material_EulerAngles,material_phase,phase_localConstitution
|
|
implicit none
|
|
|
|
integer(pInt) g,i,e, gMax,iMax,eMax, myNgrains
|
|
|
|
gMax = homogenization_maxNgrains
|
|
iMax = mesh_maxNips
|
|
eMax = mesh_NcpElems
|
|
|
|
allocate(crystallite_Fe(3,3,gMax,iMax,eMax)); crystallite_Fe = 0.0_pReal
|
|
allocate(crystallite_Fp(3,3,gMax,iMax,eMax)); crystallite_Fp = 0.0_pReal
|
|
allocate(crystallite_Lp(3,3,gMax,iMax,eMax)); crystallite_Lp = 0.0_pReal
|
|
allocate(crystallite_F0(3,3,gMax,iMax,eMax)); crystallite_F0 = 0.0_pReal
|
|
allocate(crystallite_Fp0(3,3,gMax,iMax,eMax)); crystallite_Fp0 = 0.0_pReal
|
|
allocate(crystallite_Lp0(3,3,gMax,iMax,eMax)); crystallite_Lp0 = 0.0_pReal
|
|
allocate(crystallite_partionedF(3,3,gMax,iMax,eMax)); crystallite_partionedF0 = 0.0_pReal
|
|
allocate(crystallite_partionedF0(3,3,gMax,iMax,eMax)); crystallite_partionedF0 = 0.0_pReal
|
|
allocate(crystallite_partionedFp0(3,3,gMax,iMax,eMax)); crystallite_partionedFp0 = 0.0_pReal
|
|
allocate(crystallite_partionedLp0(3,3,gMax,iMax,eMax)); crystallite_partionedLp0 = 0.0_pReal
|
|
allocate(crystallite_subF(3,3,gMax,iMax,eMax)); crystallite_subF = 0.0_pReal
|
|
allocate(crystallite_subF0(3,3,gMax,iMax,eMax)); crystallite_subF0 = 0.0_pReal
|
|
allocate(crystallite_subFp0(3,3,gMax,iMax,eMax)); crystallite_subFp0 = 0.0_pReal
|
|
allocate(crystallite_subLp0(3,3,gMax,iMax,eMax)); crystallite_subLp0 = 0.0_pReal
|
|
allocate(crystallite_P(3,3,gMax,iMax,eMax)); crystallite_P = 0.0_pReal
|
|
allocate(crystallite_Tstar_v(6,gMax,iMax,eMax)); crystallite_Tstar_v = 0.0_pReal
|
|
allocate(crystallite_dPdF(3,3,3,3,gMax,iMax,eMax)); crystallite_dPdF = 0.0_pReal
|
|
allocate(crystallite_fallbackdPdF(3,3,3,3,gMax,iMax,eMax)); crystallite_fallbackdPdF = 0.0_pReal
|
|
allocate(crystallite_dt(gMax,iMax,eMax)); crystallite_dt = 0.0_pReal
|
|
allocate(crystallite_subdt(gMax,iMax,eMax)); crystallite_subdt = 0.0_pReal
|
|
allocate(crystallite_subFrac(gMax,iMax,eMax)); crystallite_subFrac = 0.0_pReal
|
|
allocate(crystallite_subStep(gMax,iMax,eMax)); crystallite_subStep = 0.0_pReal
|
|
allocate(crystallite_Temperature(gMax,iMax,eMax)); crystallite_Temperature = 0.0_pReal
|
|
allocate(crystallite_localConstitution(gMax,iMax,eMax));
|
|
allocate(crystallite_requested(gMax,iMax,eMax)); crystallite_requested = .false.
|
|
allocate(crystallite_onTrack(gMax,iMax,eMax)); crystallite_onTrack = .false.
|
|
allocate(crystallite_converged(gMax,iMax,eMax)); crystallite_converged = .true.
|
|
|
|
!$OMP PARALLEL DO
|
|
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over all cp elements
|
|
myNgrains = homogenization_Ngrains(mesh_element(3,e))
|
|
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element
|
|
do g = 1,myNgrains
|
|
crystallite_Fp0(:,:,g,i,e) = math_EulerToR(material_EulerAngles(:,g,i,e)) ! plastic def gradient reflects init orientation
|
|
crystallite_F0(:,:,g,i,e) = math_I3
|
|
crystallite_partionedFp0(:,:,g,i,e) = crystallite_Fp0(:,:,g,i,e)
|
|
crystallite_partionedF0(:,:,g,i,e) = crystallite_F0(:,:,g,i,e)
|
|
crystallite_partionedF(:,:,g,i,e) = crystallite_F0(:,:,g,i,e)
|
|
crystallite_requested(g,i,e) = .true.
|
|
crystallite_localConstitution(g,i,e) = phase_localConstitution(material_phase(g,i,e))
|
|
enddo
|
|
enddo
|
|
enddo
|
|
!$OMP END PARALLEL DO
|
|
|
|
call crystallite_stressAndItsTangent(.true.) ! request elastic answers
|
|
crystallite_fallbackdPdF = crystallite_dPdF ! use initial elastic stiffness as fallback
|
|
|
|
! *** Output to MARC output file ***
|
|
!$OMP CRITICAL (write2out)
|
|
write(6,*)
|
|
write(6,*) '<<<+- crystallite init -+>>>'
|
|
write(6,*)
|
|
write(6,'(a32,x,7(i5,x))') 'crystallite_Nresults: ', crystallite_Nresults
|
|
write(6,'(a32,x,7(i5,x))') 'crystallite_Fe: ', shape(crystallite_Fe)
|
|
write(6,'(a32,x,7(i5,x))') 'crystallite_Fp: ', shape(crystallite_Fp)
|
|
write(6,'(a32,x,7(i5,x))') 'crystallite_Lp: ', shape(crystallite_Lp)
|
|
write(6,'(a32,x,7(i5,x))') 'crystallite_F0: ', shape(crystallite_F0)
|
|
write(6,'(a32,x,7(i5,x))') 'crystallite_Fp0: ', shape(crystallite_Fp0)
|
|
write(6,'(a32,x,7(i5,x))') 'crystallite_Lp0: ', shape(crystallite_Lp0)
|
|
write(6,'(a32,x,7(i5,x))') 'crystallite_partionedF: ', shape(crystallite_partionedF)
|
|
write(6,'(a32,x,7(i5,x))') 'crystallite_partionedF0: ', shape(crystallite_partionedF0)
|
|
write(6,'(a32,x,7(i5,x))') 'crystallite_partionedFp0: ', shape(crystallite_partionedFp0)
|
|
write(6,'(a32,x,7(i5,x))') 'crystallite_partionedLp0: ', shape(crystallite_partionedLp0)
|
|
write(6,'(a32,x,7(i5,x))') 'crystallite_subF: ', shape(crystallite_subF)
|
|
write(6,'(a32,x,7(i5,x))') 'crystallite_subF0: ', shape(crystallite_subF0)
|
|
write(6,'(a32,x,7(i5,x))') 'crystallite_subFp0: ', shape(crystallite_subFp0)
|
|
write(6,'(a32,x,7(i5,x))') 'crystallite_subLp0: ', shape(crystallite_subLp0)
|
|
write(6,'(a32,x,7(i5,x))') 'crystallite_P: ', shape(crystallite_P)
|
|
write(6,'(a32,x,7(i5,x))') 'crystallite_Tstar_v: ', shape(crystallite_Tstar_v)
|
|
write(6,'(a32,x,7(i5,x))') 'crystallite_dPdF: ', shape(crystallite_dPdF)
|
|
write(6,'(a32,x,7(i5,x))') 'crystallite_fallbackdPdF: ', shape(crystallite_fallbackdPdF)
|
|
write(6,'(a32,x,7(i5,x))') 'crystallite_dt: ', shape(crystallite_dt)
|
|
write(6,'(a32,x,7(i5,x))') 'crystallite_subdt: ', shape(crystallite_subdt)
|
|
write(6,'(a32,x,7(i5,x))') 'crystallite_subFrac: ', shape(crystallite_subFrac)
|
|
write(6,'(a32,x,7(i5,x))') 'crystallite_subStep: ', shape(crystallite_subStep)
|
|
write(6,'(a32,x,7(i5,x))') 'crystallite_Temperature: ', shape(crystallite_Temperature)
|
|
write(6,'(a32,x,7(i5,x))') 'crystallite_localConstitution: ', shape(crystallite_localConstitution)
|
|
write(6,'(a32,x,7(i5,x))') 'crystallite_requested: ', shape(crystallite_requested)
|
|
write(6,'(a32,x,7(i5,x))') 'crystallite_onTrack: ', shape(crystallite_onTrack)
|
|
write(6,'(a32,x,7(i5,x))') 'crystallite_converged: ', shape(crystallite_converged)
|
|
write(6,*)
|
|
write(6,*) 'Number of non-local grains: ',count(.not. crystallite_localConstitution)
|
|
call flush(6)
|
|
!$OMP END CRITICAL (write2out)
|
|
|
|
call debug_info()
|
|
call debug_reset()
|
|
|
|
return
|
|
|
|
end subroutine
|
|
|
|
|
|
!********************************************************************
|
|
! calculate stress (P) and tangent (dPdF) for crystallites
|
|
!********************************************************************
|
|
subroutine crystallite_stressAndItsTangent(updateJaco)
|
|
|
|
use prec, only: pInt,pReal,subStepMin,nCryst
|
|
use debug
|
|
use IO, only: IO_warning
|
|
use math
|
|
use FEsolving, only: FEsolving_execElem, FEsolving_execIP
|
|
use mesh, only: mesh_element
|
|
use material, only: homogenization_Ngrains
|
|
use constitutive
|
|
implicit none
|
|
|
|
logical, intent(in) :: updateJaco
|
|
real(pReal), dimension(3,3) :: invFp,Fe_guess,PK2,myF,myFp,myFe,myLp,myP
|
|
real(pReal), dimension(constitutive_maxSizeState) :: myState
|
|
integer(pInt) crystallite_Niteration
|
|
integer(pInt) g,i,e,k,l, myNgrains, mySizeState
|
|
logical, dimension(2) :: doneAndHappy
|
|
|
|
! ------ initialize to starting condition ------
|
|
|
|
write (6,*)
|
|
write (6,*) 'Crystallite request from Materialpoint'
|
|
write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedF0 of 1 8 1',crystallite_partionedF0(1:3,:,1,8,1)
|
|
write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedFp0 of 1 8 1',crystallite_partionedFp0(1:3,:,1,8,1)
|
|
write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedF of 1 8 1',crystallite_partionedF(1:3,:,1,8,1)
|
|
write (6,'(a,/,3(3(f12.7,x)/))') 'crystallite_partionedLp0 of 1 8 1',crystallite_partionedLp0(1:3,:,1,8,1)
|
|
|
|
|
|
!$OMP PARALLEL DO
|
|
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
|
|
myNgrains = homogenization_Ngrains(mesh_element(3,e))
|
|
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
|
|
do g = 1,myNgrains
|
|
if (crystallite_requested(g,i,e)) then ! initialize restoration point of ...
|
|
constitutive_subState0(g,i,e)%p = constitutive_partionedState0(g,i,e)%p ! ...microstructure
|
|
crystallite_subFp0(:,:,g,i,e) = crystallite_partionedFp0(:,:,g,i,e) ! ...plastic def grad
|
|
crystallite_subLp0(:,:,g,i,e) = crystallite_partionedLp0(:,:,g,i,e) ! ...plastic velocity grad
|
|
crystallite_subF0(:,:,g,i,e) = crystallite_partionedF0(:,:,g,i,e) ! ...def grad
|
|
|
|
crystallite_subFrac(g,i,e) = 0.0_pReal
|
|
crystallite_subStep(g,i,e) = 2.0_pReal
|
|
crystallite_onTrack(g,i,e) = .true.
|
|
crystallite_converged(g,i,e) = .false. ! pretend failed step of twice the required size
|
|
endif
|
|
enddo
|
|
enddo
|
|
enddo
|
|
!$OMP END PARALLEL DO
|
|
|
|
! ------ cutback loop ------
|
|
|
|
do while (any(crystallite_subStep(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) > subStepMin)) ! cutback loop for crystallites
|
|
write(6,*)
|
|
write(6,*) 'entering cutback at crystallite_stress'
|
|
if (any(.not. crystallite_converged(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) .and. & ! any non-converged grain
|
|
.not. crystallite_localConstitution(:,:,FEsolving_execELem(1):FEsolving_execElem(2))) ) & ! has non-local constitution?
|
|
crystallite_converged(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) = &
|
|
crystallite_converged(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) .and. &
|
|
crystallite_localConstitution(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) ! reset non-local grains' convergence status
|
|
!$OMP PARALLEL DO
|
|
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
|
|
myNgrains = homogenization_Ngrains(mesh_element(3,e))
|
|
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
|
|
do g = 1,myNgrains
|
|
if (crystallite_converged(g,i,e)) then
|
|
crystallite_subFrac(g,i,e) = crystallite_subFrac(g,i,e) + crystallite_subStep(g,i,e)
|
|
crystallite_subStep(g,i,e) = min(1.0_pReal-crystallite_subFrac(g,i,e), 2.0_pReal * crystallite_subStep(g,i,e))
|
|
if (crystallite_subStep(g,i,e) > subStepMin) then
|
|
crystallite_subF0(:,:,g,i,e) = crystallite_subF(:,:,g,i,e) ! wind forward...
|
|
crystallite_subFp0(:,:,g,i,e) = crystallite_Fp(:,:,g,i,e) ! ...plastic def grad
|
|
crystallite_subLp0(:,:,g,i,e) = crystallite_Lp(:,:,g,i,e) ! ...plastic velocity gradient
|
|
constitutive_subState0(g,i,e)%p = constitutive_state(g,i,e)%p ! ...microstructure
|
|
endif
|
|
else
|
|
crystallite_subStep(g,i,e) = 0.5_pReal * crystallite_subStep(g,i,e) ! cut step in half and restore...
|
|
crystallite_Fp(:,:,g,i,e) = crystallite_subFp0(:,:,g,i,e) ! ...plastic def grad
|
|
crystallite_Lp(:,:,g,i,e) = crystallite_subLp0(:,:,g,i,e) ! ...plastic velocity grad
|
|
constitutive_state(g,i,e)%p = constitutive_subState0(g,i,e)%p ! ...microstructure
|
|
endif
|
|
|
|
crystallite_onTrack(g,i,e) = crystallite_subStep(g,i,e) > subStepMin ! still on track or already done (beyond repair)
|
|
if (crystallite_onTrack(g,i,e)) then ! specify task (according to substep)
|
|
crystallite_subF(:,:,g,i,e) = crystallite_subF0(:,:,g,i,e) + &
|
|
crystallite_subStep(g,i,e) * &
|
|
(crystallite_partionedF(:,:,g,i,e) - crystallite_partionedF0(:,:,g,i,e))
|
|
crystallite_subdt(g,i,e) = crystallite_subStep(g,i,e) * crystallite_dt(g,i,e)
|
|
crystallite_converged(g,i,e) = .false. ! start out non-converged
|
|
endif
|
|
enddo
|
|
enddo
|
|
enddo
|
|
!$OMP END PARALLEL DO
|
|
|
|
! ------ convergence loop for stress and state ------
|
|
|
|
crystallite_Niteration = 0_pInt
|
|
|
|
do while (any( crystallite_requested(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) &
|
|
.and. crystallite_onTrack(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) &
|
|
.and. .not. crystallite_converged(:,:,FEsolving_execELem(1):FEsolving_execElem(2)) &
|
|
) .and. crystallite_Niteration < nCryst) ! convergence loop for crystallite
|
|
crystallite_Niteration = crystallite_Niteration + 1
|
|
|
|
|
|
! --+>> stress integration <<+--
|
|
!
|
|
! incrementing by crystallite_subdt
|
|
! based on crystallite_subF0,.._subFp0,.._subLp0
|
|
! constitutive_state is internally interpolated with .._subState0
|
|
! to account for substepping within _integrateStress
|
|
! results in crystallite_Fp,.._Lp
|
|
|
|
!$OMP PARALLEL DO
|
|
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
|
|
myNgrains = homogenization_Ngrains(mesh_element(3,e))
|
|
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
|
|
do g = 1,myNgrains
|
|
if (crystallite_requested(g,i,e) .and. &
|
|
crystallite_onTrack(g,i,e)) & ! all undone crystallites
|
|
crystallite_onTrack(g,i,e) = crystallite_integrateStress(g,i,e)
|
|
enddo
|
|
enddo
|
|
enddo
|
|
!$OMP END PARALLEL DO
|
|
|
|
write(6,'(i2,x,a10,x,16(l,x))') crystallite_Niteration,'cryst_onT',crystallite_onTrack
|
|
write (6,'(a,/,3(3(f12.7,x)/))') 'Lp of 1 8 1',crystallite_Lp(1:3,:,1,8,1)
|
|
|
|
! --+>> state integration <<+--
|
|
!
|
|
! incrementing by crystallite_subdt
|
|
! based on constitutive_subState0
|
|
! results in constitutive_state
|
|
|
|
!$OMP PARALLEL DO
|
|
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
|
|
myNgrains = homogenization_Ngrains(mesh_element(3,e))
|
|
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
|
|
do g = 1,myNgrains
|
|
if (crystallite_requested(g,i,e) .and. &
|
|
crystallite_onTrack(g,i,e)) then ! all undone crystallites
|
|
crystallite_converged(g,i,e) = crystallite_updateState(g,i,e)
|
|
if (crystallite_converged(g,i,e)) then
|
|
!$OMP CRITICAL (distributionState)
|
|
debug_StateLoopDistribution(crystallite_Niteration) = &
|
|
debug_StateLoopDistribution(crystallite_Niteration) + 1
|
|
!$OMP END CRITICAL (distributionState)
|
|
endif
|
|
endif
|
|
enddo
|
|
enddo
|
|
enddo
|
|
!$OMP END PARALLEL DO
|
|
|
|
enddo ! crystallite convergence loop
|
|
|
|
write(6,*)
|
|
write(6,'(a10,x,16(f6.4,x))') 'cryst_frac',crystallite_subFrac
|
|
write(6,'(a10,x,16(f6.4,x))') 'cryst_step',crystallite_subStep
|
|
write(6,'(a10,x,16(l,x))') 'cryst_req',crystallite_requested
|
|
write(6,'(a10,x,16(l,x))') 'cryst_onT',crystallite_onTrack
|
|
write(6,'(a10,x,16(l,x))') 'cryst_cvg',crystallite_converged
|
|
write(6,'(a10,x,16(e8.3,x))') 'cryst_dt',crystallite_subdt
|
|
enddo ! cutback loop
|
|
|
|
|
|
|
|
! ------ check for non-converged crystallites ------
|
|
|
|
!$OMP PARALLEL DO
|
|
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
|
|
myNgrains = homogenization_Ngrains(mesh_element(3,e))
|
|
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
|
|
do g = 1,myNgrains
|
|
if (.not. crystallite_converged(g,i,e)) then ! respond fully elastically
|
|
call IO_warning(600,e,i,g)
|
|
invFp = math_inv3x3(crystallite_partionedFp0(:,:,g,i,e))
|
|
Fe_guess = math_mul33x33(crystallite_partionedF(:,:,g,i,e),invFp)
|
|
PK2 = math_Mandel6to33( &
|
|
math_mul66x6( &
|
|
0.5_pReal*constitutive_homogenizedC(g,i,e), &
|
|
math_Mandel33to6( &
|
|
math_mul33x33(transpose(Fe_guess),Fe_guess) - math_I3 &
|
|
) &
|
|
) &
|
|
)
|
|
crystallite_P(:,:,g,i,e) = math_mul33x33(Fe_guess,math_mul33x33(PK2,transpose(invFp)))
|
|
endif
|
|
enddo
|
|
enddo
|
|
enddo
|
|
!$OMP END PARALLEL DO
|
|
|
|
! ------ stiffness calculation ------
|
|
|
|
if(updateJaco) then ! Jacobian required
|
|
if (debugger) then
|
|
!$OMP CRITICAL (write2out)
|
|
write (6,*) 'Jacobian calc'
|
|
write(6,'(a10,x,16(f6.4,x))') 'cryst_dt',crystallite_subdt
|
|
!$OMP END CRITICAL (write2out)
|
|
endif
|
|
|
|
!$OMP PARALLEL DO
|
|
do e = FEsolving_execElem(1),FEsolving_execElem(2) ! iterate over elements to be processed
|
|
myNgrains = homogenization_Ngrains(mesh_element(3,e))
|
|
do i = FEsolving_execIP(1,e),FEsolving_execIP(2,e) ! iterate over IPs of this element to be processed
|
|
do g = 1,myNgrains
|
|
if (crystallite_converged(g,i,e)) then ! grain converged in above iteration
|
|
mySizeState = constitutive_sizeState(g,i,e) ! number of state variables for this grain
|
|
myState(1:mySizeState) = constitutive_state(g,i,e)%p ! remember unperturbed, converged state...
|
|
myF = crystallite_subF(:,:,g,i,e) ! ... and kinematics
|
|
myFp = crystallite_Fp(:,:,g,i,e)
|
|
myFe = crystallite_Fe(:,:,g,i,e)
|
|
myLp = crystallite_Lp(:,:,g,i,e)
|
|
myP = crystallite_P(:,:,g,i,e)
|
|
do k = 1,3 ! perturbation...
|
|
do l = 1,3 ! ...components
|
|
crystallite_subF(:,:,g,i,e) = myF ! initialize perturbed F to match converged
|
|
crystallite_subF(k,l,g,i,e) = crystallite_subF(k,l,g,i,e) + pert_Fg ! perturb single component
|
|
doneAndHappy = .false.
|
|
crystallite_Niteration = 0_pInt
|
|
do while(.not. doneAndHappy(1) .and. crystallite_Niteration < nCryst) ! keep cycling until done (though unhappy)
|
|
crystallite_Niteration = crystallite_Niteration + 1_pInt
|
|
doneAndHappy = crystallite_integrateStress(g,i,e) ! stress of perturbed situation (overwrites _P,_Tstar_v,_Fp,_Lp,_Fe)
|
|
if (doneAndHappy(2)) & ! happy stress allows for state update
|
|
doneAndHappy = crystallite_updateState(g,i,e)
|
|
end do
|
|
if (doneAndHappy(2)) & ! happy outcome warrants stiffness update
|
|
crystallite_dPdF(:,:,k,l,g,i,e) = (crystallite_p(:,:,g,i,e) - myP)/pert_Fg ! tangent dP_ij/dFg_kl
|
|
!$OMP CRITICAL (out)
|
|
debug_StiffnessStateLoopDistribution(crystallite_Niteration) = &
|
|
debug_StiffnessstateLoopDistribution(crystallite_Niteration) + 1
|
|
!$OMP END CRITICAL (out)
|
|
end do
|
|
end do
|
|
constitutive_state(g,i,e)%p = myState ! restore unperturbed, converged state...
|
|
crystallite_Fp(:,:,g,i,e) = myFp ! ... and kinematics
|
|
crystallite_Fe(:,:,g,i,e) = myFe
|
|
crystallite_Lp(:,:,g,i,e) = myLp
|
|
crystallite_P(:,:,g,i,e) = myP
|
|
else ! grain has not converged
|
|
crystallite_dPdF(:,:,:,:,g,i,e) = crystallite_fallbackdPdF(:,:,:,:,g,i,e) ! use fallback
|
|
endif
|
|
enddo
|
|
enddo
|
|
enddo
|
|
!$OMP END PARALLEL DO
|
|
endif
|
|
|
|
end subroutine
|
|
|
|
|
|
|
|
|
|
!********************************************************************
|
|
! update the internal state of the constitutive law
|
|
! and tell whether state has converged
|
|
!********************************************************************
|
|
function crystallite_updateState(&
|
|
g,& ! grain number
|
|
i,& ! integration point number
|
|
e & ! element number
|
|
)
|
|
use prec, only: pReal,pInt,rTol_crystalliteState
|
|
use constitutive, only: constitutive_dotState,constitutive_sizeDotState,&
|
|
constitutive_subState0,constitutive_state
|
|
use debug
|
|
|
|
logical crystallite_updateState
|
|
|
|
integer(pLongInt) tick,tock,tickrate,maxticks
|
|
integer(pInt) g,i,e,mySize
|
|
real(pReal), dimension(6) :: Tstar_v
|
|
real(pReal), dimension(constitutive_sizeDotState(g,i,e)) :: residuum
|
|
|
|
mySize = constitutive_sizeDotState(g,i,e)
|
|
call system_clock(count=tick,count_rate=tickrate,count_max=maxticks)
|
|
residuum = constitutive_state(g,i,e)%p(1:mySize) - constitutive_subState0(g,i,e)%p(1:mySize) - &
|
|
crystallite_subdt(g,i,e)*&
|
|
constitutive_dotState(crystallite_Tstar_v(:,g,i,e),crystallite_Temperature(g,i,e),g,i,e) ! residuum from evolution of microstructure
|
|
call system_clock(count=tock,count_rate=tickrate,count_max=maxticks)
|
|
debug_cumDotStateCalls = debug_cumDotStateCalls + 1_pInt
|
|
debug_cumDotStateTicks = debug_cumDotStateTicks + tock-tick
|
|
if (tock < tick) debug_cumDotStateTicks = debug_cumDotStateTicks + maxticks
|
|
constitutive_state(g,i,e)%p(1:mySize) = constitutive_state(g,i,e)%p(1:mySize) - residuum ! update of microstructure
|
|
crystallite_updateState = maxval(abs(residuum/constitutive_state(g,i,e)%p(1:mySize)),&
|
|
constitutive_state(g,i,e)%p(1:mySize) /= 0.0_pReal) < rTol_crystalliteState
|
|
return
|
|
|
|
end function
|
|
|
|
|
|
|
|
!***********************************************************************
|
|
!*** calculation of stress (P), stiffness (dPdF), ***
|
|
!*** and announcement of any ***
|
|
!*** acceleration of the Newton-Raphson correction ***
|
|
!***********************************************************************
|
|
function crystallite_integrateStress(&
|
|
g,& ! grain number
|
|
i,& ! integration point number
|
|
e) ! element number
|
|
|
|
use prec, only: pReal,pInt,subStepMin
|
|
use debug
|
|
use constitutive, only: constitutive_state,constitutive_subState0,constitutive_sizeState
|
|
use math
|
|
|
|
implicit none
|
|
|
|
integer(pInt), intent(in) :: e,i,g
|
|
integer(pInt) mySize, Niteration
|
|
real(pReal) dt_aim,subFrac,subStep,det
|
|
real(pReal), dimension(3,3) :: inv
|
|
real(pReal), dimension(3,3) :: Fg_current,Fg_aim,deltaFg
|
|
real(pReal), dimension(3,3) :: Fp_current,Fp_new
|
|
real(pReal), dimension(3,3) :: Fe_current,Fe_new
|
|
real(pReal), dimension(constitutive_sizeState(g,i,e)) :: interpolatedState
|
|
|
|
logical crystallite_integrateStress ! still on track?
|
|
|
|
mySize = constitutive_sizeState(g,i,e)
|
|
deltaFg = crystallite_subF(:,:,g,i,e) - crystallite_subF0(:,:,g,i,e)
|
|
Fg_current = crystallite_subF0(:,:,g,i,e) ! initialize to start of inc
|
|
Fp_current = crystallite_subFp0(:,:,g,i,e)
|
|
Fe_current = math_mul33x33(Fg_current,math_inv3x3(Fp_current))
|
|
subFrac = 0.0_pReal
|
|
subStep = 1.0_pReal
|
|
Niteration = 0_pInt
|
|
crystallite_integrateStress = .false. ! be pessimisitc
|
|
|
|
! begin the cutback loop
|
|
do while (subStep > subStepMin) ! continue until finished or too much cut backing
|
|
Niteration = Niteration + 1
|
|
Fg_aim = Fg_current + subStep*deltaFg ! aim for Fg
|
|
dt_aim = subStep*crystallite_subdt(g,i,e) ! aim for dt
|
|
debugger = (g == 1 .and. i == 8 .and. e == 1)
|
|
call TimeIntegration(crystallite_integrateStress,&
|
|
crystallite_Lp(:,:,g,i,e),crystallite_Fp(:,:,g,i,e),crystallite_Fe(:,:,g,i,e),&
|
|
crystallite_Tstar_v(:,g,i,e),crystallite_P(:,:,g,i,e), &
|
|
Fg_aim,Fp_current,crystallite_Temperature(g,i,e),&
|
|
( subFrac+subStep)*constitutive_state(g,i,e)%p(1:mySize) + &
|
|
(1.0_pReal-subFrac-subStep)*constitutive_subState0(g,i,e)%p(1:mySize),& ! interpolated state
|
|
dt_aim,g,i,e)
|
|
if (crystallite_integrateStress) then ! happy with time integration
|
|
if (e == 1 .and. i == 8 .and. g == 1) then
|
|
write(6,*) '*** winding forward in IntegrateStress ***'
|
|
write(6,*) subFrac, subFrac+subStep
|
|
write (6,'(a,/,3(3(f12.7,x)/))') 'Lp of 1 8 1',crystallite_Lp(1:3,:,1,8,1)
|
|
write (6,'(a,/,3(3(f12.7,x)/))') 'Fp of 1 8 1',crystallite_Fp(1:3,:,1,8,1)
|
|
endif
|
|
Fg_current = Fg_aim ! wind forward
|
|
Fe_current = crystallite_Fe(:,:,g,i,e)
|
|
Fp_current = crystallite_Fp(:,:,g,i,e)
|
|
subFrac = subFrac + subStep
|
|
subStep = min(1.0_pReal-subFrac, subStep*2.0_pReal) ! accelerate
|
|
else ! time integration encountered trouble
|
|
subStep = 0.5_pReal * subStep ! cut time step in half
|
|
crystallite_Lp(:,:,g,i,e) = 0.5_pReal*(crystallite_Lp(:,:,g,i,e) + & ! interpolate Lp and L
|
|
(math_I3 - math_mul33x33(Fp_current,&
|
|
math_mul33x33(math_inv3x3(Fg_aim),Fe_current)))/dt_aim)
|
|
endif
|
|
enddo ! potential substepping
|
|
|
|
!$OMP CRITICAL (distributionStress)
|
|
debug_StressLoopDistribution(Niteration) = debug_StressLoopDistribution(Niteration) + 1
|
|
!$OMP END CRITICAL (distributionStress)
|
|
|
|
return ! with final happyness
|
|
|
|
end function
|
|
|
|
|
|
!***********************************************************************
|
|
!*** fully-implicit two-level time integration ***
|
|
!*** based on a residuum in Lp and intermediate ***
|
|
!*** acceleration of the Newton-Raphson correction ***
|
|
!***********************************************************************
|
|
subroutine TimeIntegration(&
|
|
happy,& ! return status
|
|
Lpguess,& ! guess of plastic velocity gradient
|
|
Fp_new,& ! new plastic deformation gradient
|
|
Fe_new,& ! new "elastic" deformation gradient
|
|
Tstar_v,& ! Stress vector
|
|
P,& ! 1st PK stress (taken as initial guess if /= 0)
|
|
!
|
|
Fg_new,& ! new total def gradient
|
|
Fp_old,& ! former plastic def gradient
|
|
Temperature,& ! temperature
|
|
state,& ! microstructural state
|
|
dt,& ! time increment
|
|
grain,& ! grain number
|
|
ip,& ! integration point number
|
|
cp_en & ! element number
|
|
)
|
|
|
|
use prec
|
|
use debug
|
|
use mesh, only: mesh_element
|
|
use constitutive, only: constitutive_microstructure,constitutive_homogenizedC,constitutive_LpAndItsTangent,&
|
|
constitutive_sizeState
|
|
use math
|
|
use IO
|
|
implicit none
|
|
|
|
logical, intent(out) :: happy
|
|
real(pReal), dimension(3,3), intent(inout) :: Lpguess
|
|
real(pReal), dimension(3,3), intent(out) :: Fp_new,Fe_new,P
|
|
real(pReal), dimension(6), intent(out) :: Tstar_v
|
|
real(pReal), dimension(3,3), intent(in) :: Fg_new,Fp_old
|
|
real(pReal), intent(in) :: Temperature,dt
|
|
integer(pInt), intent(in) :: cp_en, ip, grain
|
|
real(pReal), dimension(constitutive_sizeState(grain,ip,cp_en)), intent(in) :: state
|
|
|
|
logical error
|
|
integer(pInt) Niteration,dummy, i,j,k,l,m,n
|
|
integer(pLongInt) tick,tock,tickrate,maxticks
|
|
real(pReal) p_hydro,det, leapfrog,maxleap
|
|
real(pReal), dimension(3,3,3,3) :: C
|
|
real(pReal), dimension(9,9) :: dLp,dTdLp,dRdLp,invdRdLp,eye2
|
|
real(pReal), dimension(6,6) :: C_66
|
|
real(pReal), dimension(3,3) :: invFp_new,invFp_old,Lp,Lpguess_old,Rinner,Rinner_old,A,B,BT,AB,BTA
|
|
|
|
happy = .false. ! be pessimistic
|
|
eye2 = math_identity2nd(9)
|
|
|
|
invFp_old = math_inv3x3(Fp_old) ! inversion of Fp_old
|
|
if (all(invFp_old == 0.0_pReal)) return ! failed
|
|
! write (6,'(a,/,3(3(f12.7,x)/))') 'Fp old',Fp_old
|
|
! write (6,'(a,/,3(3(f12.7,x)/))') 'Fp old inv',invFp_old
|
|
|
|
A = math_mul33x33(transpose(invFp_old), math_mul33x33(transpose(Fg_new),math_mul33x33(Fg_new,invFp_old)))
|
|
|
|
!$OMP CRITICAL (write2out)
|
|
! if (debugger) write (6,'(a,/,3(3(f12.7,x)/))') 'Fg to be calculated',Fg_new
|
|
!$OMP END CRITICAL (write2out)
|
|
|
|
call constitutive_microstructure(Temperature,grain,ip,cp_en)
|
|
C_66 = constitutive_homogenizedC(grain,ip,cp_en)
|
|
C = math_Mandel66to3333(C_66) ! 4th rank elasticity tensor
|
|
|
|
Niteration = 0_pInt
|
|
leapfrog = 1.0_pReal ! correction as suggested by invdRdLp-step
|
|
maxleap = 1024.0_pReal ! preassign maximum acceleration level
|
|
|
|
Lpguess_old = Lpguess ! consider present Lpguess good (i.e. worth remembering)
|
|
|
|
Inner: do ! inner iteration: Lp
|
|
Niteration = Niteration + 1
|
|
if (Niteration > nLp) then ! too many loops required
|
|
Lpguess = Lpguess_old ! do not trust the last update but resort to former one
|
|
return
|
|
endif
|
|
|
|
! write(6,*) 'iteration',Niteration
|
|
! write (6,'(a,/,3(3(f12.7,x)/))') 'Lpguess',Lpguess
|
|
|
|
B = math_i3 - dt*Lpguess
|
|
BT = transpose(B)
|
|
AB = math_mul33x33(A,B)
|
|
BTA = math_mul33x33(BT,A)
|
|
Tstar_v = 0.5_pReal*math_mul66x6(C_66,math_mandel33to6(math_mul33x33(BT,AB)-math_I3))
|
|
p_hydro = (Tstar_v(1) + Tstar_v(2) + Tstar_v(3))/3.0_pReal
|
|
forall(i=1:3) Tstar_v(i) = Tstar_v(i) - p_hydro ! subtract hydrostatic pressure
|
|
! write (6,'(a,/,6(f12.7,x))') 'Tstar',Tstar_v
|
|
call system_clock(count=tick,count_rate=tickrate,count_max=maxticks)
|
|
call constitutive_LpAndItsTangent(Lp,dLp, Tstar_v,Temperature,grain,ip,cp_en)
|
|
call system_clock(count=tock,count_rate=tickrate,count_max=maxticks)
|
|
debug_cumLpCalls = debug_cumLpCalls + 1_pInt
|
|
debug_cumLpTicks = debug_cumLpTicks + tock-tick
|
|
if (tock < tick) debug_cumLpTicks = debug_cumLpTicks + maxticks
|
|
|
|
Rinner = Lpguess - Lp ! update current residuum
|
|
! write (6,'(a,/,3(3(f12.7,x)/))') 'Lp',Lp
|
|
! write (6,'(a,/,3(3(f12.7,x)/))') 'Residuum',Rinner
|
|
|
|
if (.not.(any(Rinner/=Rinner)) .and. & ! exclude any NaN in residuum
|
|
( ( maxval(abs(Rinner)) < aTol_crystalliteStress) .or. & ! below abs tol .or.
|
|
( any(abs(dt*Lpguess) > relevantStrain) .and. & ! worth checking? .and.
|
|
maxval(abs(Rinner/Lpguess),abs(dt*Lpguess) > relevantStrain) < rTol_crystalliteStress & ! below rel tol
|
|
) &
|
|
) &
|
|
) &
|
|
exit Inner ! convergence
|
|
!
|
|
! check for acceleration/deceleration in Newton--Raphson correction
|
|
!
|
|
if (any(Rinner/=Rinner) .and. & ! NaN occured at regular speed
|
|
leapfrog == 1.0) then
|
|
Lpguess = Lpguess_old ! restore known good guess and croak for cutback
|
|
return
|
|
|
|
elseif (leapfrog > 1.0_pReal .and. & ! at fast pace ?
|
|
(sum(Rinner*Rinner) > sum(Rinner_old*Rinner_old) .or. & ! worse residuum
|
|
sum(Rinner*Rinner_old) < 0.0_pReal) .or. & ! residuum changed sign (overshoot)
|
|
any(Rinner/=Rinner) ) then ! NaN
|
|
maxleap = 0.5_pReal * leapfrog ! limit next acceleration
|
|
leapfrog = 1.0_pReal ! grinding halt
|
|
|
|
else ! better residuum
|
|
dTdLp = 0.0_pReal ! calc dT/dLp
|
|
forall (i=1:3,j=1:3,k=1:3,l=1:3,m=1:3,n=1:3) &
|
|
dTdLp(3*(i-1)+j,3*(k-1)+l) = dTdLp(3*(i-1)+j,3*(k-1)+l) + &
|
|
C(i,j,l,n)*AB(k,n)+C(i,j,m,l)*BTA(m,k)
|
|
dTdLp = -0.5_pReal*dt*dTdLp
|
|
dRdLp = eye2 - math_mul99x99(dLp,dTdLp) ! calc dR/dLp
|
|
invdRdLp = 0.0_pReal
|
|
call math_invert(9,dRdLp,invdRdLp,dummy,error) ! invert dR/dLp --> dLp/dR
|
|
if (error) then
|
|
if (debugger) then
|
|
!$OMP CRITICAL (write2out)
|
|
write (6,*) 'inversion dR/dLp failed',grain,ip,cp_en
|
|
! write (6,'(a,/,9(9(e9.3,x)/))') 'dRdLp', dRdLp(1:9,:)
|
|
! write (6,'(a,/,3(4(f9.3,x)/))') 'state_new / MPa',state/1e6_pReal
|
|
write (6,'(a,/,3(3(f12.7,x)/))') 'Lpguess',Lpguess(1:3,:)
|
|
write (6,'(a,/,3(3(e12.7,x)/))') 'Lp',Lp(1:3,:)
|
|
write (6,'(a,/,6(f9.3,x))') 'Tstar / MPa',Tstar_v/1e6_pReal
|
|
!$OMP END CRITICAL (write2out)
|
|
endif
|
|
return
|
|
endif
|
|
|
|
Rinner_old = Rinner ! remember current residuum
|
|
Lpguess_old = Lpguess ! remember current Lp guess
|
|
if (Niteration > 1 .and. leapfrog < maxleap) leapfrog = 2.0_pReal * leapfrog ! accelerate if ok
|
|
endif
|
|
|
|
Lpguess = Lpguess_old ! start from current guess
|
|
Rinner = Rinner_old ! use current residuum
|
|
forall (i=1:3,j=1:3,k=1:3,l=1:3) & ! leapfrog to updated Lpguess
|
|
Lpguess(i,j) = Lpguess(i,j) - leapfrog*invdRdLp(3*(i-1)+j,3*(k-1)+l)*Rinner(k,l)
|
|
enddo Inner
|
|
|
|
!$OMP CRITICAL (distributionLp)
|
|
debug_LpLoopDistribution(Niteration) = debug_LpLoopDistribution(Niteration) + 1
|
|
!$OMP END CRITICAL (distributionLp)
|
|
invFp_new = math_mul33x33(invFp_old,B)
|
|
if (debugger) then
|
|
write (6,'(a,x,f10.6,/,3(3(f12.7,x)/))') 'Lp(guess)',dt,Lpguess(1:3,:)
|
|
write (6,'(a,/,3(3(f12.7,x)/))') 'invFp_old',invFp_old(1:3,:)
|
|
write (6,'(a,/,3(3(f12.7,x)/))') 'B',B(1:3,:)
|
|
write (6,'(a,/,3(3(f12.7,x)/))') 'invFp_new',invFp_new(1:3,:)
|
|
endif
|
|
call math_invert3x3(invFp_new,Fp_new,det,error)
|
|
if (debugger) then
|
|
write (6,'(a,/,3(3(f12.7,x)/))') 'invFp_new',invFp_new(1:3,:)
|
|
write (6,'(a,/,3(3(f12.7,x)/))') 'Fp_new',Fp_new(1:3,:)
|
|
write (6,'(a,x,l,x,a,f10.6)') 'with inversion error:',error,'and determinant:',det
|
|
endif
|
|
if (error) return ! inversion failed
|
|
|
|
Fp_new = Fp_new*det**(1.0_pReal/3.0_pReal) ! regularize Fp by det = det(InvFp_new) !!
|
|
forall (i=1:3) Tstar_v(i) = Tstar_v(i) + p_hydro ! add hydrostatic component back
|
|
Fe_new = math_mul33x33(Fg_new,invFp_new) ! calc resulting Fe
|
|
P = math_mul33x33(Fe_new,math_mul33x33(math_Mandel6to33(Tstar_v),transpose(invFp_new))) ! first PK stress
|
|
|
|
happy = .true. ! now smile...
|
|
|
|
return
|
|
|
|
end subroutine
|
|
|
|
|
|
|
|
!********************************************************************
|
|
! return results of particular grain
|
|
!********************************************************************
|
|
function crystallite_postResults(&
|
|
Tstar_v,& ! stress
|
|
Temperature, & ! temperature
|
|
dt,& ! time increment
|
|
g,& ! grain number
|
|
i,& ! integration point number
|
|
e & ! element number
|
|
)
|
|
|
|
use prec, only: pInt,pReal
|
|
use math, only: math_pDecomposition,math_RtoEuler, inDeg
|
|
use IO, only: IO_warning
|
|
use material, only: material_phase,material_volfrac
|
|
use constitutive, only: constitutive_sizePostResults, constitutive_postResults
|
|
implicit none
|
|
|
|
integer(pInt), intent(in) :: g,i,e
|
|
real(pReal), intent(in) :: Temperature,dt
|
|
real(pReal), dimension(6), intent(in) :: Tstar_v
|
|
real(pReal), dimension(3,3) :: U,R
|
|
logical error
|
|
|
|
real(pReal), dimension(crystallite_Nresults + constitutive_sizePostResults(g,i,e)) :: crystallite_postResults
|
|
|
|
if (crystallite_Nresults >= 2) then
|
|
crystallite_postResults(1) = material_phase(g,i,e)
|
|
crystallite_postResults(2) = material_volfrac(g,i,e)
|
|
endif
|
|
if (crystallite_Nresults >= 5) then
|
|
call math_pDecomposition(crystallite_Fe(:,:,g,i,e),U,R,error) ! polar decomposition of Fe
|
|
if (error) then
|
|
call IO_warning(650,e,i,g)
|
|
crystallite_postResults(3:5) = (/400.0,400.0,400.0/) ! fake orientation
|
|
else
|
|
crystallite_postResults(3:5) = math_RtoEuler(transpose(R))*inDeg ! orientation
|
|
endif
|
|
endif
|
|
|
|
crystallite_postResults(crystallite_Nresults+1:crystallite_Nresults+constitutive_sizePostResults(g,i,e)) = &
|
|
constitutive_postResults(Tstar_v,Temperature,dt,g,i,e)
|
|
return
|
|
|
|
end function
|
|
|
|
|
|
END MODULE
|
|
!############################################################## |