Merge branch 'development' into new-gmsh-version
This commit is contained in:
Sharan Roongta
commit
1611b74c79
|
|
@ -74,9 +74,23 @@ end subroutine CPFEM_initAll
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine CPFEM_init
|
||||
|
||||
integer(HID_T) :: fileHandle
|
||||
character(len=pStringLen) :: fileName
|
||||
|
||||
|
||||
print'(/,a)', ' <<<+- CPFEM init -+>>>'; flush(IO_STDOUT)
|
||||
|
||||
if (interface_restartInc > 0) call crystallite_restartRead
|
||||
|
||||
if (interface_restartInc > 0) then
|
||||
print'(/,a,i0,a)', ' reading restart information of increment from file'; flush(IO_STDOUT)
|
||||
write(fileName,'(a,i0,a)') trim(getSolverJobName())//'_',worldrank,'.hdf5'
|
||||
fileHandle = HDF5_openFile(fileName)
|
||||
|
||||
call homogenization_restartRead(fileHandle)
|
||||
call constitutive_restartRead(fileHandle)
|
||||
|
||||
call HDF5_closeFile(fileHandle)
|
||||
endif
|
||||
|
||||
end subroutine CPFEM_init
|
||||
|
||||
|
|
@ -86,7 +100,19 @@ end subroutine CPFEM_init
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine CPFEM_restartWrite
|
||||
|
||||
call crystallite_restartWrite
|
||||
integer(HID_T) :: fileHandle
|
||||
character(len=pStringLen) :: fileName
|
||||
|
||||
|
||||
print*, ' writing field and constitutive data required for restart to file';flush(IO_STDOUT)
|
||||
|
||||
write(fileName,'(a,i0,a)') trim(getSolverJobName())//'_',worldrank,'.hdf5'
|
||||
fileHandle = HDF5_openFile(fileName,'a')
|
||||
|
||||
call homogenization_restartWrite(fileHandle)
|
||||
call constitutive_restartWrite(fileHandle)
|
||||
|
||||
call HDF5_closeFile(fileHandle)
|
||||
|
||||
end subroutine CPFEM_restartWrite
|
||||
|
||||
|
|
|
|||
|
|
@ -15,7 +15,6 @@ module constitutive
|
|||
use discretization
|
||||
use parallelization
|
||||
use HDF5_utilities
|
||||
use DAMASK_interface
|
||||
use results
|
||||
|
||||
implicit none
|
||||
|
|
@ -42,42 +41,14 @@ module constitutive
|
|||
KINEMATICS_SLIPPLANE_OPENING_ID, &
|
||||
KINEMATICS_THERMAL_EXPANSION_ID
|
||||
end enum
|
||||
real(pReal), dimension(:,:,:), allocatable :: &
|
||||
crystallite_subdt !< substepped time increment of each grain
|
||||
|
||||
type(rotation), dimension(:,:,:), allocatable :: &
|
||||
crystallite_orientation !< current orientation
|
||||
real(pReal), dimension(:,:,:,:,:), allocatable :: &
|
||||
crystallite_F0, & !< def grad at start of FE inc
|
||||
crystallite_Fe, & !< current "elastic" def grad (end of converged time step)
|
||||
crystallite_subFp0,& !< plastic def grad at start of crystallite inc
|
||||
crystallite_subFi0,& !< intermediate def grad at start of crystallite inc
|
||||
crystallite_Lp0, & !< plastic velocitiy grad at start of FE inc
|
||||
crystallite_partitionedLp0, & !< plastic velocity grad at start of homog inc
|
||||
crystallite_S0, & !< 2nd Piola-Kirchhoff stress vector at start of FE inc
|
||||
crystallite_partitionedS0 !< 2nd Piola-Kirchhoff stress vector at start of homog inc
|
||||
real(pReal), dimension(:,:,:,:,:), allocatable, public :: &
|
||||
crystallite_P, & !< 1st Piola-Kirchhoff stress per grain
|
||||
crystallite_Lp, & !< current plastic velocitiy grad (end of converged time step)
|
||||
crystallite_S, & !< current 2nd Piola-Kirchhoff stress vector (end of converged time step)
|
||||
crystallite_partitionedF0, & !< def grad at start of homog inc
|
||||
crystallite_F !< def grad to be reached at end of homog inc
|
||||
|
||||
type :: tTensorContainer
|
||||
real(pReal), dimension(:,:,:), allocatable :: data
|
||||
end type
|
||||
|
||||
type(tTensorContainer), dimension(:), allocatable :: &
|
||||
constitutive_mech_Fi, &
|
||||
constitutive_mech_Fi0, &
|
||||
constitutive_mech_partitionedFi0, &
|
||||
constitutive_mech_Li, &
|
||||
constitutive_mech_Li0, &
|
||||
constitutive_mech_partitionedLi0, &
|
||||
constitutive_mech_Fp, &
|
||||
constitutive_mech_Fp0, &
|
||||
constitutive_mech_partitionedFp0
|
||||
|
||||
|
||||
type :: tNumerics
|
||||
integer :: &
|
||||
iJacoLpresiduum, & !< frequency of Jacobian update of residuum in Lp
|
||||
|
|
@ -162,10 +133,6 @@ module constitutive
|
|||
end subroutine damage_results
|
||||
|
||||
|
||||
module subroutine mech_restart_read(fileHandle)
|
||||
integer(HID_T), intent(in) :: fileHandle
|
||||
end subroutine mech_restart_read
|
||||
|
||||
module subroutine mech_initializeRestorationPoints(ph,me)
|
||||
integer, intent(in) :: ph, me
|
||||
end subroutine mech_initializeRestorationPoints
|
||||
|
|
@ -185,6 +152,61 @@ module constitutive
|
|||
includeL
|
||||
end subroutine mech_restore
|
||||
|
||||
module function constitutive_mech_dPdF(dt,co,ip,el) result(dPdF)
|
||||
real(pReal), intent(in) :: dt
|
||||
integer, intent(in) :: &
|
||||
co, & !< counter in constituent loop
|
||||
ip, & !< counter in integration point loop
|
||||
el !< counter in element loop
|
||||
real(pReal), dimension(3,3,3,3) :: dPdF
|
||||
end function constitutive_mech_dPdF
|
||||
|
||||
module subroutine mech_restartWrite(groupHandle,ph)
|
||||
integer(HID_T), intent(in) :: groupHandle
|
||||
integer, intent(in) :: ph
|
||||
end subroutine mech_restartWrite
|
||||
|
||||
module subroutine mech_restartRead(groupHandle,ph)
|
||||
integer(HID_T), intent(in) :: groupHandle
|
||||
integer, intent(in) :: ph
|
||||
end subroutine mech_restartRead
|
||||
|
||||
|
||||
module function constitutive_mech_getS(co,ip,el) result(S)
|
||||
integer, intent(in) :: co, ip, el
|
||||
real(pReal), dimension(3,3) :: S
|
||||
end function constitutive_mech_getS
|
||||
|
||||
module function constitutive_mech_getLp(co,ip,el) result(Lp)
|
||||
integer, intent(in) :: co, ip, el
|
||||
real(pReal), dimension(3,3) :: Lp
|
||||
end function constitutive_mech_getLp
|
||||
|
||||
module function constitutive_mech_getF(co,ip,el) result(F)
|
||||
integer, intent(in) :: co, ip, el
|
||||
real(pReal), dimension(3,3) :: F
|
||||
end function constitutive_mech_getF
|
||||
|
||||
module function constitutive_mech_getF_e(co,ip,el) result(F_e)
|
||||
integer, intent(in) :: co, ip, el
|
||||
real(pReal), dimension(3,3) :: F_e
|
||||
end function constitutive_mech_getF_e
|
||||
|
||||
module function constitutive_mech_getP(co,ip,el) result(P)
|
||||
integer, intent(in) :: co, ip, el
|
||||
real(pReal), dimension(3,3) :: P
|
||||
end function constitutive_mech_getP
|
||||
|
||||
module function constitutive_thermal_T(co,ip,el) result(T)
|
||||
integer, intent(in) :: co, ip, el
|
||||
real(pReal) :: T
|
||||
end function constitutive_thermal_T
|
||||
|
||||
module subroutine constitutive_mech_setF(F,co,ip,el)
|
||||
real(pReal), dimension(3,3), intent(in) :: F
|
||||
integer, intent(in) :: co, ip, el
|
||||
end subroutine constitutive_mech_setF
|
||||
|
||||
! == cleaned:end ===================================================================================
|
||||
|
||||
module function crystallite_stress(dt,co,ip,el) result(converged_)
|
||||
|
|
@ -238,15 +260,12 @@ module constitutive
|
|||
dPhiDot_dPhi
|
||||
end subroutine constitutive_damage_getRateAndItsTangents
|
||||
|
||||
module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, S, Lp, ip, el)
|
||||
module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, ip, el)
|
||||
integer, intent(in) :: &
|
||||
ip, & !< integration point number
|
||||
el !< element number
|
||||
real(pReal), intent(in) :: &
|
||||
T
|
||||
real(pReal), intent(in), dimension(:,:,:,:,:) :: &
|
||||
S, & !< current 2nd Piola Kitchoff stress vector
|
||||
Lp !< plastic velocity gradient
|
||||
real(pReal), intent(inout) :: &
|
||||
TDot, &
|
||||
dTDot_dT
|
||||
|
|
@ -342,13 +361,11 @@ module constitutive
|
|||
end subroutine constitutive_plastic_LpAndItsTangents
|
||||
|
||||
|
||||
module subroutine constitutive_plastic_dependentState(F, co, ip, el)
|
||||
module subroutine constitutive_plastic_dependentState(co,ip,el)
|
||||
integer, intent(in) :: &
|
||||
co, & !< component-ID of integration point
|
||||
ip, & !< integration point
|
||||
el !< element
|
||||
real(pReal), intent(in), dimension(3,3) :: &
|
||||
F !< elastic deformation gradient
|
||||
end subroutine constitutive_plastic_dependentState
|
||||
|
||||
|
||||
|
|
@ -390,12 +407,17 @@ module constitutive
|
|||
converged, &
|
||||
crystallite_init, &
|
||||
crystallite_stress, &
|
||||
crystallite_stressTangent, &
|
||||
constitutive_mech_dPdF, &
|
||||
crystallite_orientations, &
|
||||
crystallite_push33ToRef, &
|
||||
crystallite_restartWrite, &
|
||||
constitutive_restartWrite, &
|
||||
constitutive_restartRead, &
|
||||
integrateSourceState, &
|
||||
crystallite_restartRead, &
|
||||
constitutive_mech_setF, &
|
||||
constitutive_mech_getP, &
|
||||
constitutive_mech_getLp, &
|
||||
constitutive_mech_getF, &
|
||||
constitutive_mech_getS, &
|
||||
constitutive_initializeRestorationPoints, &
|
||||
constitutive_windForward, &
|
||||
PLASTICITY_UNDEFINED_ID, &
|
||||
|
|
@ -615,7 +637,7 @@ end subroutine constitutive_LiAndItsTangents
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief contains the constitutive equation for calculating the rate of change of microstructure
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
function constitutive_damage_collectDotState(S, co, ip, el,ph,of) result(broken)
|
||||
function constitutive_damage_collectDotState(co,ip,el,ph,of) result(broken)
|
||||
|
||||
integer, intent(in) :: &
|
||||
co, & !< component-ID of integration point
|
||||
|
|
@ -623,8 +645,6 @@ function constitutive_damage_collectDotState(S, co, ip, el,ph,of) result(broken)
|
|||
el, & !< element
|
||||
ph, &
|
||||
of
|
||||
real(pReal), intent(in), dimension(3,3) :: &
|
||||
S !< 2nd Piola Kirchhoff stress (vector notation)
|
||||
integer :: &
|
||||
so !< counter in source loop
|
||||
logical :: broken
|
||||
|
|
@ -637,7 +657,7 @@ function constitutive_damage_collectDotState(S, co, ip, el,ph,of) result(broken)
|
|||
sourceType: select case (phase_source(so,ph))
|
||||
|
||||
case (SOURCE_damage_anisoBrittle_ID) sourceType
|
||||
call source_damage_anisoBrittle_dotState(S, co, ip, el) ! correct stress?
|
||||
call source_damage_anisoBrittle_dotState(constitutive_mech_getS(co,ip,el), co, ip, el) ! correct stress?
|
||||
|
||||
case (SOURCE_damage_isoDuctile_ID) sourceType
|
||||
call source_damage_isoDuctile_dotState(co, ip, el)
|
||||
|
|
@ -748,7 +768,6 @@ subroutine constitutive_allocateState(state, &
|
|||
allocate(state%atol (sizeState), source=0.0_pReal)
|
||||
allocate(state%state0 (sizeState,Nconstituents), source=0.0_pReal)
|
||||
allocate(state%partitionedState0(sizeState,Nconstituents), source=0.0_pReal)
|
||||
allocate(state%subState0 (sizeState,Nconstituents), source=0.0_pReal)
|
||||
allocate(state%state (sizeState,Nconstituents), source=0.0_pReal)
|
||||
|
||||
allocate(state%dotState (sizeDotState,Nconstituents), source=0.0_pReal)
|
||||
|
|
@ -792,17 +811,14 @@ end subroutine constitutive_restore
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine constitutive_forward
|
||||
|
||||
integer :: i, j
|
||||
integer :: ph, so
|
||||
|
||||
crystallite_F0 = crystallite_F
|
||||
crystallite_Lp0 = crystallite_Lp
|
||||
crystallite_S0 = crystallite_S
|
||||
|
||||
call constitutive_mech_forward()
|
||||
|
||||
do i = 1, size(sourceState)
|
||||
do j = 1,phase_Nsources(i)
|
||||
sourceState(i)%p(j)%state0 = sourceState(i)%p(j)%state
|
||||
do ph = 1, size(sourceState)
|
||||
do so = 1,phase_Nsources(ph)
|
||||
sourceState(ph)%p(so)%state0 = sourceState(ph)%p(so)%state
|
||||
enddo; enddo
|
||||
|
||||
end subroutine constitutive_forward
|
||||
|
|
@ -838,12 +854,12 @@ end subroutine constitutive_results
|
|||
subroutine crystallite_init
|
||||
|
||||
integer :: &
|
||||
Nconstituents, &
|
||||
ph, &
|
||||
me, &
|
||||
co, & !< counter in integration point component loop
|
||||
ip, & !< counter in integration point loop
|
||||
el, & !< counter in element loop
|
||||
so, &
|
||||
cMax, & !< maximum number of integration point components
|
||||
iMax, & !< maximum number of integration points
|
||||
eMax !< maximum number of elements
|
||||
|
|
@ -866,19 +882,6 @@ subroutine crystallite_init
|
|||
iMax = discretization_nIPs
|
||||
eMax = discretization_Nelems
|
||||
|
||||
allocate(crystallite_F(3,3,cMax,iMax,eMax),source=0.0_pReal)
|
||||
|
||||
allocate(crystallite_S0, &
|
||||
crystallite_F0,crystallite_Lp0, &
|
||||
crystallite_partitionedS0, &
|
||||
crystallite_partitionedF0,&
|
||||
crystallite_partitionedLp0, &
|
||||
crystallite_S,crystallite_P, &
|
||||
crystallite_Fe,crystallite_Lp, &
|
||||
crystallite_subFp0,crystallite_subFi0, &
|
||||
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)
|
||||
|
|
@ -914,27 +917,10 @@ subroutine crystallite_init
|
|||
|
||||
phases => config_material%get('phase')
|
||||
|
||||
allocate(constitutive_mech_Fi(phases%length))
|
||||
allocate(constitutive_mech_Fi0(phases%length))
|
||||
allocate(constitutive_mech_partitionedFi0(phases%length))
|
||||
allocate(constitutive_mech_Fp(phases%length))
|
||||
allocate(constitutive_mech_Fp0(phases%length))
|
||||
allocate(constitutive_mech_partitionedFp0(phases%length))
|
||||
allocate(constitutive_mech_Li(phases%length))
|
||||
allocate(constitutive_mech_Li0(phases%length))
|
||||
allocate(constitutive_mech_partitionedLi0(phases%length))
|
||||
do ph = 1, phases%length
|
||||
Nconstituents = count(material_phaseAt == ph) * discretization_nIPs
|
||||
|
||||
allocate(constitutive_mech_Fi(ph)%data(3,3,Nconstituents))
|
||||
allocate(constitutive_mech_Fi0(ph)%data(3,3,Nconstituents))
|
||||
allocate(constitutive_mech_partitionedFi0(ph)%data(3,3,Nconstituents))
|
||||
allocate(constitutive_mech_Fp(ph)%data(3,3,Nconstituents))
|
||||
allocate(constitutive_mech_Fp0(ph)%data(3,3,Nconstituents))
|
||||
allocate(constitutive_mech_partitionedFp0(ph)%data(3,3,Nconstituents))
|
||||
allocate(constitutive_mech_Li(ph)%data(3,3,Nconstituents))
|
||||
allocate(constitutive_mech_Li0(ph)%data(3,3,Nconstituents))
|
||||
allocate(constitutive_mech_partitionedLi0(ph)%data(3,3,Nconstituents))
|
||||
do so = 1, phase_Nsources(ph)
|
||||
allocate(sourceState(ph)%p(so)%subState0,source=sourceState(ph)%p(so)%state0) ! ToDo: hack
|
||||
enddo
|
||||
enddo
|
||||
|
||||
print'(a42,1x,i10)', ' # of elements: ', eMax
|
||||
|
|
@ -942,34 +928,6 @@ subroutine crystallite_init
|
|||
print'(a42,1x,i10)', 'max # of constituents/integration point: ', cMax
|
||||
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))
|
||||
|
||||
ph = material_phaseAt(co,el)
|
||||
me = material_phaseMemberAt(co,ip,el)
|
||||
constitutive_mech_Fp0(ph)%data(1:3,1:3,me) = material_orientation0(co,ip,el)%asMatrix() ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005)
|
||||
constitutive_mech_Fp0(ph)%data(1:3,1:3,me) = constitutive_mech_Fp0(ph)%data(1:3,1:3,me) &
|
||||
/ math_det33(constitutive_mech_Fp0(ph)%data(1:3,1:3,me))**(1.0_pReal/3.0_pReal)
|
||||
constitutive_mech_Fi0(ph)%data(1:3,1:3,me) = math_I3
|
||||
|
||||
crystallite_F0(1:3,1:3,co,ip,el) = math_I3
|
||||
|
||||
crystallite_Fe(1:3,1:3,co,ip,el) = math_inv33(matmul(constitutive_mech_Fi0(ph)%data(1:3,1:3,me), &
|
||||
constitutive_mech_Fp0(ph)%data(1:3,1:3,me))) ! assuming that euler angles are given in internal strain free configuration
|
||||
constitutive_mech_Fp(ph)%data(1:3,1:3,me) = constitutive_mech_Fp0(ph)%data(1:3,1:3,me)
|
||||
constitutive_mech_Fi(ph)%data(1:3,1:3,me) = constitutive_mech_Fi0(ph)%data(1:3,1:3,me)
|
||||
|
||||
constitutive_mech_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
|
||||
!$OMP END PARALLEL DO
|
||||
|
||||
crystallite_partitionedF0 = crystallite_F0
|
||||
crystallite_F = crystallite_F0
|
||||
|
||||
|
||||
!$OMP PARALLEL DO PRIVATE(ph,me)
|
||||
do el = 1, size(material_phaseMemberAt,3)
|
||||
|
|
@ -978,7 +936,7 @@ subroutine crystallite_init
|
|||
ph = material_phaseAt(co,el)
|
||||
me = material_phaseMemberAt(co,ip,el)
|
||||
call crystallite_orientations(co,ip,el)
|
||||
call constitutive_plastic_dependentState(crystallite_partitionedF0(1:3,1:3,co,ip,el),co,ip,el) ! update dependent state variables to be consistent with basic states
|
||||
call constitutive_plastic_dependentState(co,ip,el) ! update dependent state variables to be consistent with basic states
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
|
@ -1003,15 +961,11 @@ subroutine constitutive_initializeRestorationPoints(ip,el)
|
|||
do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
|
||||
ph = material_phaseAt(co,el)
|
||||
me = material_phaseMemberAt(co,ip,el)
|
||||
crystallite_partitionedLp0(1:3,1:3,co,ip,el) = crystallite_Lp0(1:3,1:3,co,ip,el)
|
||||
crystallite_partitionedF0(1:3,1:3,co,ip,el) = crystallite_F0(1:3,1:3,co,ip,el)
|
||||
crystallite_partitionedS0(1:3,1:3,co,ip,el) = crystallite_S0(1:3,1:3,co,ip,el)
|
||||
|
||||
call mech_initializeRestorationPoints(ph,me)
|
||||
|
||||
do 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))
|
||||
sourceState(ph)%p(so)%partitionedState0(:,me) = sourceState(ph)%p(so)%state0(:,me)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
|
|
@ -1035,9 +989,6 @@ subroutine constitutive_windForward(ip,el)
|
|||
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_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 so = 1, phase_Nsources(material_phaseAt(co,el))
|
||||
|
|
@ -1048,136 +999,6 @@ subroutine constitutive_windForward(ip,el)
|
|||
end subroutine constitutive_windForward
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief Calculate tangent (dPdF).
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
function crystallite_stressTangent(co,ip,el) result(dPdF)
|
||||
|
||||
real(pReal), dimension(3,3,3,3) :: dPdF
|
||||
integer, intent(in) :: &
|
||||
co, & !< counter in constituent loop
|
||||
ip, & !< counter in integration point loop
|
||||
el !< counter in element loop
|
||||
|
||||
integer :: &
|
||||
o, &
|
||||
p, ph, me
|
||||
real(pReal), dimension(3,3) :: devNull, &
|
||||
invSubFp0,invSubFi0,invFp,invFi, &
|
||||
temp_33_1, temp_33_2, temp_33_3
|
||||
real(pReal), dimension(3,3,3,3) :: dSdFe, &
|
||||
dSdF, &
|
||||
dSdFi, &
|
||||
dLidS, & ! tangent in lattice configuration
|
||||
dLidFi, &
|
||||
dLpdS, &
|
||||
dLpdFi, &
|
||||
dFidS, &
|
||||
dFpinvdF, &
|
||||
rhs_3333, &
|
||||
lhs_3333, &
|
||||
temp_3333
|
||||
real(pReal), dimension(9,9):: temp_99
|
||||
logical :: error
|
||||
|
||||
|
||||
ph = material_phaseAt(co,el)
|
||||
me = material_phaseMemberAt(co,ip,el)
|
||||
|
||||
call constitutive_hooke_SandItsTangents(devNull,dSdFe,dSdFi, &
|
||||
crystallite_Fe(1:3,1:3,co,ip,el), &
|
||||
constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el)
|
||||
call constitutive_LiAndItsTangents(devNull,dLidS,dLidFi, &
|
||||
crystallite_S (1:3,1:3,co,ip,el), &
|
||||
constitutive_mech_Fi(ph)%data(1:3,1:3,me), &
|
||||
co,ip,el)
|
||||
|
||||
invFp = math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,me))
|
||||
invFi = math_inv33(constitutive_mech_Fi(ph)%data(1:3,1:3,me))
|
||||
invSubFp0 = math_inv33(crystallite_subFp0(1:3,1:3,co,ip,el))
|
||||
invSubFi0 = math_inv33(crystallite_subFi0(1:3,1:3,co,ip,el))
|
||||
|
||||
if (sum(abs(dLidS)) < tol_math_check) then
|
||||
dFidS = 0.0_pReal
|
||||
else
|
||||
lhs_3333 = 0.0_pReal; rhs_3333 = 0.0_pReal
|
||||
do o=1,3; do p=1,3
|
||||
lhs_3333(1:3,1:3,o,p) = lhs_3333(1:3,1:3,o,p) &
|
||||
+ crystallite_subdt(co,ip,el)*matmul(invSubFi0,dLidFi(1:3,1:3,o,p))
|
||||
lhs_3333(1:3,o,1:3,p) = lhs_3333(1:3,o,1:3,p) &
|
||||
+ invFi*invFi(p,o)
|
||||
rhs_3333(1:3,1:3,o,p) = rhs_3333(1:3,1:3,o,p) &
|
||||
- crystallite_subdt(co,ip,el)*matmul(invSubFi0,dLidS(1:3,1:3,o,p))
|
||||
enddo; enddo
|
||||
call math_invert(temp_99,error,math_3333to99(lhs_3333))
|
||||
if (error) then
|
||||
call IO_warning(warning_ID=600,el=el,ip=ip,g=co, &
|
||||
ext_msg='inversion error in analytic tangent calculation')
|
||||
dFidS = 0.0_pReal
|
||||
else
|
||||
dFidS = math_mul3333xx3333(math_99to3333(temp_99),rhs_3333)
|
||||
endif
|
||||
dLidS = math_mul3333xx3333(dLidFi,dFidS) + dLidS
|
||||
endif
|
||||
|
||||
call constitutive_plastic_LpAndItsTangents(devNull,dLpdS,dLpdFi, &
|
||||
crystallite_S (1:3,1:3,co,ip,el), &
|
||||
constitutive_mech_Fi(ph)%data(1:3,1:3,me),co,ip,el)
|
||||
dLpdS = math_mul3333xx3333(dLpdFi,dFidS) + dLpdS
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! calculate dSdF
|
||||
temp_33_1 = transpose(matmul(invFp,invFi))
|
||||
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)
|
||||
temp_3333(1:3,1:3,p,o) = matmul(matmul(temp_33_2,dLpdS(1:3,1:3,p,o)), invFi) &
|
||||
+ matmul(temp_33_3,dLidS(1:3,1:3,p,o))
|
||||
enddo; enddo
|
||||
lhs_3333 = crystallite_subdt(co,ip,el)*math_mul3333xx3333(dSdFe,temp_3333) &
|
||||
+ math_mul3333xx3333(dSdFi,dFidS)
|
||||
|
||||
call math_invert(temp_99,error,math_eye(9)+math_3333to99(lhs_3333))
|
||||
if (error) then
|
||||
call IO_warning(warning_ID=600,el=el,ip=ip,g=co, &
|
||||
ext_msg='inversion error in analytic tangent calculation')
|
||||
dSdF = rhs_3333
|
||||
else
|
||||
dSdF = math_mul3333xx3333(math_99to3333(temp_99),rhs_3333)
|
||||
endif
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! calculate dFpinvdF
|
||||
temp_3333 = math_mul3333xx3333(dLpdS,dSdF)
|
||||
do o=1,3; do p=1,3
|
||||
dFpinvdF(1:3,1:3,p,o) = -crystallite_subdt(co,ip,el) &
|
||||
* matmul(invSubFp0, matmul(temp_3333(1:3,1:3,p,o),invFi))
|
||||
enddo; enddo
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! assemble dPdF
|
||||
temp_33_1 = matmul(crystallite_S(1:3,1:3,co,ip,el),transpose(invFp))
|
||||
temp_33_2 = matmul(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(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_F(1:3,1:3,co,ip,el), &
|
||||
dFpinvdF(1:3,1:3,p,o)),temp_33_1) &
|
||||
+ matmul(matmul(temp_33_2,dSdF(1:3,1:3,p,o)), &
|
||||
transpose(invFp)) &
|
||||
+ matmul(temp_33_3,transpose(dFpinvdF(1:3,1:3,p,o)))
|
||||
enddo; enddo
|
||||
|
||||
end function crystallite_stressTangent
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief calculates orientations
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
|
|
@ -1189,7 +1010,8 @@ subroutine crystallite_orientations(co,ip,el)
|
|||
el !< counter in element loop
|
||||
|
||||
|
||||
call crystallite_orientation(co,ip,el)%fromMatrix(transpose(math_rotationalPart(crystallite_Fe(1:3,1:3,co,ip,el))))
|
||||
call crystallite_orientation(co,ip,el)%fromMatrix(transpose(math_rotationalPart(&
|
||||
constitutive_mech_getF_e(co,ip,el))))
|
||||
|
||||
if (plasticState(material_phaseAt(1,el))%nonlocal) &
|
||||
call plastic_nonlocal_updateCompatibility(crystallite_orientation, &
|
||||
|
|
@ -1205,17 +1027,17 @@ end subroutine crystallite_orientations
|
|||
function crystallite_push33ToRef(co,ip,el, tensor33)
|
||||
|
||||
real(pReal), dimension(3,3), intent(in) :: tensor33
|
||||
real(pReal), dimension(3,3) :: T
|
||||
integer, intent(in):: &
|
||||
el, &
|
||||
ip, &
|
||||
co
|
||||
|
||||
real(pReal), dimension(3,3) :: crystallite_push33ToRef
|
||||
|
||||
real(pReal), dimension(3,3) :: T
|
||||
|
||||
|
||||
T = matmul(material_orientation0(co,ip,el)%asMatrix(),transpose(math_inv33(constitutive_mech_getF(co,ip,el)))) ! ToDo: initial orientation correct?
|
||||
|
||||
T = matmul(material_orientation0(co,ip,el)%asMatrix(), & ! ToDo: initial orientation correct?
|
||||
transpose(math_inv33(crystallite_F(1:3,1:3,co,ip,el))))
|
||||
crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T))
|
||||
|
||||
end function crystallite_push33ToRef
|
||||
|
|
@ -1254,7 +1076,7 @@ function integrateSourceState(dt,co,ip,el) result(broken)
|
|||
|
||||
converged_ = .true.
|
||||
broken = constitutive_thermal_collectDotState(ph,me)
|
||||
broken = broken .or. constitutive_damage_collectDotState(crystallite_S(1:3,1:3,co,ip,el), co,ip,el,ph,me)
|
||||
broken = broken .or. constitutive_damage_collectDotState(co,ip,el,ph,me)
|
||||
if(broken) return
|
||||
|
||||
do so = 1, phase_Nsources(ph)
|
||||
|
|
@ -1272,7 +1094,7 @@ function integrateSourceState(dt,co,ip,el) result(broken)
|
|||
enddo
|
||||
|
||||
broken = constitutive_thermal_collectDotState(ph,me)
|
||||
broken = broken .or. constitutive_damage_collectDotState(crystallite_S(1:3,1:3,co,ip,el), co,ip,el,ph,me)
|
||||
broken = broken .or. constitutive_damage_collectDotState(co,ip,el,ph,me)
|
||||
if(broken) exit iteration
|
||||
|
||||
do so = 1, phase_Nsources(ph)
|
||||
|
|
@ -1292,7 +1114,7 @@ function integrateSourceState(dt,co,ip,el) result(broken)
|
|||
enddo
|
||||
|
||||
if(converged_) then
|
||||
broken = constitutive_damage_deltaState(crystallite_Fe(1:3,1:3,co,ip,el),co,ip,el,ph,me)
|
||||
broken = constitutive_damage_deltaState(constitutive_mech_getF_e(co,ip,el),co,ip,el,ph,me)
|
||||
exit iteration
|
||||
endif
|
||||
|
||||
|
|
@ -1345,90 +1167,60 @@ end function converged
|
|||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief Write current restart information (Field and constitutive data) to file.
|
||||
! ToDo: Merge data into one file for MPI, move state to constitutive and homogenization, respectively
|
||||
! ToDo: Merge data into one file for MPI
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine crystallite_restartWrite
|
||||
subroutine constitutive_restartWrite(fileHandle)
|
||||
|
||||
integer(HID_T), intent(in) :: fileHandle
|
||||
|
||||
integer(HID_T), dimension(2) :: groupHandle
|
||||
integer :: ph
|
||||
integer(HID_T) :: fileHandle, groupHandle
|
||||
character(len=pStringLen) :: fileName, datasetName
|
||||
|
||||
print*, ' writing field and constitutive data required for restart to file';flush(IO_STDOUT)
|
||||
|
||||
write(fileName,'(a,i0,a)') trim(getSolverJobName())//'_',worldrank,'.hdf5'
|
||||
fileHandle = HDF5_openFile(fileName,'a')
|
||||
groupHandle(1) = HDF5_addGroup(fileHandle,'phase')
|
||||
|
||||
call HDF5_write(fileHandle,crystallite_F,'F')
|
||||
call HDF5_write(fileHandle,crystallite_Lp, 'L_p')
|
||||
call HDF5_write(fileHandle,crystallite_S, 'S')
|
||||
do ph = 1, size(material_name_phase)
|
||||
|
||||
groupHandle(2) = HDF5_addGroup(groupHandle(1),material_name_phase(ph))
|
||||
|
||||
call mech_restartWrite(groupHandle(2),ph)
|
||||
|
||||
call HDF5_closeGroup(groupHandle(2))
|
||||
|
||||
groupHandle = HDF5_addGroup(fileHandle,'phase')
|
||||
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 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)
|
||||
call HDF5_closeGroup(groupHandle(1))
|
||||
|
||||
call HDF5_closeFile(fileHandle)
|
||||
|
||||
end subroutine crystallite_restartWrite
|
||||
end subroutine constitutive_restartWrite
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief Read data for restart
|
||||
! ToDo: Merge data into one file for MPI, move state to constitutive and homogenization, respectively
|
||||
! ToDo: Merge data into one file for MPI
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine crystallite_restartRead
|
||||
subroutine constitutive_restartRead(fileHandle)
|
||||
|
||||
integer(HID_T), intent(in) :: fileHandle
|
||||
|
||||
integer(HID_T), dimension(2) :: groupHandle
|
||||
integer :: ph
|
||||
integer(HID_T) :: fileHandle, groupHandle
|
||||
character(len=pStringLen) :: fileName, datasetName
|
||||
|
||||
print'(/,a,i0,a)', ' reading restart information of increment from file'
|
||||
|
||||
write(fileName,'(a,i0,a)') trim(getSolverJobName())//'_',worldrank,'.hdf5'
|
||||
fileHandle = HDF5_openFile(fileName)
|
||||
groupHandle(1) = HDF5_openGroup(fileHandle,'phase')
|
||||
|
||||
call HDF5_read(fileHandle,crystallite_F0, 'F')
|
||||
call HDF5_read(fileHandle,crystallite_Lp0,'L_p')
|
||||
call HDF5_read(fileHandle,crystallite_S0, 'S')
|
||||
do ph = 1, size(material_name_phase)
|
||||
|
||||
groupHandle(2) = HDF5_openGroup(groupHandle(1),material_name_phase(ph))
|
||||
|
||||
call mech_restartRead(groupHandle(2),ph)
|
||||
|
||||
call HDF5_closeGroup(groupHandle(2))
|
||||
|
||||
groupHandle = HDF5_openGroup(fileHandle,'phase')
|
||||
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 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)
|
||||
call HDF5_closeGroup(groupHandle(1))
|
||||
|
||||
call HDF5_closeFile(fileHandle)
|
||||
|
||||
end subroutine crystallite_restartRead
|
||||
end subroutine constitutive_restartRead
|
||||
|
||||
|
||||
end module constitutive
|
||||
|
|
|
|||
File diff suppressed because it is too large
Load Diff
|
|
@ -552,10 +552,8 @@ end function plastic_nonlocal_init
|
|||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief calculates quantities characterizing the microstructure
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
module subroutine plastic_nonlocal_dependentState(F, instance, of, ip, el)
|
||||
module subroutine plastic_nonlocal_dependentState(instance, of, ip, el)
|
||||
|
||||
real(pReal), dimension(3,3), intent(in) :: &
|
||||
F
|
||||
integer, intent(in) :: &
|
||||
instance, &
|
||||
of, &
|
||||
|
|
@ -647,7 +645,7 @@ module subroutine plastic_nonlocal_dependentState(F, instance, of, ip, el)
|
|||
ph = material_phaseAt(1,el)
|
||||
me = material_phaseMemberAt(1,ip,el)
|
||||
invFp = math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,me))
|
||||
invFe = matmul(constitutive_mech_Fp(ph)%data(1:3,1:3,me),math_inv33(F))
|
||||
invFe = math_inv33(constitutive_mech_Fe(ph)%data(1:3,1:3,me))
|
||||
|
||||
rho_edg_delta = rho0(:,mob_edg_pos) - rho0(:,mob_edg_neg)
|
||||
rho_scr_delta = rho0(:,mob_scr_pos) - rho0(:,mob_scr_neg)
|
||||
|
|
@ -976,13 +974,11 @@ end subroutine plastic_nonlocal_deltaState
|
|||
!---------------------------------------------------------------------------------------------------
|
||||
!> @brief calculates the rate of change of microstructure
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
module subroutine plastic_nonlocal_dotState(Mp, F, Temperature,timestep, &
|
||||
module subroutine plastic_nonlocal_dotState(Mp, Temperature,timestep, &
|
||||
instance,of,ip,el)
|
||||
|
||||
real(pReal), dimension(3,3), intent(in) :: &
|
||||
Mp !< MandelStress
|
||||
real(pReal), dimension(3,3,homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems), intent(in) :: &
|
||||
F !< Deformation gradient
|
||||
real(pReal), intent(in) :: &
|
||||
Temperature, & !< temperature
|
||||
timestep !< substepped crystallite time increment
|
||||
|
|
@ -1149,7 +1145,7 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Temperature,timestep, &
|
|||
- rhoDip(s,1) / timestep - rhoDotAthermalAnnihilation(s,9) &
|
||||
- rhoDotSingle2DipoleGlide(s,9)) ! make sure that we do not annihilate more dipoles than we have
|
||||
|
||||
rhoDot = rhoDotFlux(F,timestep, instance,of,ip,el) &
|
||||
rhoDot = rhoDotFlux(timestep, instance,of,ip,el) &
|
||||
+ rhoDotMultiplication &
|
||||
+ rhoDotSingle2DipoleGlide &
|
||||
+ rhoDotAthermalAnnihilation &
|
||||
|
|
@ -1178,10 +1174,8 @@ end subroutine plastic_nonlocal_dotState
|
|||
!---------------------------------------------------------------------------------------------------
|
||||
!> @brief calculates the rate of change of microstructure
|
||||
!---------------------------------------------------------------------------------------------------
|
||||
function rhoDotFlux(F,timestep, instance,of,ip,el)
|
||||
function rhoDotFlux(timestep,instance,of,ip,el)
|
||||
|
||||
real(pReal), dimension(3,3,homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems), intent(in) :: &
|
||||
F !< Deformation gradient
|
||||
real(pReal), intent(in) :: &
|
||||
timestep !< substepped crystallite time increment
|
||||
integer, intent(in) :: &
|
||||
|
|
@ -1293,7 +1287,7 @@ function rhoDotFlux(F,timestep, instance,of,ip,el)
|
|||
m(1:3,:,3) = -prm%slip_transverse
|
||||
m(1:3,:,4) = prm%slip_transverse
|
||||
|
||||
my_F = F(1:3,1:3,1,ip,el)
|
||||
my_F = constitutive_mech_F(ph)%data(1:3,1:3,of)
|
||||
my_Fe = matmul(my_F, math_inv33(constitutive_mech_Fp(ph)%data(1:3,1:3,of)))
|
||||
|
||||
neighbors: do n = 1,nIPneighbors
|
||||
|
|
@ -1311,7 +1305,7 @@ function rhoDotFlux(F,timestep, instance,of,ip,el)
|
|||
|
||||
if (neighbor_n > 0) then ! if neighbor exists, average deformation gradient
|
||||
neighbor_instance = phase_plasticityInstance(material_phaseAt(1,neighbor_el))
|
||||
neighbor_F = F(1:3,1:3,1,neighbor_ip,neighbor_el)
|
||||
neighbor_F = constitutive_mech_F(np)%data(1:3,1:3,no)
|
||||
neighbor_Fe = matmul(neighbor_F, math_inv33(constitutive_mech_Fp(np)%data(1:3,1:3,no)))
|
||||
Favg = 0.5_pReal * (my_F + neighbor_F)
|
||||
else ! if no neighbor, take my value as average
|
||||
|
|
|
|||
|
|
@ -68,15 +68,13 @@ end subroutine thermal_init
|
|||
!----------------------------------------------------------------------------------------------
|
||||
!< @brief calculates thermal dissipation rate
|
||||
!----------------------------------------------------------------------------------------------
|
||||
module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, S, Lp, ip, el)
|
||||
module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, ip, el)
|
||||
|
||||
integer, intent(in) :: &
|
||||
ip, & !< integration point number
|
||||
el !< element number
|
||||
real(pReal), intent(in) :: &
|
||||
T
|
||||
real(pReal), intent(in), dimension(:,:,:,:,:) :: &
|
||||
S, & !< current 2nd Piola Kirchhoff stress
|
||||
Lp !< plastic velocity gradient
|
||||
T !< plastic velocity gradient
|
||||
real(pReal), intent(inout) :: &
|
||||
TDot, &
|
||||
dTDot_dT
|
||||
|
|
@ -84,6 +82,7 @@ module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T,
|
|||
real(pReal) :: &
|
||||
my_Tdot, &
|
||||
my_dTdot_dT
|
||||
real(pReal), dimension(3,3) :: Lp, S
|
||||
integer :: &
|
||||
phase, &
|
||||
homog, &
|
||||
|
|
@ -101,10 +100,10 @@ module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T,
|
|||
do source = 1, phase_Nsources(phase)
|
||||
select case(phase_source(source,phase))
|
||||
case (SOURCE_thermal_dissipation_ID)
|
||||
Lp = constitutive_mech_getLp(grain,ip,el)
|
||||
S = constitutive_mech_getS(grain,ip,el)
|
||||
call source_thermal_dissipation_getRateAndItsTangent(my_Tdot, my_dTdot_dT, &
|
||||
S(1:3,1:3,grain,ip,el), &
|
||||
Lp(1:3,1:3,grain,ip,el), &
|
||||
phase)
|
||||
S, Lp, phase)
|
||||
|
||||
case (SOURCE_thermal_externalheat_ID)
|
||||
call source_thermal_externalheat_getRateAndItsTangent(my_Tdot, my_dTdot_dT, &
|
||||
|
|
@ -122,4 +121,22 @@ module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T,
|
|||
end subroutine constitutive_thermal_getRateAndItsTangents
|
||||
|
||||
|
||||
|
||||
|
||||
! getter for non-thermal (e.g. mech)
|
||||
module function constitutive_thermal_T(co,ip,el) result(T)
|
||||
|
||||
integer, intent(in) :: co, ip, el
|
||||
real(pReal) :: T
|
||||
|
||||
integer :: ho, tme
|
||||
|
||||
ho = material_homogenizationAt(el)
|
||||
tme = material_homogenizationMemberAt(ip,el)
|
||||
|
||||
T = temperature(ho)%p(tme)
|
||||
|
||||
end function constitutive_thermal_T
|
||||
|
||||
|
||||
end submodule constitutive_thermal
|
||||
|
|
|
|||
|
|
@ -16,6 +16,7 @@ module homogenization
|
|||
use thermal_conduction
|
||||
use damage_none
|
||||
use damage_nonlocal
|
||||
use HDF5_utilities
|
||||
use results
|
||||
|
||||
implicit none
|
||||
|
|
@ -63,7 +64,8 @@ module homogenization
|
|||
el !< element number
|
||||
end subroutine mech_partition
|
||||
|
||||
module subroutine mech_homogenize(ip,el)
|
||||
module subroutine mech_homogenize(dt,ip,el)
|
||||
real(pReal), intent(in) :: dt
|
||||
integer, intent(in) :: &
|
||||
ip, & !< integration point
|
||||
el !< element number
|
||||
|
|
@ -91,7 +93,9 @@ module homogenization
|
|||
homogenization_init, &
|
||||
materialpoint_stressAndItsTangent, &
|
||||
homogenization_forward, &
|
||||
homogenization_results
|
||||
homogenization_results, &
|
||||
homogenization_restartRead, &
|
||||
homogenization_restartWrite
|
||||
|
||||
contains
|
||||
|
||||
|
|
@ -168,10 +172,8 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
|
|||
converged = .false. ! pretend failed step ...
|
||||
subStep = 1.0_pReal/num%subStepSizeHomog ! ... larger then the requested calculation
|
||||
|
||||
if (homogState(ho)%sizeState > 0) &
|
||||
homogState(ho)%subState0(:,me) = homogState(ho)%State0(:,me)
|
||||
if (damageState(ho)%sizeState > 0) &
|
||||
damageState(ho)%subState0(:,me) = damageState(ho)%State0(:,me)
|
||||
if (homogState(ho)%sizeState > 0) homogState(ho)%subState0(:,me) = homogState(ho)%State0(:,me)
|
||||
if (damageState(ho)%sizeState > 0) damageState(ho)%subState0(:,me) = damageState(ho)%State0(:,me)
|
||||
|
||||
cutBackLooping: do while (.not. terminallyIll .and. subStep > num%subStepMinHomog)
|
||||
|
||||
|
|
@ -184,10 +186,8 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
|
|||
! wind forward grain starting point
|
||||
call constitutive_windForward(ip,el)
|
||||
|
||||
if(homogState(ho)%sizeState > 0) &
|
||||
homogState(ho)%subState0(:,me) = homogState(ho)%State(:,me)
|
||||
if(damageState(ho)%sizeState > 0) &
|
||||
damageState(ho)%subState0(:,me) = damageState(ho)%State(:,me)
|
||||
if(homogState(ho)%sizeState > 0) homogState(ho)%subState0(:,me) = homogState(ho)%State(:,me)
|
||||
if(damageState(ho)%sizeState > 0) damageState(ho)%subState0(:,me) = damageState(ho)%State(:,me)
|
||||
|
||||
endif steppingNeeded
|
||||
elseif ( (myNgrains == 1 .and. subStep <= 1.0 ) .or. & ! single grain already tried internal subStepping in crystallite
|
||||
|
|
@ -201,10 +201,8 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
|
|||
|
||||
call constitutive_restore(ip,el,subStep < 1.0_pReal)
|
||||
|
||||
if(homogState(ho)%sizeState > 0) &
|
||||
homogState(ho)%State(:,me) = homogState(ho)%subState0(:,me)
|
||||
if(damageState(ho)%sizeState > 0) &
|
||||
damageState(ho)%State(:,me) = damageState(ho)%subState0(:,me)
|
||||
if(homogState(ho)%sizeState > 0) homogState(ho)%State(:,me) = homogState(ho)%subState0(:,me)
|
||||
if(damageState(ho)%sizeState > 0) damageState(ho)%State(:,me) = damageState(ho)%subState0(:,me)
|
||||
endif
|
||||
|
||||
if (subStep > num%subStepMinHomog) doneAndHappy = [.false.,.true.]
|
||||
|
|
@ -257,7 +255,7 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE
|
|||
do co = 1, myNgrains
|
||||
call crystallite_orientations(co,ip,el)
|
||||
enddo
|
||||
call mech_homogenize(ip,el)
|
||||
call mech_homogenize(dt,ip,el)
|
||||
enddo IpLooping3
|
||||
enddo elementLooping3
|
||||
!$OMP END PARALLEL DO
|
||||
|
|
@ -320,4 +318,59 @@ subroutine homogenization_forward
|
|||
|
||||
end subroutine homogenization_forward
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine homogenization_restartWrite(fileHandle)
|
||||
|
||||
integer(HID_T), intent(in) :: fileHandle
|
||||
|
||||
integer(HID_T), dimension(2) :: groupHandle
|
||||
integer :: ho
|
||||
|
||||
|
||||
groupHandle(1) = HDF5_addGroup(fileHandle,'homogenization')
|
||||
|
||||
do ho = 1, size(material_name_homogenization)
|
||||
|
||||
groupHandle(2) = HDF5_addGroup(groupHandle(1),material_name_homogenization(ho))
|
||||
|
||||
call HDF5_read(groupHandle(2),homogState(ho)%state,'omega') ! ToDo: should be done by mech
|
||||
|
||||
call HDF5_closeGroup(groupHandle(2))
|
||||
|
||||
enddo
|
||||
|
||||
call HDF5_closeGroup(groupHandle(1))
|
||||
|
||||
end subroutine homogenization_restartWrite
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
subroutine homogenization_restartRead(fileHandle)
|
||||
|
||||
integer(HID_T), intent(in) :: fileHandle
|
||||
|
||||
integer(HID_T), dimension(2) :: groupHandle
|
||||
integer :: ho
|
||||
|
||||
|
||||
groupHandle(1) = HDF5_openGroup(fileHandle,'homogenization')
|
||||
|
||||
do ho = 1, size(material_name_homogenization)
|
||||
|
||||
groupHandle(2) = HDF5_openGroup(groupHandle(1),material_name_homogenization(ho))
|
||||
|
||||
call HDF5_write(groupHandle(2),homogState(ho)%state,'omega') ! ToDo: should be done by mech
|
||||
|
||||
call HDF5_closeGroup(groupHandle(2))
|
||||
|
||||
enddo
|
||||
|
||||
call HDF5_closeGroup(groupHandle(1))
|
||||
|
||||
end subroutine homogenization_restartRead
|
||||
|
||||
|
||||
end module homogenization
|
||||
|
|
|
|||
|
|
@ -52,12 +52,11 @@ submodule(homogenization) homogenization_mech
|
|||
end subroutine mech_RGC_averageStressAndItsTangent
|
||||
|
||||
|
||||
module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) result(doneAndHappy)
|
||||
module function mech_RGC_updateState(P,F,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
|
||||
F !< partitioned 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
|
||||
|
|
@ -113,67 +112,73 @@ module subroutine mech_partition(subF,ip,el)
|
|||
ip, & !< integration point
|
||||
el !< element number
|
||||
|
||||
integer :: co
|
||||
real(pReal) :: F(3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
|
||||
|
||||
|
||||
chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
|
||||
|
||||
case (HOMOGENIZATION_NONE_ID) chosenHomogenization
|
||||
crystallite_F(1:3,1:3,1,ip,el) = subF
|
||||
F(1:3,1:3,1) = subF
|
||||
|
||||
case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
|
||||
call mech_isostrain_partitionDeformation(&
|
||||
crystallite_F(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
|
||||
subF)
|
||||
call mech_isostrain_partitionDeformation(F,subF)
|
||||
|
||||
case (HOMOGENIZATION_RGC_ID) chosenHomogenization
|
||||
call mech_RGC_partitionDeformation(&
|
||||
crystallite_F(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
|
||||
subF,&
|
||||
ip, &
|
||||
el)
|
||||
call mech_RGC_partitionDeformation(F,subF,ip,el)
|
||||
|
||||
end select chosenHomogenization
|
||||
|
||||
do co = 1,homogenization_Nconstituents(material_homogenizationAt(el))
|
||||
call constitutive_mech_setF(F(1:3,1:3,co),co,ip,el)
|
||||
enddo
|
||||
|
||||
|
||||
end subroutine mech_partition
|
||||
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
!> @brief Average P and dPdF from the individual constituents.
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
module subroutine mech_homogenize(ip,el)
|
||||
module subroutine mech_homogenize(dt,ip,el)
|
||||
|
||||
real(pReal), intent(in) :: dt
|
||||
integer, intent(in) :: &
|
||||
ip, & !< integration point
|
||||
el !< element number
|
||||
|
||||
integer :: co,ce
|
||||
real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
|
||||
real(pReal) :: Ps(3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
|
||||
|
||||
|
||||
ce = (el-1)* discretization_nIPs + ip
|
||||
chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
|
||||
|
||||
case (HOMOGENIZATION_NONE_ID) chosenHomogenization
|
||||
homogenization_P(1:3,1:3,ce) = crystallite_P(1:3,1:3,1,ip,el)
|
||||
homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = crystallite_stressTangent(1,ip,el)
|
||||
homogenization_P(1:3,1:3,ce) = constitutive_mech_getP(1,ip,el)
|
||||
homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = constitutive_mech_dPdF(dt,1,ip,el)
|
||||
|
||||
case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
|
||||
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
|
||||
dPdFs(:,:,:,:,co) = crystallite_stressTangent(co,ip,el)
|
||||
dPdFs(:,:,:,:,co) = constitutive_mech_dPdF(dt,co,ip,el)
|
||||
Ps(:,:,co) = constitutive_mech_getP(co,ip,el)
|
||||
enddo
|
||||
call mech_isostrain_averageStressAndItsTangent(&
|
||||
homogenization_P(1:3,1:3,ce), &
|
||||
homogenization_dPdF(1:3,1:3,1:3,1:3,ce),&
|
||||
crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
|
||||
dPdFs, &
|
||||
Ps,dPdFs, &
|
||||
homogenization_typeInstance(material_homogenizationAt(el)))
|
||||
|
||||
case (HOMOGENIZATION_RGC_ID) chosenHomogenization
|
||||
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
|
||||
dPdFs(:,:,:,:,co) = crystallite_stressTangent(co,ip,el)
|
||||
dPdFs(:,:,:,:,co) = constitutive_mech_dPdF(dt,co,ip,el)
|
||||
Ps(:,:,co) = constitutive_mech_getP(co,ip,el)
|
||||
enddo
|
||||
call mech_RGC_averageStressAndItsTangent(&
|
||||
homogenization_P(1:3,1:3,ce), &
|
||||
homogenization_dPdF(1:3,1:3,1:3,1:3,ce),&
|
||||
crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
|
||||
dPdFs, &
|
||||
Ps,dPdFs, &
|
||||
homogenization_typeInstance(material_homogenizationAt(el)))
|
||||
|
||||
end select chosenHomogenization
|
||||
|
|
@ -198,21 +203,17 @@ module function mech_updateState(subdt,subF,ip,el) result(doneAndHappy)
|
|||
|
||||
integer :: co
|
||||
real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
|
||||
real(pReal) :: Fs(3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
|
||||
real(pReal) :: Ps(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)
|
||||
dPdFs(:,:,:,:,co) = constitutive_mech_dPdF(subdt,co,ip,el)
|
||||
Fs(:,:,co) = constitutive_mech_getF(co,ip,el)
|
||||
Ps(:,:,co) = constitutive_mech_getP(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)
|
||||
doneAndHappy = mech_RGC_updateState(Ps,Fs,subF,subdt,dPdFs,ip,el)
|
||||
else
|
||||
doneAndHappy = .true.
|
||||
endif
|
||||
|
|
|
|||
|
|
@ -24,9 +24,6 @@ submodule(homogenization:homogenization_mech) homogenization_mech_RGC
|
|||
end type tParameters
|
||||
|
||||
type :: tRGCstate
|
||||
real(pReal), pointer, dimension(:) :: &
|
||||
work, &
|
||||
penaltyEnergy
|
||||
real(pReal), pointer, dimension(:,:) :: &
|
||||
relaxationVector
|
||||
end type tRGCstate
|
||||
|
|
@ -170,8 +167,7 @@ module subroutine mech_RGC_init(num_homogMech)
|
|||
nIntFaceTot = 3*( (prm%N_constituents(1)-1)*prm%N_constituents(2)*prm%N_constituents(3) &
|
||||
+ prm%N_constituents(1)*(prm%N_constituents(2)-1)*prm%N_constituents(3) &
|
||||
+ prm%N_constituents(1)*prm%N_constituents(2)*(prm%N_constituents(3)-1))
|
||||
sizeState = nIntFaceTot &
|
||||
+ size(['avg constitutive work ','average penalty energy'])
|
||||
sizeState = nIntFaceTot
|
||||
|
||||
homogState(h)%sizeState = sizeState
|
||||
allocate(homogState(h)%state0 (sizeState,Nmaterialpoints), source=0.0_pReal)
|
||||
|
|
@ -180,8 +176,6 @@ module subroutine mech_RGC_init(num_homogMech)
|
|||
|
||||
stt%relaxationVector => homogState(h)%state(1:nIntFaceTot,:)
|
||||
st0%relaxationVector => homogState(h)%state0(1:nIntFaceTot,:)
|
||||
stt%work => homogState(h)%state(nIntFaceTot+1,:)
|
||||
stt%penaltyEnergy => homogState(h)%state(nIntFaceTot+2,:)
|
||||
|
||||
allocate(dst%volumeDiscrepancy( Nmaterialpoints), source=0.0_pReal)
|
||||
allocate(dst%relaxationRate_avg( Nmaterialpoints), source=0.0_pReal)
|
||||
|
|
@ -243,12 +237,11 @@ end subroutine mech_RGC_partitionDeformation
|
|||
!> @brief update the internal state of the homogenization scheme and tell whether "done" and
|
||||
! "happy" with result
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) result(doneAndHappy)
|
||||
module function mech_RGC_updateState(P,F,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
|
||||
F !< partitioned 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
|
||||
|
|
@ -346,17 +339,6 @@ module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) result(doneAndHa
|
|||
if (residMax < num%rtol*stresMax .or. residMax < num%atol) then
|
||||
doneAndHappy = .true.
|
||||
|
||||
!--------------------------------------------------------------------------------------------------
|
||||
! compute/update the state for postResult, i.e., all energy densities computed by time-integration
|
||||
do iGrain = 1,product(prm%N_constituents)
|
||||
do i = 1,3;do j = 1,3
|
||||
stt%work(of) = stt%work(of) &
|
||||
+ P(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain))/real(nGrain,pReal)
|
||||
stt%penaltyEnergy(of) = stt%penaltyEnergy(of) &
|
||||
+ R(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain))/real(nGrain,pReal)
|
||||
enddo; enddo
|
||||
enddo
|
||||
|
||||
dst%mismatch(1:3,of) = sum(NN,2)/real(nGrain,pReal)
|
||||
dst%relaxationRate_avg(of) = sum(abs(drelax))/dt/real(3*nIntFaceTot,pReal)
|
||||
dst%relaxationRate_max(of) = maxval(abs(drelax))/dt
|
||||
|
|
@ -523,7 +505,6 @@ module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el) result(doneAndHa
|
|||
integer, dimension (3) :: iGrain3,iGNghb3,nGDim
|
||||
real(pReal), dimension (3,3) :: gDef,nDef
|
||||
real(pReal), dimension (3) :: nVect,surfCorr
|
||||
real(pReal), dimension (2) :: Gmoduli
|
||||
integer :: iGrain,iGNghb,iFace,i,j,k,l
|
||||
real(pReal) :: muGrain,muGNghb,nDefNorm
|
||||
real(pReal), parameter :: &
|
||||
|
|
@ -755,15 +736,9 @@ module subroutine mech_RGC_results(instance,group)
|
|||
associate(stt => state(instance), dst => dependentState(instance), prm => param(instance))
|
||||
outputsLoop: do o = 1,size(prm%output)
|
||||
select case(trim(prm%output(o)))
|
||||
case('W')
|
||||
call results_writeDataset(group,stt%work,trim(prm%output(o)), &
|
||||
'work density','J/m³')
|
||||
case('M')
|
||||
call results_writeDataset(group,dst%mismatch,trim(prm%output(o)), &
|
||||
'average mismatch tensor','1')
|
||||
case('R')
|
||||
call results_writeDataset(group,stt%penaltyEnergy,trim(prm%output(o)), &
|
||||
'mismatch penalty density','J/m³')
|
||||
case('Delta_V')
|
||||
call results_writeDataset(group,dst%volumeDiscrepancy,trim(prm%output(o)), &
|
||||
'volume discrepancy','m³')
|
||||
|
|
|
|||
|
|
@ -94,7 +94,7 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
|
|||
dTdot_dT = 0.0_pReal
|
||||
|
||||
homog = material_homogenizationAt(el)
|
||||
call constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, crystallite_S,crystallite_Lp ,ip, el)
|
||||
call constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T, ip, el)
|
||||
|
||||
Tdot = Tdot/real(homogenization_Nconstituents(homog),pReal)
|
||||
dTdot_dT = dTdot_dT/real(homogenization_Nconstituents(homog),pReal)
|
||||
|
|
@ -112,14 +112,16 @@ function thermal_conduction_getConductivity(ip,el)
|
|||
el !< element number
|
||||
real(pReal), dimension(3,3) :: &
|
||||
thermal_conduction_getConductivity
|
||||
|
||||
integer :: &
|
||||
grain
|
||||
co
|
||||
|
||||
|
||||
thermal_conduction_getConductivity = 0.0_pReal
|
||||
do grain = 1, homogenization_Nconstituents(material_homogenizationAt(el))
|
||||
|
||||
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
|
||||
thermal_conduction_getConductivity = thermal_conduction_getConductivity + &
|
||||
crystallite_push33ToRef(grain,ip,el,lattice_K(:,:,material_phaseAt(grain,el)))
|
||||
crystallite_push33ToRef(co,ip,el,lattice_K(:,:,material_phaseAt(co,el)))
|
||||
enddo
|
||||
|
||||
thermal_conduction_getConductivity = thermal_conduction_getConductivity &
|
||||
|
|
@ -138,14 +140,16 @@ function thermal_conduction_getSpecificHeat(ip,el)
|
|||
el !< element number
|
||||
real(pReal) :: &
|
||||
thermal_conduction_getSpecificHeat
|
||||
|
||||
integer :: &
|
||||
grain
|
||||
co
|
||||
|
||||
|
||||
thermal_conduction_getSpecificHeat = 0.0_pReal
|
||||
|
||||
do grain = 1, homogenization_Nconstituents(material_homogenizationAt(el))
|
||||
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
|
||||
thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat &
|
||||
+ lattice_c_p(material_phaseAt(grain,el))
|
||||
+ lattice_c_p(material_phaseAt(co,el))
|
||||
enddo
|
||||
|
||||
thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat &
|
||||
|
|
@ -164,15 +168,16 @@ function thermal_conduction_getMassDensity(ip,el)
|
|||
el !< element number
|
||||
real(pReal) :: &
|
||||
thermal_conduction_getMassDensity
|
||||
|
||||
integer :: &
|
||||
grain
|
||||
co
|
||||
|
||||
|
||||
thermal_conduction_getMassDensity = 0.0_pReal
|
||||
|
||||
|
||||
do grain = 1, homogenization_Nconstituents(material_homogenizationAt(el))
|
||||
do co = 1, homogenization_Nconstituents(material_homogenizationAt(el))
|
||||
thermal_conduction_getMassDensity = thermal_conduction_getMassDensity &
|
||||
+ lattice_rho(material_phaseAt(grain,el))
|
||||
+ lattice_rho(material_phaseAt(co,el))
|
||||
enddo
|
||||
|
||||
thermal_conduction_getMassDensity = thermal_conduction_getMassDensity &
|
||||
|
|
|
|||
Loading…
Reference in New Issue