Merge branch 'no-crystallite_subF' into 'development'
No crystallite_subF See merge request damask/DAMASK!313
This commit is contained in:
Sharan Roongta
commit
584c7cc3a6
|
|
@ -5,8 +5,8 @@ homogenization:
|
|||
|
||||
phase:
|
||||
Aluminum:
|
||||
mechanics:
|
||||
lattice: cF
|
||||
mechanics:
|
||||
output: [F, P, F_e, F_p, L_p]
|
||||
elasticity: {C_11: 106.75e9, C_12: 60.41e9, C_44: 28.34e9, type: hooke}
|
||||
plasticity:
|
||||
|
|
|
|||
|
|
@ -48,7 +48,6 @@ module constitutive
|
|||
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_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
|
||||
|
|
@ -60,9 +59,8 @@ module constitutive
|
|||
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
|
||||
crystallite_partitionedF0, & !< def grad at start of homog inc
|
||||
crystallite_F !< def grad to be reached at end of homog inc
|
||||
|
||||
type :: tTensorContainer
|
||||
real(pReal), dimension(:,:,:), allocatable :: data
|
||||
|
|
@ -179,6 +177,14 @@ module constitutive
|
|||
module subroutine constitutive_mech_forward
|
||||
end subroutine constitutive_mech_forward
|
||||
|
||||
module subroutine mech_restore(ip,el,includeL)
|
||||
integer, intent(in) :: &
|
||||
ip, &
|
||||
el
|
||||
logical, intent(in) :: &
|
||||
includeL
|
||||
end subroutine mech_restore
|
||||
|
||||
! == cleaned:end ===================================================================================
|
||||
|
||||
module function crystallite_stress(dt,co,ip,el) result(converged_)
|
||||
|
|
@ -392,7 +398,6 @@ module constitutive
|
|||
crystallite_restartRead, &
|
||||
constitutive_initializeRestorationPoints, &
|
||||
constitutive_windForward, &
|
||||
crystallite_restore, &
|
||||
PLASTICITY_UNDEFINED_ID, &
|
||||
PLASTICITY_NONE_ID, &
|
||||
PLASTICITY_ISOTROPIC_ID, &
|
||||
|
|
@ -734,6 +739,7 @@ subroutine constitutive_allocateState(state, &
|
|||
sizeDotState, &
|
||||
sizeDeltaState
|
||||
|
||||
|
||||
state%sizeState = sizeState
|
||||
state%sizeDotState = sizeDotState
|
||||
state%sizeDeltaState = sizeDeltaState
|
||||
|
|
@ -756,22 +762,27 @@ end subroutine constitutive_allocateState
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief Restore data after homog cutback.
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine constitutive_restore(ip,el)
|
||||
subroutine constitutive_restore(ip,el,includeL)
|
||||
|
||||
logical, intent(in) :: includeL
|
||||
integer, intent(in) :: &
|
||||
ip, & !< integration point number
|
||||
el !< element number
|
||||
|
||||
integer :: &
|
||||
co, & !< constituent number
|
||||
s
|
||||
so
|
||||
|
||||
|
||||
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))
|
||||
do so = 1, phase_Nsources(material_phaseAt(co,el))
|
||||
sourceState(material_phaseAt(co,el))%p(so)%state( :,material_phasememberAt(co,ip,el)) = &
|
||||
sourceState(material_phaseAt(co,el))%p(so)%partitionedState0(:,material_phasememberAt(co,ip,el))
|
||||
enddo
|
||||
enddo
|
||||
|
||||
call mech_restore(ip,el,includeL)
|
||||
|
||||
end subroutine constitutive_restore
|
||||
|
||||
|
||||
|
|
@ -783,7 +794,7 @@ subroutine constitutive_forward
|
|||
|
||||
integer :: i, j
|
||||
|
||||
crystallite_F0 = crystallite_partitionedF
|
||||
crystallite_F0 = crystallite_F
|
||||
crystallite_Lp0 = crystallite_Lp
|
||||
crystallite_S0 = crystallite_S
|
||||
|
||||
|
|
@ -837,6 +848,7 @@ subroutine crystallite_init
|
|||
iMax, & !< maximum number of integration points
|
||||
eMax !< maximum number of elements
|
||||
|
||||
|
||||
class(tNode), pointer :: &
|
||||
num_crystallite, &
|
||||
debug_crystallite, & ! pointer to debug options for crystallite
|
||||
|
|
@ -854,7 +866,7 @@ subroutine crystallite_init
|
|||
iMax = discretization_nIPs
|
||||
eMax = discretization_Nelems
|
||||
|
||||
allocate(crystallite_partitionedF(3,3,cMax,iMax,eMax),source=0.0_pReal)
|
||||
allocate(crystallite_F(3,3,cMax,iMax,eMax),source=0.0_pReal)
|
||||
|
||||
allocate(crystallite_S0, &
|
||||
crystallite_F0,crystallite_Lp0, &
|
||||
|
|
@ -863,14 +875,12 @@ subroutine crystallite_init
|
|||
crystallite_partitionedLp0, &
|
||||
crystallite_S,crystallite_P, &
|
||||
crystallite_Fe,crystallite_Lp, &
|
||||
crystallite_subF, &
|
||||
crystallite_subFp0,crystallite_subFi0, &
|
||||
source = crystallite_partitionedF)
|
||||
source = crystallite_F)
|
||||
|
||||
allocate(crystallite_subdt(cMax,iMax,eMax),source=0.0_pReal)
|
||||
allocate(crystallite_orientation(cMax,iMax,eMax))
|
||||
|
||||
|
||||
num_crystallite => config_numerics%get('crystallite',defaultVal=emptyDict)
|
||||
|
||||
num%subStepMinCryst = num_crystallite%get_asFloat ('subStepMin', defaultVal=1.0e-3_pReal)
|
||||
|
|
@ -933,8 +943,8 @@ subroutine crystallite_init
|
|||
flush(IO_STDOUT)
|
||||
|
||||
!$OMP PARALLEL DO PRIVATE(ph,me)
|
||||
do el = 1, size(material_phaseMemberAt,3)
|
||||
do ip = 1, size(material_phaseMemberAt,2); do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
|
||||
do el = 1, size(material_phaseMemberAt,3); do ip = 1, size(material_phaseMemberAt,2)
|
||||
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
|
||||
|
||||
ph = material_phaseAt(co,el)
|
||||
me = material_phaseMemberAt(co,ip,el)
|
||||
|
|
@ -953,12 +963,12 @@ subroutine crystallite_init
|
|||
constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me) = constitutive_mech_Fi0(ph)%data(1:3,1:3,me)
|
||||
constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me) = constitutive_mech_Fp0(ph)%data(1:3,1:3,me)
|
||||
|
||||
enddo; enddo
|
||||
enddo
|
||||
enddo; enddo
|
||||
!$OMP END PARALLEL DO
|
||||
|
||||
crystallite_partitionedF0 = crystallite_F0
|
||||
crystallite_partitionedF = crystallite_F0
|
||||
crystallite_F = crystallite_F0
|
||||
|
||||
|
||||
!$OMP PARALLEL DO PRIVATE(ph,me)
|
||||
|
|
@ -978,9 +988,6 @@ subroutine crystallite_init
|
|||
end subroutine crystallite_init
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief Backup data for homog cutback.
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
@ -991,7 +998,7 @@ subroutine constitutive_initializeRestorationPoints(ip,el)
|
|||
el !< element number
|
||||
integer :: &
|
||||
co, & !< constituent number
|
||||
s,ph, me
|
||||
so,ph, me
|
||||
|
||||
do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
|
||||
ph = material_phaseAt(co,el)
|
||||
|
|
@ -1002,9 +1009,9 @@ subroutine constitutive_initializeRestorationPoints(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))
|
||||
do so = 1, phase_Nsources(material_phaseAt(co,el))
|
||||
sourceState(material_phaseAt(co,el))%p(so)%partitionedState0(:,material_phasememberAt(co,ip,el)) = &
|
||||
sourceState(material_phaseAt(co,el))%p(so)%state0( :,material_phasememberAt(co,ip,el))
|
||||
enddo
|
||||
enddo
|
||||
|
||||
|
|
@ -1019,57 +1026,28 @@ subroutine constitutive_windForward(ip,el)
|
|||
integer, intent(in) :: &
|
||||
ip, & !< integration point number
|
||||
el !< element number
|
||||
|
||||
integer :: &
|
||||
co, & !< constituent number
|
||||
s, ph, me
|
||||
so, 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_partitionedF0 (1:3,1:3,co,ip,el) = crystallite_F (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)
|
||||
do so = 1, phase_Nsources(material_phaseAt(co,el))
|
||||
sourceState(ph)%p(so)%partitionedState0(:,me) = sourceState(ph)%p(so)%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_partitionedLi0(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_partitionedFp0(p)%data(1:3,1:3,m)
|
||||
constitutive_mech_Fi(p)%data(1:3,1:3,m) = constitutive_mech_partitionedFi0(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).
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
@ -1080,13 +1058,13 @@ function crystallite_stressTangent(co,ip,el) result(dPdF)
|
|||
co, & !< counter in constituent loop
|
||||
ip, & !< counter in integration point loop
|
||||
el !< counter in element loop
|
||||
|
||||
integer :: &
|
||||
o, &
|
||||
p, ph, me
|
||||
|
||||
real(pReal), dimension(3,3) :: devNull, &
|
||||
invSubFp0,invSubFi0,invFp,invFi, &
|
||||
temp_33_1, temp_33_2, temp_33_3, temp_33_4
|
||||
temp_33_1, temp_33_2, temp_33_3
|
||||
real(pReal), dimension(3,3,3,3) :: dSdFe, &
|
||||
dSdF, &
|
||||
dSdFi, &
|
||||
|
|
@ -1102,6 +1080,7 @@ function crystallite_stressTangent(co,ip,el) result(dPdF)
|
|||
real(pReal), dimension(9,9):: temp_99
|
||||
logical :: error
|
||||
|
||||
|
||||
ph = material_phaseAt(co,el)
|
||||
me = material_phaseMemberAt(co,ip,el)
|
||||
|
||||
|
|
@ -1149,8 +1128,8 @@ function crystallite_stressTangent(co,ip,el) result(dPdF)
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
! 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)
|
||||
temp_33_2 = matmul(crystallite_F(1:3,1:3,co,ip,el),invSubFp0)
|
||||
temp_33_3 = matmul(matmul(crystallite_F(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)
|
||||
|
|
@ -1180,21 +1159,20 @@ function crystallite_stressTangent(co,ip,el) result(dPdF)
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
! 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))
|
||||
temp_33_2 = matmul(crystallite_F(1:3,1:3,co,ip,el),invFp)
|
||||
temp_33_3 = matmul(temp_33_2,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)
|
||||
dPdF(p,1:3,p,1:3) = transpose(matmul(invFp,temp_33_1))
|
||||
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), &
|
||||
+ matmul(matmul(crystallite_F(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)), &
|
||||
+ matmul(matmul(temp_33_2,dSdF(1:3,1:3,p,o)), &
|
||||
transpose(invFp)) &
|
||||
+ matmul(temp_33_4,transpose(dFpinvdF(1:3,1:3,p,o)))
|
||||
+ matmul(temp_33_3,transpose(dFpinvdF(1:3,1:3,p,o)))
|
||||
enddo; enddo
|
||||
|
||||
end function crystallite_stressTangent
|
||||
|
|
@ -1237,7 +1215,7 @@ function crystallite_push33ToRef(co,ip,el, tensor33)
|
|||
|
||||
|
||||
T = matmul(material_orientation0(co,ip,el)%asMatrix(), & ! ToDo: initial orientation correct?
|
||||
transpose(math_inv33(crystallite_subF(1:3,1:3,co,ip,el))))
|
||||
transpose(math_inv33(crystallite_F(1:3,1:3,co,ip,el))))
|
||||
crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T))
|
||||
|
||||
end function crystallite_push33ToRef
|
||||
|
|
@ -1247,8 +1225,9 @@ end function crystallite_push33ToRef
|
|||
!> @brief integrate stress, state with adaptive 1st order explicit Euler method
|
||||
!> using Fixed Point Iteration to adapt the stepsize
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
function integrateSourceState(co,ip,el) result(broken)
|
||||
function integrateSourceState(dt,co,ip,el) result(broken)
|
||||
|
||||
real(pReal), intent(in) :: dt
|
||||
integer, intent(in) :: &
|
||||
el, & !< element index in element loop
|
||||
ip, & !< integration point index in ip loop
|
||||
|
|
@ -1281,8 +1260,7 @@ function integrateSourceState(co,ip,el) result(broken)
|
|||
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)
|
||||
+ sourceState(ph)%p(so)%dotState (1:size_so(so),me) * dt
|
||||
source_dotState(1:size_so(so),2,so) = 0.0_pReal
|
||||
enddo
|
||||
|
||||
|
|
@ -1305,7 +1283,7 @@ function integrateSourceState(co,ip,el) result(broken)
|
|||
+ 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)%dotState (1:size_so(so),me) * dt
|
||||
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))
|
||||
converged_ = converged_ .and. converged(r(1:size_so(so)), &
|
||||
|
|
@ -1371,7 +1349,7 @@ end function converged
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine crystallite_restartWrite
|
||||
|
||||
integer :: i
|
||||
integer :: ph
|
||||
integer(HID_T) :: fileHandle, groupHandle
|
||||
character(len=pStringLen) :: fileName, datasetName
|
||||
|
||||
|
|
@ -1380,27 +1358,27 @@ subroutine crystallite_restartWrite
|
|||
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_F,'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)
|
||||
do ph = 1,size(material_name_phase)
|
||||
write(datasetName,'(i0,a)') ph,'_omega'
|
||||
call HDF5_write(groupHandle,plasticState(ph)%state,datasetName)
|
||||
write(datasetName,'(i0,a)') ph,'_F_i'
|
||||
call HDF5_write(groupHandle,constitutive_mech_Fi(ph)%data,datasetName)
|
||||
write(datasetName,'(i0,a)') ph,'_L_i'
|
||||
call HDF5_write(groupHandle,constitutive_mech_Li(ph)%data,datasetName)
|
||||
write(datasetName,'(i0,a)') ph,'_F_p'
|
||||
call HDF5_write(groupHandle,constitutive_mech_Fp(ph)%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)
|
||||
do ph = 1, size(material_name_homogenization)
|
||||
write(datasetName,'(i0,a)') ph,'_omega'
|
||||
call HDF5_write(groupHandle,homogState(ph)%state,datasetName)
|
||||
enddo
|
||||
call HDF5_closeGroup(groupHandle)
|
||||
|
||||
|
|
@ -1415,7 +1393,7 @@ end subroutine crystallite_restartWrite
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine crystallite_restartRead
|
||||
|
||||
integer :: i
|
||||
integer :: ph
|
||||
integer(HID_T) :: fileHandle, groupHandle
|
||||
character(len=pStringLen) :: fileName, datasetName
|
||||
|
||||
|
|
@ -1429,22 +1407,22 @@ subroutine crystallite_restartRead
|
|||
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)
|
||||
do ph = 1,size(material_name_phase)
|
||||
write(datasetName,'(i0,a)') ph,'_omega'
|
||||
call HDF5_read(groupHandle,plasticState(ph)%state0,datasetName)
|
||||
write(datasetName,'(i0,a)') ph,'_F_i'
|
||||
call HDF5_read(groupHandle,constitutive_mech_Fi0(ph)%data,datasetName)
|
||||
write(datasetName,'(i0,a)') ph,'_L_i'
|
||||
call HDF5_read(groupHandle,constitutive_mech_Li0(ph)%data,datasetName)
|
||||
write(datasetName,'(i0,a)') ph,'_F_p'
|
||||
call HDF5_read(groupHandle,constitutive_mech_Fp0(ph)%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)
|
||||
do ph = 1,size(material_name_homogenization)
|
||||
write(datasetName,'(i0,a)') ph,'_omega'
|
||||
call HDF5_read(groupHandle,homogState(ph)%state0,datasetName)
|
||||
enddo
|
||||
call HDF5_closeGroup(groupHandle)
|
||||
|
||||
|
|
|
|||
|
|
@ -800,7 +800,7 @@ function integrateStress(F,Delta_t,co,ip,el) result(broken)
|
|||
|
||||
broken = .true.
|
||||
|
||||
call constitutive_plastic_dependentState(crystallite_partitionedF(1:3,1:3,co,ip,el),co,ip,el)
|
||||
call constitutive_plastic_dependentState(crystallite_F(1:3,1:3,co,ip,el),co,ip,el)
|
||||
|
||||
ph = material_phaseAt(co,el)
|
||||
me = material_phaseMemberAt(co,ip,el)
|
||||
|
|
@ -959,19 +959,21 @@ function integrateStateFPI(F_0,F,Delta_t,co,ip,el) result(broken)
|
|||
el, & !< element index in element loop
|
||||
ip, & !< integration point index in ip loop
|
||||
co !< grain index in grain loop
|
||||
logical :: &
|
||||
broken
|
||||
|
||||
integer :: &
|
||||
NiterationState, & !< number of iterations in state loop
|
||||
ph, &
|
||||
me, &
|
||||
size_pl
|
||||
sizeDotState
|
||||
real(pReal) :: &
|
||||
zeta
|
||||
real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: &
|
||||
r ! state residuum
|
||||
real(pReal), dimension(constitutive_plasticity_maxSizeDotState,2) :: &
|
||||
plastic_dotState
|
||||
logical :: &
|
||||
broken
|
||||
dotState
|
||||
|
||||
|
||||
ph = material_phaseAt(co,el)
|
||||
me = material_phaseMemberAt(co,ip,el)
|
||||
|
|
@ -979,15 +981,15 @@ function integrateStateFPI(F_0,F,Delta_t,co,ip,el) result(broken)
|
|||
broken = mech_collectDotState(Delta_t, co,ip,el,ph,me)
|
||||
if(broken) return
|
||||
|
||||
size_pl = plasticState(ph)%sizeDotState
|
||||
plasticState(ph)%state(1:size_pl,me) = plasticState(ph)%subState0(1:size_pl,me) &
|
||||
+ plasticState(ph)%dotState (1:size_pl,me) * Delta_t
|
||||
plastic_dotState(1:size_pl,2) = 0.0_pReal
|
||||
sizeDotState = plasticState(ph)%sizeDotState
|
||||
plasticState(ph)%state(1:sizeDotState,me) = plasticState(ph)%subState0(1:sizeDotState,me) &
|
||||
+ plasticState(ph)%dotState (1:sizeDotState,me) * Delta_t
|
||||
dotState(1:sizeDotState,2) = 0.0_pReal
|
||||
|
||||
iteration: do NiterationState = 1, num%nState
|
||||
|
||||
if(nIterationState > 1) plastic_dotState(1:size_pl,2) = plastic_dotState(1:size_pl,1)
|
||||
plastic_dotState(1:size_pl,1) = plasticState(ph)%dotState(:,me)
|
||||
if(nIterationState > 1) dotState(1:sizeDotState,2) = dotState(1:sizeDotState,1)
|
||||
dotState(1:sizeDotState,1) = plasticState(ph)%dotState(:,me)
|
||||
|
||||
broken = integrateStress(F,Delta_t,co,ip,el)
|
||||
if(broken) exit iteration
|
||||
|
|
@ -995,16 +997,16 @@ function integrateStateFPI(F_0,F,Delta_t,co,ip,el) result(broken)
|
|||
broken = mech_collectDotState(Delta_t, co,ip,el,ph,me)
|
||||
if(broken) exit iteration
|
||||
|
||||
zeta = damper(plasticState(ph)%dotState(:,me),plastic_dotState(1:size_pl,1),&
|
||||
plastic_dotState(1:size_pl,2))
|
||||
zeta = damper(plasticState(ph)%dotState(:,me),dotState(1:sizeDotState,1),&
|
||||
dotState(1:sizeDotState,2))
|
||||
plasticState(ph)%dotState(:,me) = plasticState(ph)%dotState(:,me) * zeta &
|
||||
+ plastic_dotState(1:size_pl,1) * (1.0_pReal - zeta)
|
||||
r(1:size_pl) = plasticState(ph)%state (1:size_pl,me) &
|
||||
- plasticState(ph)%subState0(1:size_pl,me) &
|
||||
- plasticState(ph)%dotState (1:size_pl,me) * Delta_t
|
||||
plasticState(ph)%state(1:size_pl,me) = plasticState(ph)%state(1:size_pl,me) &
|
||||
- r(1:size_pl)
|
||||
if (converged(r(1:size_pl),plasticState(ph)%state(1:size_pl,me),plasticState(ph)%atol(1:size_pl))) then
|
||||
+ dotState(1:sizeDotState,1) * (1.0_pReal - zeta)
|
||||
r(1:sizeDotState) = plasticState(ph)%state (1:sizeDotState,me) &
|
||||
- plasticState(ph)%subState0(1:sizeDotState,me) &
|
||||
- plasticState(ph)%dotState (1:sizeDotState,me) * Delta_t
|
||||
plasticState(ph)%state(1:sizeDotState,me) = plasticState(ph)%state(1:sizeDotState,me) &
|
||||
- r(1:sizeDotState)
|
||||
if (converged(r(1:sizeDotState),plasticState(ph)%state(1:sizeDotState,me),plasticState(ph)%atol(1:sizeDotState))) then
|
||||
broken = constitutive_deltaState(crystallite_S(1:3,1:3,co,ip,el), &
|
||||
constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el,ph,me)
|
||||
exit iteration
|
||||
|
|
@ -1012,6 +1014,7 @@ function integrateStateFPI(F_0,F,Delta_t,co,ip,el) result(broken)
|
|||
|
||||
enddo iteration
|
||||
|
||||
|
||||
contains
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
@ -1048,12 +1051,14 @@ function integrateStateEuler(F_0,F,Delta_t,co,ip,el) result(broken)
|
|||
el, & !< element index in element loop
|
||||
ip, & !< integration point index in ip loop
|
||||
co !< grain index in grain loop
|
||||
logical :: &
|
||||
broken
|
||||
|
||||
integer :: &
|
||||
ph, &
|
||||
me, &
|
||||
sizeDotState
|
||||
logical :: &
|
||||
broken
|
||||
|
||||
|
||||
ph = material_phaseAt(co,el)
|
||||
me = material_phaseMemberAt(co,ip,el)
|
||||
|
|
@ -1085,13 +1090,13 @@ function integrateStateAdaptiveEuler(F_0,F,Delta_t,co,ip,el) result(broken)
|
|||
el, & !< element index in element loop
|
||||
ip, & !< integration point index in ip loop
|
||||
co !< grain index in grain loop
|
||||
logical :: &
|
||||
broken
|
||||
|
||||
integer :: &
|
||||
ph, &
|
||||
me, &
|
||||
sizeDotState
|
||||
logical :: &
|
||||
broken
|
||||
|
||||
real(pReal), dimension(constitutive_plasticity_maxSizeDotState) :: residuum_plastic
|
||||
|
||||
|
||||
|
|
@ -1145,6 +1150,7 @@ function integrateStateRK4(F_0,F,Delta_t,co,ip,el) result(broken)
|
|||
real(pReal), dimension(4), parameter :: &
|
||||
B = [1.0_pReal/6.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/3.0_pReal, 1.0_pReal/6.0_pReal]
|
||||
|
||||
|
||||
broken = integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C)
|
||||
|
||||
end function integrateStateRK4
|
||||
|
|
@ -1178,6 +1184,7 @@ function integrateStateRKCK45(F_0,F,Delta_t,co,ip,el) result(broken)
|
|||
[2825.0_pReal/27648.0_pReal, .0_pReal, 18575.0_pReal/48384.0_pReal,&
|
||||
13525.0_pReal/55296.0_pReal, 277.0_pReal/14336.0_pReal, 1._pReal/4._pReal]
|
||||
|
||||
|
||||
broken = integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB)
|
||||
|
||||
end function integrateStateRKCK45
|
||||
|
|
@ -1215,18 +1222,18 @@ function integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB) result(broken)
|
|||
broken = mech_collectDotState(Delta_t,co,ip,el,ph,me)
|
||||
if(broken) return
|
||||
|
||||
do stage = 1, size(A,1)
|
||||
sizeDotState = plasticState(ph)%sizeDotState
|
||||
|
||||
do stage = 1, size(A,1)
|
||||
|
||||
plastic_RKdotState(1:sizeDotState,stage) = plasticState(ph)%dotState(:,me)
|
||||
plasticState(ph)%dotState(:,me) = A(1,stage) * plastic_RKdotState(1:sizeDotState,1)
|
||||
|
||||
do n = 2, stage
|
||||
sizeDotState = plasticState(ph)%sizeDotState
|
||||
plasticState(ph)%dotState(:,me) = plasticState(ph)%dotState(:,me) &
|
||||
+ A(n,stage) * plastic_RKdotState(1:sizeDotState,n)
|
||||
enddo
|
||||
|
||||
sizeDotState = plasticState(ph)%sizeDotState
|
||||
plasticState(ph)%state(1:sizeDotState,me) = plasticState(ph)%subState0(1:sizeDotState,me) &
|
||||
+ plasticState(ph)%dotState (1:sizeDotState,me) * Delta_t
|
||||
|
||||
|
|
@ -1239,7 +1246,6 @@ function integrateStateRK(F_0,F,Delta_t,co,ip,el,A,B,C,DB) result(broken)
|
|||
enddo
|
||||
if(broken) return
|
||||
|
||||
sizeDotState = plasticState(ph)%sizeDotState
|
||||
|
||||
plastic_RKdotState(1:sizeDotState,size(B)) = plasticState (ph)%dotState(:,me)
|
||||
plasticState(ph)%dotState(:,me) = matmul(plastic_RKdotState(1:sizeDotState,1:size(B)),B)
|
||||
|
|
@ -1282,7 +1288,7 @@ subroutine crystallite_results(group,ph)
|
|||
|
||||
select case (output_constituent(ph)%label(ou))
|
||||
case('F')
|
||||
selected_tensors = select_tensors(crystallite_partitionedF,ph)
|
||||
selected_tensors = select_tensors(crystallite_F,ph)
|
||||
call results_writeDataset(group//'/mechanics/',selected_tensors,output_constituent(ph)%label(ou),&
|
||||
'deformation gradient','1')
|
||||
case('F_e')
|
||||
|
|
@ -1482,25 +1488,24 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
|
|||
formerSubStep
|
||||
integer :: &
|
||||
NiterationCrystallite, & ! number of iterations in crystallite loop
|
||||
s, ph, me
|
||||
so, 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
|
||||
subF0
|
||||
subF0, &
|
||||
subF
|
||||
|
||||
|
||||
ph = material_phaseAt(co,el)
|
||||
me = material_phaseMemberAt(co,ip,el)
|
||||
subLi0 = constitutive_mech_partitionedLi0(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))
|
||||
plasticState(ph)%subState0(:,me) = plasticState(ph)%partitionedState0(:,me)
|
||||
do so = 1, phase_Nsources(ph)
|
||||
sourceState(ph)%p(so)%subState0(:,me) = sourceState(ph)%p(so)%partitionedState0(:,me)
|
||||
enddo
|
||||
crystallite_subFp0(1:3,1:3,co,ip,el) = constitutive_mech_partitionedFp0(ph)%data(1:3,1:3,me)
|
||||
crystallite_subFi0(1:3,1:3,co,ip,el) = constitutive_mech_partitionedFi0(ph)%data(1:3,1:3,me)
|
||||
|
|
@ -1525,16 +1530,14 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
|
|||
todo = subStep > 0.0_pReal ! still time left to integrate on?
|
||||
|
||||
if (todo) then
|
||||
subF0 = crystallite_subF(1:3,1:3,co,ip,el)
|
||||
subF0 = subF
|
||||
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))
|
||||
plasticState(ph)%subState0(:,me) = plasticState(ph)%state(:,me)
|
||||
do so = 1, phase_Nsources(ph)
|
||||
sourceState(ph)%p(so)%subState0(:,me) = sourceState(ph)%p(so)%state(:,me)
|
||||
enddo
|
||||
endif
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
@ -1548,11 +1551,9 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
|
|||
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))
|
||||
plasticState(ph)%state(:,me) = plasticState(ph)%subState0(:,me)
|
||||
do so = 1, phase_Nsources(ph)
|
||||
sourceState(ph)%p(so)%state(:,me) = sourceState(ph)%p(so)%subState0(:,me)
|
||||
enddo
|
||||
|
||||
todo = subStep > num%subStepMinCryst ! still on track or already done (beyond repair)
|
||||
|
|
@ -1561,21 +1562,50 @@ module function crystallite_stress(dt,co,ip,el) result(converged_)
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
! prepare for integration
|
||||
if (todo) then
|
||||
crystallite_subF(1:3,1:3,co,ip,el) = subF0 &
|
||||
+ 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), &
|
||||
subF = subF0 &
|
||||
+ subStep * (crystallite_F(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(subF,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
|
||||
converged_ = .not. integrateState(subF0,crystallite_subF(1:3,1:3,co,ip,el),&
|
||||
crystallite_subdt(co,ip,el),co,ip,el)
|
||||
converged_ = converged_ .and. .not. integrateSourceState(co,ip,el)
|
||||
converged_ = .not. integrateState(subF0,subF,subStep * dt,co,ip,el)
|
||||
converged_ = converged_ .and. .not. integrateSourceState(subStep * dt,co,ip,el)
|
||||
endif
|
||||
|
||||
enddo cutbackLooping
|
||||
|
||||
end function crystallite_stress
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief Restore data after homog cutback.
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
module subroutine mech_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_partitionedLi0(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_partitionedFp0(p)%data(1:3,1:3,m)
|
||||
constitutive_mech_Fi(p)%data(1:3,1:3,m) = constitutive_mech_partitionedFi0(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 mech_restore
|
||||
|
||||
end submodule constitutive_mech
|
||||
|
||||
|
|
|
|||
|
|
@ -70,29 +70,22 @@ module homogenization
|
|||
end subroutine mech_homogenize
|
||||
|
||||
module subroutine mech_results(group_base,h)
|
||||
|
||||
character(len=*), intent(in) :: group_base
|
||||
integer, intent(in) :: h
|
||||
|
||||
end subroutine mech_results
|
||||
|
||||
! -------- ToDo ---------------------------------------------------------
|
||||
module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
|
||||
logical, dimension(2) :: mech_RGC_updateState
|
||||
real(pReal), dimension(:,:,:), intent(in) :: &
|
||||
P,& !< partitioned stresses
|
||||
F,& !< partitioned deformation gradients
|
||||
F0 !< partitioned initial deformation gradients
|
||||
real(pReal), dimension(:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses
|
||||
real(pReal), dimension(3,3), intent(in) :: avgF !< average F
|
||||
real(pReal), intent(in) :: dt !< time increment
|
||||
module function mech_updateState(subdt,subF,ip,el) result(doneAndHappy)
|
||||
real(pReal), intent(in) :: &
|
||||
subdt !< current time step
|
||||
real(pReal), intent(in), dimension(3,3) :: &
|
||||
subF
|
||||
integer, intent(in) :: &
|
||||
ip, & !< integration point number
|
||||
ip, & !< integration point
|
||||
el !< element number
|
||||
end function mech_RGC_updateState
|
||||
logical, dimension(2) :: doneAndHappy
|
||||
end function mech_updateState
|
||||
|
||||
end interface
|
||||
! -----------------------------------------------------------------------
|
||||
|
||||
public :: &
|
||||
homogenization_init, &
|
||||
|
|
@ -148,11 +141,10 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
|
|||
real(pReal), intent(in) :: dt !< time increment
|
||||
integer, dimension(2), intent(in) :: FEsolving_execElem, FEsolving_execIP
|
||||
integer :: &
|
||||
NiterationHomog, &
|
||||
NiterationMPstate, &
|
||||
ip, & !< integration point number
|
||||
el, & !< element number
|
||||
myNgrains, co, ce, ho
|
||||
myNgrains, co, ce, ho, me
|
||||
real(pReal) :: &
|
||||
subFrac, &
|
||||
subStep
|
||||
|
|
@ -162,12 +154,12 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
|
|||
doneAndHappy
|
||||
|
||||
|
||||
!$OMP PARALLEL DO PRIVATE(ce,ho,myNgrains,NiterationMPstate,NiterationHomog,subFrac,converged,subStep,doneAndHappy)
|
||||
!$OMP PARALLEL DO PRIVATE(ce,me,ho,myNgrains,NiterationMPstate,subFrac,converged,subStep,doneAndHappy)
|
||||
do el = FEsolving_execElem(1),FEsolving_execElem(2)
|
||||
ho = material_homogenizationAt(el)
|
||||
myNgrains = homogenization_Nconstituents(ho)
|
||||
do ip = FEsolving_execIP(1),FEsolving_execIP(2)
|
||||
|
||||
me = material_homogenizationMemberAt(ip,el)
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! initialize restoration points
|
||||
call constitutive_initializeRestorationPoints(ip,el)
|
||||
|
|
@ -177,15 +169,10 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
|
|||
subStep = 1.0_pReal/num%subStepSizeHomog ! ... larger then the requested calculation
|
||||
|
||||
if (homogState(ho)%sizeState > 0) &
|
||||
homogState(ho)%subState0(:,material_homogenizationMemberAt(ip,el)) = &
|
||||
homogState(ho)%State0( :,material_homogenizationMemberAt(ip,el))
|
||||
|
||||
homogState(ho)%subState0(:,me) = homogState(ho)%State0(:,me)
|
||||
if (damageState(ho)%sizeState > 0) &
|
||||
damageState(ho)%subState0(:,material_homogenizationMemberAt(ip,el)) = &
|
||||
damageState(ho)%State0( :,material_homogenizationMemberAt(ip,el))
|
||||
damageState(ho)%subState0(:,me) = damageState(ho)%State0(:,me)
|
||||
|
||||
|
||||
NiterationHomog = 0
|
||||
cutBackLooping: do while (.not. terminallyIll .and. subStep > num%subStepMinHomog)
|
||||
|
||||
if (converged) then
|
||||
|
|
@ -198,15 +185,12 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
|
|||
call constitutive_windForward(ip,el)
|
||||
|
||||
if(homogState(ho)%sizeState > 0) &
|
||||
homogState(ho)%subState0(:,material_homogenizationMemberAt(ip,el)) = &
|
||||
homogState(ho)%State (:,material_homogenizationMemberAt(ip,el))
|
||||
homogState(ho)%subState0(:,me) = homogState(ho)%State(:,me)
|
||||
if(damageState(ho)%sizeState > 0) &
|
||||
damageState(ho)%subState0(:,material_homogenizationMemberAt(ip,el)) = &
|
||||
damageState(ho)%State (:,material_homogenizationMemberAt(ip,el))
|
||||
damageState(ho)%subState0(:,me) = damageState(ho)%State(:,me)
|
||||
|
||||
endif steppingNeeded
|
||||
else
|
||||
if ( (myNgrains == 1 .and. subStep <= 1.0 ) .or. & ! single grain already tried internal subStepping in crystallite
|
||||
elseif ( (myNgrains == 1 .and. subStep <= 1.0 ) .or. & ! single grain already tried internal subStepping in crystallite
|
||||
num%subStepSizeHomog * subStep <= num%subStepMinHomog ) then ! would require too small subStep
|
||||
! cutback makes no sense
|
||||
if (.not. terminallyIll) & ! so first signals terminally ill...
|
||||
|
|
@ -215,16 +199,12 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
|
|||
else ! cutback makes sense
|
||||
subStep = num%subStepSizeHomog * subStep ! crystallite had severe trouble, so do a significant cutback
|
||||
|
||||
call crystallite_restore(ip,el,subStep < 1.0_pReal)
|
||||
call constitutive_restore(ip,el)
|
||||
call constitutive_restore(ip,el,subStep < 1.0_pReal)
|
||||
|
||||
if(homogState(ho)%sizeState > 0) &
|
||||
homogState(ho)%State( :,material_homogenizationMemberAt(ip,el)) = &
|
||||
homogState(ho)%subState0(:,material_homogenizationMemberAt(ip,el))
|
||||
homogState(ho)%State(:,me) = homogState(ho)%subState0(:,me)
|
||||
if(damageState(ho)%sizeState > 0) &
|
||||
damageState(ho)%State( :,material_homogenizationMemberAt(ip,el)) = &
|
||||
damageState(ho)%subState0(:,material_homogenizationMemberAt(ip,el))
|
||||
endif
|
||||
damageState(ho)%State(:,me) = damageState(ho)%subState0(:,me)
|
||||
endif
|
||||
|
||||
if (subStep > num%subStepMinHomog) doneAndHappy = [.false.,.true.]
|
||||
|
|
@ -253,7 +233,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
|
|||
doneAndHappy = [.true.,.false.]
|
||||
else
|
||||
ce = (el-1)*discretization_nIPs + ip
|
||||
doneAndHappy = updateState(dt*subStep, &
|
||||
doneAndHappy = mech_updateState(dt*subStep, &
|
||||
homogenization_F0(1:3,1:3,ce) &
|
||||
+ (homogenization_F(1:3,1:3,ce)-homogenization_F0(1:3,1:3,ce)) &
|
||||
*(subStep+subFrac), &
|
||||
|
|
@ -263,8 +243,6 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
|
|||
endif
|
||||
|
||||
enddo convergenceLooping
|
||||
NiterationHomog = NiterationHomog + 1
|
||||
|
||||
enddo cutBackLooping
|
||||
enddo
|
||||
enddo
|
||||
|
|
@ -290,74 +268,35 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
|
|||
end subroutine materialpoint_stressAndItsTangent
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief update the internal state of the homogenization scheme and tell whether "done" and
|
||||
!> "happy" with result
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
function updateState(subdt,subF,ip,el)
|
||||
|
||||
real(pReal), intent(in) :: &
|
||||
subdt !< current time step
|
||||
real(pReal), intent(in), dimension(3,3) :: &
|
||||
subF
|
||||
integer, intent(in) :: &
|
||||
ip, & !< integration point
|
||||
el !< element number
|
||||
integer :: c
|
||||
logical, dimension(2) :: updateState
|
||||
real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
|
||||
|
||||
updateState = .true.
|
||||
chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
|
||||
case (HOMOGENIZATION_RGC_ID) chosenHomogenization
|
||||
do c=1,homogenization_Nconstituents(material_homogenizationAt(el))
|
||||
dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el)
|
||||
enddo
|
||||
updateState = &
|
||||
updateState .and. &
|
||||
mech_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
|
||||
crystallite_partitionedF(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
|
||||
crystallite_partitionedF0(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el),&
|
||||
subF,&
|
||||
subdt, &
|
||||
dPdFs, &
|
||||
ip, &
|
||||
el)
|
||||
end select chosenHomogenization
|
||||
|
||||
end function updateState
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief writes homogenization results to HDF5 output file
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine homogenization_results
|
||||
use material, only: &
|
||||
material_homogenization_type => homogenization_type
|
||||
|
||||
integer :: p
|
||||
integer :: ho
|
||||
character(len=:), allocatable :: group_base,group
|
||||
|
||||
|
||||
call results_closeGroup(results_addGroup('current/homogenization/'))
|
||||
|
||||
do p=1,size(material_name_homogenization)
|
||||
group_base = 'current/homogenization/'//trim(material_name_homogenization(p))
|
||||
do ho=1,size(material_name_homogenization)
|
||||
group_base = 'current/homogenization/'//trim(material_name_homogenization(ho))
|
||||
call results_closeGroup(results_addGroup(group_base))
|
||||
|
||||
call mech_results(group_base,p)
|
||||
call mech_results(group_base,ho)
|
||||
|
||||
group = trim(group_base)//'/damage'
|
||||
call results_closeGroup(results_addGroup(group))
|
||||
select case(damage_type(p))
|
||||
select case(damage_type(ho))
|
||||
case(DAMAGE_NONLOCAL_ID)
|
||||
call damage_nonlocal_results(p,group)
|
||||
call damage_nonlocal_results(ho,group)
|
||||
end select
|
||||
|
||||
group = trim(group_base)//'/thermal'
|
||||
call results_closeGroup(results_addGroup(group))
|
||||
select case(thermal_type(p))
|
||||
select case(thermal_type(ho))
|
||||
case(THERMAL_CONDUCTION_ID)
|
||||
call thermal_conduction_results(p,group)
|
||||
call thermal_conduction_results(ho,group)
|
||||
end select
|
||||
|
||||
enddo
|
||||
|
|
@ -373,6 +312,7 @@ subroutine homogenization_forward
|
|||
|
||||
integer :: ho
|
||||
|
||||
|
||||
do ho = 1, size(material_name_homogenization)
|
||||
homogState (ho)%state0 = homogState (ho)%state
|
||||
damageState(ho)%state0 = damageState(ho)%state
|
||||
|
|
|
|||
|
|
@ -52,6 +52,21 @@ submodule(homogenization) homogenization_mech
|
|||
end subroutine mech_RGC_averageStressAndItsTangent
|
||||
|
||||
|
||||
module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) result(doneAndHappy)
|
||||
logical, dimension(2) :: doneAndHappy
|
||||
real(pReal), dimension(:,:,:), intent(in) :: &
|
||||
P,& !< partitioned stresses
|
||||
F,& !< partitioned deformation gradients
|
||||
F0 !< partitioned initial deformation gradients
|
||||
real(pReal), dimension(:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses
|
||||
real(pReal), dimension(3,3), intent(in) :: avgF !< average F
|
||||
real(pReal), intent(in) :: dt !< time increment
|
||||
integer, intent(in) :: &
|
||||
ip, & !< integration point number
|
||||
el !< element number
|
||||
end function mech_RGC_updateState
|
||||
|
||||
|
||||
module subroutine mech_RGC_results(instance,group)
|
||||
integer, intent(in) :: instance !< homogenization instance
|
||||
character(len=*), intent(in) :: group !< group name in HDF5 file
|
||||
|
|
@ -101,16 +116,16 @@ module subroutine mech_partition(subF,ip,el)
|
|||
chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
|
||||
|
||||
case (HOMOGENIZATION_NONE_ID) chosenHomogenization
|
||||
crystallite_partitionedF(1:3,1:3,1,ip,el) = subF
|
||||
crystallite_F(1:3,1:3,1,ip,el) = subF
|
||||
|
||||
case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
|
||||
call mech_isostrain_partitionDeformation(&
|
||||
crystallite_partitionedF(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
|
||||
crystallite_F(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
|
||||
subF)
|
||||
|
||||
case (HOMOGENIZATION_RGC_ID) chosenHomogenization
|
||||
call mech_RGC_partitionDeformation(&
|
||||
crystallite_partitionedF(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
|
||||
crystallite_F(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
|
||||
subF,&
|
||||
ip, &
|
||||
el)
|
||||
|
|
@ -166,6 +181,45 @@ module subroutine mech_homogenize(ip,el)
|
|||
end subroutine mech_homogenize
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief update the internal state of the homogenization scheme and tell whether "done" and
|
||||
!> "happy" with result
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
module function mech_updateState(subdt,subF,ip,el) result(doneAndHappy)
|
||||
|
||||
real(pReal), intent(in) :: &
|
||||
subdt !< current time step
|
||||
real(pReal), intent(in), dimension(3,3) :: &
|
||||
subF
|
||||
integer, intent(in) :: &
|
||||
ip, & !< integration point
|
||||
el !< element number
|
||||
logical, dimension(2) :: doneAndHappy
|
||||
|
||||
integer :: co
|
||||
real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
|
||||
|
||||
|
||||
if (homogenization_type(material_homogenizationAt(el)) == HOMOGENIZATION_RGC_ID) then
|
||||
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
|
||||
dPdFs(:,:,:,:,co) = crystallite_stressTangent(co,ip,el)
|
||||
enddo
|
||||
doneAndHappy = &
|
||||
mech_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
|
||||
crystallite_F(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
|
||||
crystallite_partitionedF0(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el),&
|
||||
subF,&
|
||||
subdt, &
|
||||
dPdFs, &
|
||||
ip, &
|
||||
el)
|
||||
else
|
||||
doneAndHappy = .true.
|
||||
endif
|
||||
|
||||
end function mech_updateState
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief Write results to file.
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
|
|||
|
|
@ -8,6 +8,7 @@
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
submodule(homogenization:homogenization_mech) homogenization_mech_RGC
|
||||
use rotations
|
||||
use lattice
|
||||
|
||||
type :: tParameters
|
||||
integer, dimension(:), allocatable :: &
|
||||
|
|
@ -242,7 +243,18 @@ end subroutine mech_RGC_partitionDeformation
|
|||
!> @brief update the internal state of the homogenization scheme and tell whether "done" and
|
||||
! "happy" with result
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
module procedure mech_RGC_updateState
|
||||
module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) result(doneAndHappy)
|
||||
logical, dimension(2) :: doneAndHappy
|
||||
real(pReal), dimension(:,:,:), intent(in) :: &
|
||||
P,& !< partitioned stresses
|
||||
F,& !< partitioned deformation gradients
|
||||
F0 !< partitioned initial deformation gradients
|
||||
real(pReal), dimension(:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses
|
||||
real(pReal), dimension(3,3), intent(in) :: avgF !< average F
|
||||
real(pReal), intent(in) :: dt !< time increment
|
||||
integer, intent(in) :: &
|
||||
ip, & !< integration point number
|
||||
el !< element number
|
||||
|
||||
integer, dimension(4) :: intFaceN,intFaceP,faceID
|
||||
integer, dimension(3) :: nGDim,iGr3N,iGr3P
|
||||
|
|
@ -256,7 +268,7 @@ module procedure mech_RGC_updateState
|
|||
real(pReal), dimension(:), allocatable :: resid,relax,p_relax,p_resid,drelax
|
||||
|
||||
zeroTimeStep: if(dEq0(dt)) then
|
||||
mech_RGC_updateState = .true. ! pretend everything is fine and return
|
||||
doneAndHappy = .true. ! pretend everything is fine and return
|
||||
return
|
||||
endif zeroTimeStep
|
||||
|
||||
|
|
@ -327,12 +339,12 @@ module procedure mech_RGC_updateState
|
|||
stresMax = maxval(abs(P)) ! get the maximum of first Piola-Kirchhoff (material) stress
|
||||
residMax = maxval(abs(tract)) ! get the maximum of the residual
|
||||
|
||||
mech_RGC_updateState = .false.
|
||||
doneAndHappy = .false.
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! If convergence reached => done and happy
|
||||
if (residMax < num%rtol*stresMax .or. residMax < num%atol) then
|
||||
mech_RGC_updateState = .true.
|
||||
doneAndHappy = .true.
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! compute/update the state for postResult, i.e., all energy densities computed by time-integration
|
||||
|
|
@ -354,7 +366,7 @@ module procedure mech_RGC_updateState
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
! if residual blows-up => done but unhappy
|
||||
elseif (residMax > num%relMax*stresMax .or. residMax > num%absMax) then ! try to restart when residual blows up exceeding maximum bound
|
||||
mech_RGC_updateState = [.true.,.false.] ! with direct cut-back
|
||||
doneAndHappy = [.true.,.false.] ! with direct cut-back
|
||||
return
|
||||
endif
|
||||
|
||||
|
|
@ -484,7 +496,7 @@ module procedure mech_RGC_updateState
|
|||
enddo; enddo
|
||||
stt%relaxationVector(:,of) = relax + drelax ! Updateing the state variable for the next iteration
|
||||
if (any(abs(drelax) > num%maxdRelax)) then ! Forcing cutback when the incremental change of relaxation vector becomes too large
|
||||
mech_RGC_updateState = [.true.,.false.]
|
||||
doneAndHappy = [.true.,.false.]
|
||||
!$OMP CRITICAL (write2out)
|
||||
print'(a,i3,a,i3,a)',' RGC_updateState: ip ',ip,' | el ',el,' enforces cutback'
|
||||
print'(a,e15.8)',' due to large relaxation change = ',maxval(abs(drelax))
|
||||
|
|
@ -513,8 +525,10 @@ module procedure mech_RGC_updateState
|
|||
real(pReal), dimension (3) :: nVect,surfCorr
|
||||
real(pReal), dimension (2) :: Gmoduli
|
||||
integer :: iGrain,iGNghb,iFace,i,j,k,l
|
||||
real(pReal) :: muGrain,muGNghb,nDefNorm,bgGrain,bgGNghb
|
||||
real(pReal), parameter :: nDefToler = 1.0e-10_pReal
|
||||
real(pReal) :: muGrain,muGNghb,nDefNorm
|
||||
real(pReal), parameter :: &
|
||||
nDefToler = 1.0e-10_pReal, &
|
||||
b = 2.5e-10_pReal ! Length of Burgers vector
|
||||
|
||||
nGDim = param(instance)%N_constituents
|
||||
rPen = 0.0_pReal
|
||||
|
|
@ -532,9 +546,7 @@ module procedure mech_RGC_updateState
|
|||
!-----------------------------------------------------------------------------------------------
|
||||
! computing the mismatch and penalty stress tensor of all grains
|
||||
grainLoop: do iGrain = 1,product(prm%N_constituents)
|
||||
Gmoduli = equivalentModuli(iGrain,ip,el)
|
||||
muGrain = Gmoduli(1) ! collecting the equivalent shear modulus of grain
|
||||
bgGrain = Gmoduli(2) ! and the lengthh of Burgers vector
|
||||
muGrain = equivalentMu(iGrain,ip,el)
|
||||
iGrain3 = grain1to3(iGrain,prm%N_constituents) ! get the grain ID in local 3-dimensional index (x,y,z)-position
|
||||
|
||||
interfaceLoop: do iFace = 1,6
|
||||
|
|
@ -546,9 +558,7 @@ module procedure mech_RGC_updateState
|
|||
where(iGNghb3 < 1) iGNghb3 = nGDim
|
||||
where(iGNghb3 >nGDim) iGNghb3 = 1
|
||||
iGNghb = grain3to1(iGNghb3,prm%N_constituents) ! get the ID of the neighboring grain
|
||||
Gmoduli = equivalentModuli(iGNghb,ip,el) ! collect the shear modulus and Burgers vector of the neighbor
|
||||
muGNghb = Gmoduli(1)
|
||||
bgGNghb = Gmoduli(2)
|
||||
muGNghb = equivalentMu(iGNghb,ip,el)
|
||||
gDef = 0.5_pReal*(fDef(1:3,1:3,iGNghb) - fDef(1:3,1:3,iGrain)) ! difference/jump in deformation gradeint across the neighbor
|
||||
|
||||
!-------------------------------------------------------------------------------------------
|
||||
|
|
@ -568,7 +578,7 @@ module procedure mech_RGC_updateState
|
|||
!-------------------------------------------------------------------------------------------
|
||||
! compute the stress penalty of all interfaces
|
||||
do i = 1,3; do j = 1,3; do k = 1,3; do l = 1,3
|
||||
rPen(i,j,iGrain) = rPen(i,j,iGrain) + 0.5_pReal*(muGrain*bgGrain + muGNghb*bgGNghb)*prm%xi_alpha &
|
||||
rPen(i,j,iGrain) = rPen(i,j,iGrain) + 0.5_pReal*(muGrain*b + muGNghb*b)*prm%xi_alpha &
|
||||
*surfCorr(abs(intFace(1)))/prm%D_alpha(abs(intFace(1))) &
|
||||
*cosh(prm%c_alpha*nDefNorm) &
|
||||
*0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_LeviCivita(k,l,j) &
|
||||
|
|
@ -655,44 +665,26 @@ module procedure mech_RGC_updateState
|
|||
end function surfaceCorrection
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!-------------------------------------------------------------------------------------------------
|
||||
!> @brief compute the equivalent shear and bulk moduli from the elasticity tensor
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
function equivalentModuli(grainID,ip,el)
|
||||
|
||||
real(pReal), dimension(2) :: equivalentModuli
|
||||
!-------------------------------------------------------------------------------------------------
|
||||
real(pReal) function equivalentMu(grainID,ip,el)
|
||||
|
||||
integer, intent(in) :: &
|
||||
grainID,&
|
||||
ip, & !< integration point number
|
||||
el !< element number
|
||||
real(pReal), dimension(6,6) :: elasTens
|
||||
real(pReal) :: &
|
||||
cEquiv_11, &
|
||||
cEquiv_12, &
|
||||
cEquiv_44
|
||||
|
||||
elasTens = constitutive_homogenizedC(grainID,ip,el)
|
||||
|
||||
!----------------------------------------------------------------------------------------------
|
||||
! compute the equivalent shear modulus after Turterltaub and Suiker, JMPS (2005)
|
||||
cEquiv_11 = (elasTens(1,1) + elasTens(2,2) + elasTens(3,3))/3.0_pReal
|
||||
cEquiv_12 = (elasTens(1,2) + elasTens(2,3) + elasTens(3,1) + &
|
||||
elasTens(1,3) + elasTens(2,1) + elasTens(3,2))/6.0_pReal
|
||||
cEquiv_44 = (elasTens(4,4) + elasTens(5,5) + elasTens(6,6))/3.0_pReal
|
||||
equivalentModuli(1) = 0.2_pReal*(cEquiv_11 - cEquiv_12) + 0.6_pReal*cEquiv_44
|
||||
|
||||
!----------------------------------------------------------------------------------------------
|
||||
! obtain the length of Burgers vector (could be model dependend)
|
||||
equivalentModuli(2) = 2.5e-10_pReal
|
||||
|
||||
end function equivalentModuli
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
equivalentMu = lattice_equivalent_mu(constitutive_homogenizedC(grainID,ip,el),'voigt')
|
||||
|
||||
end function equivalentMu
|
||||
|
||||
|
||||
!-------------------------------------------------------------------------------------------------
|
||||
!> @brief calculating the grain deformation gradient (the same with
|
||||
! homogenization_RGC_partitionDeformation, but used only for perturbation scheme)
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!-------------------------------------------------------------------------------------------------
|
||||
subroutine grainDeformation(F, avgF, instance, of)
|
||||
|
||||
real(pReal), dimension(:,:,:), intent(out) :: F !< partitioned F per grain
|
||||
|
|
@ -707,7 +699,7 @@ module procedure mech_RGC_updateState
|
|||
integer, dimension(3) :: iGrain3
|
||||
integer :: iGrain,iFace,i,j
|
||||
|
||||
!-------------------------------------------------------------------------------------------------
|
||||
!-----------------------------------------------------------------------------------------------
|
||||
! compute the deformation gradient of individual grains due to relaxations
|
||||
|
||||
associate(prm => param(instance))
|
||||
|
|
@ -729,7 +721,7 @@ module procedure mech_RGC_updateState
|
|||
|
||||
end subroutine grainDeformation
|
||||
|
||||
end procedure mech_RGC_updateState
|
||||
end function mech_RGC_updateState
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
|
|||
|
|
@ -421,6 +421,8 @@ module lattice
|
|||
lattice_BCT_ID, &
|
||||
lattice_HEX_ID, &
|
||||
lattice_ORT_ID, &
|
||||
lattice_equivalent_nu, &
|
||||
lattice_equivalent_mu, &
|
||||
lattice_applyLatticeSymmetry33, &
|
||||
lattice_SchmidMatrix_slip, &
|
||||
lattice_SchmidMatrix_twin, &
|
||||
|
|
@ -508,8 +510,8 @@ subroutine lattice_init
|
|||
|
||||
lattice_C66(1:6,1:6,p) = applyLatticeSymmetryC66(lattice_C66(1:6,1:6,p),phase%get_asString('lattice'))
|
||||
|
||||
lattice_mu(p) = equivalent_mu(lattice_C66(1:6,1:6,p),'voigt')
|
||||
lattice_nu(p) = equivalent_nu(lattice_C66(1:6,1:6,p),'voigt')
|
||||
lattice_nu(p) = lattice_equivalent_nu(lattice_C66(1:6,1:6,p),'voigt')
|
||||
lattice_mu(p) = lattice_equivalent_mu(lattice_C66(1:6,1:6,p),'voigt')
|
||||
|
||||
lattice_C66(1:6,1:6,p) = math_sym3333to66(math_Voigt66to3333(lattice_C66(1:6,1:6,p))) ! Literature data is in Voigt notation
|
||||
do i = 1, 6
|
||||
|
|
@ -2188,15 +2190,16 @@ end function getlabels
|
|||
!> @brief Equivalent Poisson's ratio (ν)
|
||||
!> @details https://doi.org/10.1143/JPSJ.20.635
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
function equivalent_nu(C,assumption) result(nu)
|
||||
function lattice_equivalent_nu(C,assumption) result(nu)
|
||||
|
||||
real(pReal), dimension(6,6), intent(in) :: C !< Stiffness tensor (Voigt notation)
|
||||
character(len=*), intent(in) :: assumption !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress)
|
||||
|
||||
real(pReal) :: K, mu, nu
|
||||
|
||||
logical :: error
|
||||
real(pReal), dimension(6,6) :: S
|
||||
|
||||
|
||||
if (IO_lc(assumption) == 'voigt') then
|
||||
K = (C(1,1)+C(2,2)+C(3,3) +2.0_pReal*(C(1,2)+C(2,3)+C(1,3))) &
|
||||
/ 9.0_pReal
|
||||
|
|
@ -2210,25 +2213,26 @@ function equivalent_nu(C,assumption) result(nu)
|
|||
K = 0.0_pReal
|
||||
endif
|
||||
|
||||
mu = equivalent_mu(C,assumption)
|
||||
mu = lattice_equivalent_mu(C,assumption)
|
||||
nu = (1.5_pReal*K -mu)/(3.0_pReal*K+mu)
|
||||
|
||||
end function equivalent_nu
|
||||
end function lattice_equivalent_nu
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief Equivalent shear modulus (μ)
|
||||
!> @details https://doi.org/10.1143/JPSJ.20.635
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
function equivalent_mu(C,assumption) result(mu)
|
||||
function lattice_equivalent_mu(C,assumption) result(mu)
|
||||
|
||||
real(pReal), dimension(6,6), intent(in) :: C !< Stiffness tensor (Voigt notation)
|
||||
character(len=*), intent(in) :: assumption !< Assumption ('Voigt' = isostrain, 'Reuss' = isostress)
|
||||
|
||||
real(pReal) :: mu
|
||||
|
||||
logical :: error
|
||||
real(pReal), dimension(6,6) :: S
|
||||
|
||||
|
||||
if (IO_lc(assumption) == 'voigt') then
|
||||
mu = (1.0_pReal*(C(1,1)+C(2,2)+C(3,3)) -1.0_pReal*(C(1,2)+C(2,3)+C(1,3)) +3.0_pReal*(C(4,4)+C(5,5)+C(6,6))) &
|
||||
/ 15.0_pReal
|
||||
|
|
@ -2242,7 +2246,7 @@ function equivalent_mu(C,assumption) result(mu)
|
|||
mu = 0.0_pReal
|
||||
endif
|
||||
|
||||
end function equivalent_mu
|
||||
end function lattice_equivalent_mu
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
@ -2266,14 +2270,14 @@ subroutine selfTest
|
|||
call random_number(C)
|
||||
C(1,1) = C(1,1) + 1.0_pReal
|
||||
C = applyLatticeSymmetryC66(C,'aP')
|
||||
if(dNeq(C(6,6),equivalent_mu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_mu/voigt'
|
||||
if(dNeq(C(6,6),equivalent_mu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_mu/reuss'
|
||||
if(dNeq(C(6,6),lattice_equivalent_mu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_mu/voigt'
|
||||
if(dNeq(C(6,6),lattice_equivalent_mu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_mu/reuss'
|
||||
|
||||
lambda = C(1,2)
|
||||
if(dNeq(lambda*0.5_pReal/(lambda+equivalent_mu(C,'voigt')),equivalent_nu(C,'voigt'),1.0e-12_pReal)) &
|
||||
error stop 'equivalent_nu/voigt'
|
||||
if(dNeq(lambda*0.5_pReal/(lambda+equivalent_mu(C,'reuss')),equivalent_nu(C,'reuss'),1.0e-12_pReal)) &
|
||||
error stop 'equivalent_nu/reuss'
|
||||
if(dNeq(lambda*0.5_pReal/(lambda+lattice_equivalent_mu(C,'voigt')), &
|
||||
lattice_equivalent_nu(C,'voigt'),1.0e-12_pReal)) error stop 'equivalent_nu/voigt'
|
||||
if(dNeq(lambda*0.5_pReal/(lambda+lattice_equivalent_mu(C,'reuss')), &
|
||||
lattice_equivalent_nu(C,'reuss'),1.0e-12_pReal)) error stop 'equivalent_nu/reuss'
|
||||
|
||||
end subroutine selfTest
|
||||
|
||||
|
|
|
|||
|
|
@ -279,9 +279,12 @@ real(pReal) pure function math_LeviCivita(i,j,k)
|
|||
|
||||
integer, intent(in) :: i,j,k
|
||||
|
||||
if (all([i,j,k] == [1,2,3]) .or. all([i,j,k] == [2,3,1]) .or. all([i,j,k] == [3,1,2])) then
|
||||
integer :: o
|
||||
|
||||
|
||||
if (any([(all(cshift([i,j,k],o) == [1,2,3]),o=0,2)])) then
|
||||
math_LeviCivita = +1.0_pReal
|
||||
elseif (all([i,j,k] == [3,2,1]) .or. all([i,j,k] == [2,1,3]) .or. all([i,j,k] == [1,3,2])) then
|
||||
elseif (any([(all(cshift([i,j,k],o) == [3,2,1]),o=0,2)])) then
|
||||
math_LeviCivita = -1.0_pReal
|
||||
else
|
||||
math_LeviCivita = 0.0_pReal
|
||||
|
|
|
|||
26
src/prec.f90
26
src/prec.f90
|
|
@ -108,8 +108,10 @@ logical elemental pure function dEq(a,b,tol)
|
|||
|
||||
real(pReal), intent(in) :: a,b
|
||||
real(pReal), intent(in), optional :: tol
|
||||
|
||||
real(pReal) :: eps
|
||||
|
||||
|
||||
if (present(tol)) then
|
||||
eps = tol
|
||||
else
|
||||
|
|
@ -132,11 +134,8 @@ logical elemental pure function dNeq(a,b,tol)
|
|||
real(pReal), intent(in) :: a,b
|
||||
real(pReal), intent(in), optional :: tol
|
||||
|
||||
if (present(tol)) then
|
||||
|
||||
dNeq = .not. dEq(a,b,tol)
|
||||
else
|
||||
dNeq = .not. dEq(a,b)
|
||||
endif
|
||||
|
||||
end function dNeq
|
||||
|
||||
|
|
@ -151,8 +150,10 @@ logical elemental pure function dEq0(a,tol)
|
|||
|
||||
real(pReal), intent(in) :: a
|
||||
real(pReal), intent(in), optional :: tol
|
||||
|
||||
real(pReal) :: eps
|
||||
|
||||
|
||||
if (present(tol)) then
|
||||
eps = tol
|
||||
else
|
||||
|
|
@ -175,11 +176,8 @@ logical elemental pure function dNeq0(a,tol)
|
|||
real(pReal), intent(in) :: a
|
||||
real(pReal), intent(in), optional :: tol
|
||||
|
||||
if (present(tol)) then
|
||||
|
||||
dNeq0 = .not. dEq0(a,tol)
|
||||
else
|
||||
dNeq0 = .not. dEq0(a)
|
||||
endif
|
||||
|
||||
end function dNeq0
|
||||
|
||||
|
|
@ -195,8 +193,10 @@ logical elemental pure function cEq(a,b,tol)
|
|||
|
||||
complex(pReal), intent(in) :: a,b
|
||||
real(pReal), intent(in), optional :: tol
|
||||
|
||||
real(pReal) :: eps
|
||||
|
||||
|
||||
if (present(tol)) then
|
||||
eps = tol
|
||||
else
|
||||
|
|
@ -220,11 +220,8 @@ logical elemental pure function cNeq(a,b,tol)
|
|||
complex(pReal), intent(in) :: a,b
|
||||
real(pReal), intent(in), optional :: tol
|
||||
|
||||
if (present(tol)) then
|
||||
|
||||
cNeq = .not. cEq(a,b,tol)
|
||||
else
|
||||
cNeq = .not. cEq(a,b)
|
||||
endif
|
||||
|
||||
end function cNeq
|
||||
|
||||
|
|
@ -238,6 +235,7 @@ pure function prec_bytesToC_FLOAT(bytes)
|
|||
real(C_FLOAT), dimension(size(bytes,kind=pI64)/(storage_size(0._C_FLOAT,pI64)/8_pI64)) :: &
|
||||
prec_bytesToC_FLOAT
|
||||
|
||||
|
||||
prec_bytesToC_FLOAT = transfer(bytes,prec_bytesToC_FLOAT,size(prec_bytesToC_FLOAT))
|
||||
|
||||
end function prec_bytesToC_FLOAT
|
||||
|
|
@ -252,6 +250,7 @@ pure function prec_bytesToC_DOUBLE(bytes)
|
|||
real(C_DOUBLE), dimension(size(bytes,kind=pI64)/(storage_size(0._C_DOUBLE,pI64)/8_pI64)) :: &
|
||||
prec_bytesToC_DOUBLE
|
||||
|
||||
|
||||
prec_bytesToC_DOUBLE = transfer(bytes,prec_bytesToC_DOUBLE,size(prec_bytesToC_DOUBLE))
|
||||
|
||||
end function prec_bytesToC_DOUBLE
|
||||
|
|
@ -266,6 +265,7 @@ pure function prec_bytesToC_INT32_T(bytes)
|
|||
integer(C_INT32_T), dimension(size(bytes,kind=pI64)/(storage_size(0_C_INT32_T,pI64)/8_pI64)) :: &
|
||||
prec_bytesToC_INT32_T
|
||||
|
||||
|
||||
prec_bytesToC_INT32_T = transfer(bytes,prec_bytesToC_INT32_T,size(prec_bytesToC_INT32_T))
|
||||
|
||||
end function prec_bytesToC_INT32_T
|
||||
|
|
@ -280,6 +280,7 @@ pure function prec_bytesToC_INT64_T(bytes)
|
|||
integer(C_INT64_T), dimension(size(bytes,kind=pI64)/(storage_size(0_C_INT64_T,pI64)/8_pI64)) :: &
|
||||
prec_bytesToC_INT64_T
|
||||
|
||||
|
||||
prec_bytesToC_INT64_T = transfer(bytes,prec_bytesToC_INT64_T,size(prec_bytesToC_INT64_T))
|
||||
|
||||
end function prec_bytesToC_INT64_T
|
||||
|
|
@ -295,6 +296,7 @@ subroutine selfTest
|
|||
integer(pInt), dimension(1) :: i
|
||||
real(pReal), dimension(2) :: r
|
||||
|
||||
|
||||
realloc_lhs_test = [1,2]
|
||||
if (any(realloc_lhs_test/=[1,2])) error stop 'LHS allocation'
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue