223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

Merge branch 'clean-material.yaml-parse' into development

This commit is contained in:
Sharan Roongta 2020-10-29 12:50:40 +01:00
parent 76a57af101 568d1a020c
commit 4b25097a9a
42 changed files with 554 additions and 539 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
src/CPFEM.f90
View File

@ -107,9 +107,9 @@ subroutine CPFEM_init
print'(/,a)', ' <<<+- CPFEM init -+>>>'; flush(IO_STDOUT) print'(/,a)', ' <<<+- CPFEM init -+>>>'; flush(IO_STDOUT)
allocate(CPFEM_cs( 6,discretization_nIP,discretization_nElem), source= 0.0_pReal) allocate(CPFEM_cs( 6,discretization_nIPs,discretization_Nelems), source= 0.0_pReal)
allocate(CPFEM_dcsdE( 6,6,discretization_nIP,discretization_nElem), source= 0.0_pReal) allocate(CPFEM_dcsdE( 6,6,discretization_nIPs,discretization_Nelems), source= 0.0_pReal)
allocate(CPFEM_dcsdE_knownGood(6,6,discretization_nIP,discretization_nElem), source= 0.0_pReal) allocate(CPFEM_dcsdE_knownGood(6,6,discretization_nIPs,discretization_Nelems), source= 0.0_pReal)
!------------------------------------------------------------------------------ !------------------------------------------------------------------------------
! read debug options ! read debug options
@ -184,8 +184,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
temperature(material_homogenizationAt(elCP))%p(thermalMapping(material_homogenizationAt(elCP))%p(ip,elCP)) = & temperature(material_homogenizationAt(elCP))%p(thermalMapping(material_homogenizationAt(elCP))%p(ip,elCP)) = &
temperature_inp temperature_inp
end select chosenThermal1 end select chosenThermal1
materialpoint_F0(1:3,1:3,ip,elCP) = ffn homogenization_F0(1:3,1:3,ip,elCP) = ffn
materialpoint_F(1:3,1:3,ip,elCP) = ffn1 homogenization_F(1:3,1:3,ip,elCP) = ffn1
if (iand(mode, CPFEM_CALCRESULTS) /= 0_pInt) then if (iand(mode, CPFEM_CALCRESULTS) /= 0_pInt) then
@ -212,17 +212,17 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
else terminalIllness else terminalIllness
! translate from P to sigma ! translate from P to sigma
Kirchhoff = matmul(materialpoint_P(1:3,1:3,ip,elCP), transpose(materialpoint_F(1:3,1:3,ip,elCP))) Kirchhoff = matmul(homogenization_P(1:3,1:3,ip,elCP), transpose(homogenization_F(1:3,1:3,ip,elCP)))
J_inverse = 1.0_pReal / math_det33(materialpoint_F(1:3,1:3,ip,elCP)) J_inverse = 1.0_pReal / math_det33(homogenization_F(1:3,1:3,ip,elCP))
CPFEM_cs(1:6,ip,elCP) = math_sym33to6(J_inverse * Kirchhoff,weighted=.false.) CPFEM_cs(1:6,ip,elCP) = math_sym33to6(J_inverse * Kirchhoff,weighted=.false.)
! translate from dP/dF to dCS/dE ! translate from dP/dF to dCS/dE
H = 0.0_pReal H = 0.0_pReal
do i=1,3; do j=1,3; do k=1,3; do l=1,3; do m=1,3; do n=1,3 do i=1,3; do j=1,3; do k=1,3; do l=1,3; do m=1,3; do n=1,3
H(i,j,k,l) = H(i,j,k,l) & H(i,j,k,l) = H(i,j,k,l) &
+ materialpoint_F(j,m,ip,elCP) * materialpoint_F(l,n,ip,elCP) & + homogenization_F(j,m,ip,elCP) * homogenization_F(l,n,ip,elCP) &
* materialpoint_dPdF(i,m,k,n,ip,elCP) & * homogenization_dPdF(i,m,k,n,ip,elCP) &
- math_delta(j,l) * materialpoint_F(i,m,ip,elCP) * materialpoint_P(k,m,ip,elCP) & - math_delta(j,l) * homogenization_F(i,m,ip,elCP) * homogenization_P(k,m,ip,elCP) &
+ 0.5_pReal * ( Kirchhoff(j,l)*math_delta(i,k) + Kirchhoff(i,k)*math_delta(j,l) & + 0.5_pReal * ( Kirchhoff(j,l)*math_delta(i,k) + Kirchhoff(i,k)*math_delta(j,l) &
+ Kirchhoff(j,k)*math_delta(i,l) + Kirchhoff(i,l)*math_delta(j,k)) + Kirchhoff(j,k)*math_delta(i,l) + Kirchhoff(i,l)*math_delta(j,k))
enddo; enddo; enddo; enddo; enddo; enddo enddo; enddo; enddo; enddo; enddo; enddo

16
src/YAML_types.f90
View File

@ -583,9 +583,9 @@ function tNode_get_byIndex_asStrings(self,i) result(nodeAsStrings)
end function tNode_get_byIndex_asStrings end function tNode_get_byIndex_asStrings
!------------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Returns the key in a dictionary as a string !> @brief Returns the key in a dictionary as a string
!------------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function tNode_getKey_byIndex(self,i) result(key) function tNode_getKey_byIndex(self,i) result(key)
class(tNode), intent(in), target :: self class(tNode), intent(in), target :: self
@ -1169,17 +1169,21 @@ function tList_asStrings(self)
type(tScalar), pointer :: scalar type(tScalar), pointer :: scalar
len_max = 0 len_max = 0
allocate(character(len=pStringLen) :: tList_asStrings(self%length)) item => self%first
do i = 1, self%length
scalar => item%node%asScalar()
len_max = max(len_max, len_trim(scalar%asString()))
item => item%next
enddo
allocate(character(len=len_max) :: tList_asStrings(self%length))
item => self%first item => self%first
do i = 1, self%length do i = 1, self%length
scalar => item%node%asScalar() scalar => item%node%asScalar()
tList_asStrings(i) = scalar%asString() tList_asStrings(i) = scalar%asString()
len_max = max(len_max, len_trim(tList_asStrings(i)))
item => item%next item => item%next
enddo enddo
!ToDo: trim to len_max
end function tList_asStrings end function tList_asStrings

22
src/constitutive.f90
View File

@ -127,7 +127,7 @@ module constitutive
instance,of,ip,el) instance,of,ip,el)
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
Mp !< MandelStress Mp !< MandelStress
real(pReal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: & real(pReal), dimension(3,3,homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems), intent(in) :: &
F, & !< deformation gradient F, & !< deformation gradient
Fp !< plastic deformation gradient Fp !< plastic deformation gradient
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
@ -218,7 +218,7 @@ module constitutive
instance, & instance, &
i, & i, &
e e
type(rotation), dimension(1,discretization_nIP,discretization_nElem), intent(in) :: & type(rotation), dimension(1,discretization_nIPs,discretization_Nelems), intent(in) :: &
orientation !< crystal orientation orientation !< crystal orientation
end subroutine plastic_nonlocal_updateCompatibility end subroutine plastic_nonlocal_updateCompatibility
@ -753,7 +753,7 @@ function constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el
of of
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
subdt !< timestep subdt !< timestep
real(pReal), intent(in), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: & real(pReal), intent(in), dimension(3,3,homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems) :: &
FArray, & !< elastic deformation gradient FArray, & !< elastic deformation gradient
FpArray !< plastic deformation gradient FpArray !< plastic deformation gradient
real(pReal), intent(in), dimension(3,3) :: & real(pReal), intent(in), dimension(3,3) :: &
@ -905,12 +905,12 @@ end function constitutive_deltaState
!> @brief Allocate the components of the state structure for a given phase !> @brief Allocate the components of the state structure for a given phase
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine constitutive_allocateState(state, & subroutine constitutive_allocateState(state, &
NipcMyPhase,sizeState,sizeDotState,sizeDeltaState) Nconstituents,sizeState,sizeDotState,sizeDeltaState)
class(tState), intent(out) :: & class(tState), intent(out) :: &
state state
integer, intent(in) :: & integer, intent(in) :: &
NipcMyPhase, & Nconstituents, &
sizeState, & sizeState, &
sizeDotState, & sizeDotState, &
sizeDeltaState sizeDeltaState
@ -921,14 +921,14 @@ subroutine constitutive_allocateState(state, &
state%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition state%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition
allocate(state%atol (sizeState), source=0.0_pReal) allocate(state%atol (sizeState), source=0.0_pReal)
allocate(state%state0 (sizeState,NipcMyPhase), source=0.0_pReal) allocate(state%state0 (sizeState,Nconstituents), source=0.0_pReal)
allocate(state%partitionedState0(sizeState,NipcMyPhase), source=0.0_pReal) allocate(state%partitionedState0(sizeState,Nconstituents), source=0.0_pReal)
allocate(state%subState0 (sizeState,NipcMyPhase), source=0.0_pReal) allocate(state%subState0 (sizeState,Nconstituents), source=0.0_pReal)
allocate(state%state (sizeState,NipcMyPhase), source=0.0_pReal) allocate(state%state (sizeState,Nconstituents), source=0.0_pReal)
allocate(state%dotState (sizeDotState,NipcMyPhase), source=0.0_pReal) allocate(state%dotState (sizeDotState,Nconstituents), source=0.0_pReal)
allocate(state%deltaState(sizeDeltaState,NipcMyPhase), source=0.0_pReal) allocate(state%deltaState(sizeDeltaState,Nconstituents), source=0.0_pReal)
end subroutine constitutive_allocateState end subroutine constitutive_allocateState

2
src/constitutive_damage.f90
View File

@ -184,7 +184,7 @@ module subroutine constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi
phiDot = 0.0_pReal phiDot = 0.0_pReal
dPhiDot_dPhi = 0.0_pReal dPhiDot_dPhi = 0.0_pReal
do grain = 1, homogenization_Ngrains(material_homogenizationAt(el)) do grain = 1, homogenization_Nconstituents(material_homogenizationAt(el))
phase = material_phaseAt(grain,el) phase = material_phaseAt(grain,el)
constituent = material_phasememberAt(grain,ip,el) constituent = material_phasememberAt(grain,ip,el)
do source = 1, phase_Nsources(phase) do source = 1, phase_Nsources(phase)

30
src/constitutive_plastic_disloTungsten.f90
View File

@ -78,9 +78,9 @@ module function plastic_disloTungsten_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity logical, dimension(:), allocatable :: myPlasticity
integer :: & integer :: &
Ninstance, & Ninstances, &
p, i, & p, i, &
NipcMyPhase, & Nconstituents, &
sizeState, sizeDotState, & sizeState, sizeDotState, &
startIndex, endIndex startIndex, endIndex
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
@ -99,17 +99,17 @@ module function plastic_disloTungsten_init() result(myPlasticity)
print'(/,a)', ' <<<+- plastic_dislotungsten init -+>>>' print'(/,a)', ' <<<+- plastic_dislotungsten init -+>>>'
myPlasticity = plastic_active('disloTungsten') myPlasticity = plastic_active('disloTungsten')
Ninstance = count(myPlasticity) Ninstances = count(myPlasticity)
print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstances == 0) return
print*, 'Cereceda et al., International Journal of Plasticity 78:242–256, 2016' print*, 'Cereceda et al., International Journal of Plasticity 78:242–256, 2016'
print*, 'https://dx.doi.org/10.1016/j.ijplas.2015.09.002' print*, 'https://dx.doi.org/10.1016/j.ijplas.2015.09.002'
allocate(param(Ninstance)) allocate(param(Ninstances))
allocate(state(Ninstance)) allocate(state(Ninstances))
allocate(dotState(Ninstance)) allocate(dotState(Ninstances))
allocate(dependentState(Ninstance)) allocate(dependentState(Ninstances))
phases => config_material%get('phase') phases => config_material%get('phase')
i = 0 i = 0
@ -221,18 +221,18 @@ module function plastic_disloTungsten_init() result(myPlasticity)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
NipcMyPhase = count(material_phaseAt == p) * discretization_nIP Nconstituents = count(material_phaseAt == p) * discretization_nIPs
sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl
sizeState = sizeDotState sizeState = sizeDotState
call constitutive_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,0) call constitutive_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! state aliases and initialization ! state aliases and initialization
startIndex = 1 startIndex = 1
endIndex = prm%sum_N_sl endIndex = prm%sum_N_sl
stt%rho_mob => plasticState(p)%state(startIndex:endIndex,:) stt%rho_mob => plasticState(p)%state(startIndex:endIndex,:)
stt%rho_mob = spread(rho_mob_0,2,NipcMyPhase) stt%rho_mob = spread(rho_mob_0,2,Nconstituents)
dot%rho_mob => plasticState(p)%dotState(startIndex:endIndex,:) dot%rho_mob => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal) plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
if (any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho' if (any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho'
@ -240,7 +240,7 @@ module function plastic_disloTungsten_init() result(myPlasticity)
startIndex = endIndex + 1 startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl endIndex = endIndex + prm%sum_N_sl
stt%rho_dip => plasticState(p)%state(startIndex:endIndex,:) stt%rho_dip => plasticState(p)%state(startIndex:endIndex,:)
stt%rho_dip = spread(rho_dip_0,2,NipcMyPhase) stt%rho_dip = spread(rho_dip_0,2,Nconstituents)
dot%rho_dip => plasticState(p)%dotState(startIndex:endIndex,:) dot%rho_dip => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal) plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
@ -252,8 +252,8 @@ module function plastic_disloTungsten_init() result(myPlasticity)
! global alias ! global alias
plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:) plasticState(p)%slipRate => plasticState(p)%dotState(startIndex:endIndex,:)
allocate(dst%Lambda_sl(prm%sum_N_sl,NipcMyPhase), source=0.0_pReal) allocate(dst%Lambda_sl(prm%sum_N_sl,Nconstituents), source=0.0_pReal)
allocate(dst%threshold_stress(prm%sum_N_sl,NipcMyPhase), source=0.0_pReal) allocate(dst%threshold_stress(prm%sum_N_sl,Nconstituents), source=0.0_pReal)
plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally

46
src/constitutive_plastic_dislotwin.f90
View File

@ -126,9 +126,9 @@ module function plastic_dislotwin_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity logical, dimension(:), allocatable :: myPlasticity
integer :: & integer :: &
Ninstance, & Ninstances, &
p, i, & p, i, &
NipcMyPhase, & Nconstituents, &
sizeState, sizeDotState, & sizeState, sizeDotState, &
startIndex, endIndex startIndex, endIndex
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
@ -146,9 +146,9 @@ module function plastic_dislotwin_init() result(myPlasticity)
print'(/,a)', ' <<<+- plastic_dislotwin init -+>>>' print'(/,a)', ' <<<+- plastic_dislotwin init -+>>>'
myPlasticity = plastic_active('dislotwin') myPlasticity = plastic_active('dislotwin')
Ninstance = count(myPlasticity) Ninstances = count(myPlasticity)
print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstances == 0) return
print*, 'Ma and Roters, Acta Materialia 52(12):3603–3612, 2004' print*, 'Ma and Roters, Acta Materialia 52(12):3603–3612, 2004'
print*, 'https://doi.org/10.1016/j.actamat.2004.04.012'//IO_EOL print*, 'https://doi.org/10.1016/j.actamat.2004.04.012'//IO_EOL
@ -159,10 +159,10 @@ module function plastic_dislotwin_init() result(myPlasticity)
print*, 'Wong et al., Acta Materialia 118:140–151, 2016' print*, 'Wong et al., Acta Materialia 118:140–151, 2016'
print*, 'https://doi.org/10.1016/j.actamat.2016.07.032' print*, 'https://doi.org/10.1016/j.actamat.2016.07.032'
allocate(param(Ninstance)) allocate(param(Ninstances))
allocate(state(Ninstance)) allocate(state(Ninstances))
allocate(dotState(Ninstance)) allocate(dotState(Ninstances))
allocate(dependentState(Ninstance)) allocate(dependentState(Ninstances))
phases => config_material%get('phase') phases => config_material%get('phase')
i = 0 i = 0
@ -407,21 +407,21 @@ module function plastic_dislotwin_init() result(myPlasticity)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
NipcMyPhase = count(material_phaseAt == p) * discretization_nIP Nconstituents = count(material_phaseAt == p) * discretization_nIPs
sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl & sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl &
+ size(['f_tw']) * prm%sum_N_tw & + size(['f_tw']) * prm%sum_N_tw &
+ size(['f_tr']) * prm%sum_N_tr + size(['f_tr']) * prm%sum_N_tr
sizeState = sizeDotState sizeState = sizeDotState
call constitutive_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,0) call constitutive_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! locally defined state aliases and initialization of state0 and atol ! locally defined state aliases and initialization of state0 and atol
startIndex = 1 startIndex = 1
endIndex = prm%sum_N_sl endIndex = prm%sum_N_sl
stt%rho_mob=>plasticState(p)%state(startIndex:endIndex,:) stt%rho_mob=>plasticState(p)%state(startIndex:endIndex,:)
stt%rho_mob= spread(rho_mob_0,2,NipcMyPhase) stt%rho_mob= spread(rho_mob_0,2,Nconstituents)
dot%rho_mob=>plasticState(p)%dotState(startIndex:endIndex,:) dot%rho_mob=>plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal) plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
if (any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho' if (any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_rho'
@ -429,7 +429,7 @@ module function plastic_dislotwin_init() result(myPlasticity)
startIndex = endIndex + 1 startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_sl endIndex = endIndex + prm%sum_N_sl
stt%rho_dip=>plasticState(p)%state(startIndex:endIndex,:) stt%rho_dip=>plasticState(p)%state(startIndex:endIndex,:)
stt%rho_dip= spread(rho_dip_0,2,NipcMyPhase) stt%rho_dip= spread(rho_dip_0,2,Nconstituents)
dot%rho_dip=>plasticState(p)%dotState(startIndex:endIndex,:) dot%rho_dip=>plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal) plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_rho',defaultVal=1.0_pReal)
@ -455,18 +455,18 @@ module function plastic_dislotwin_init() result(myPlasticity)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('f_trans',defaultVal=1.0e-6_pReal) plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('f_trans',defaultVal=1.0e-6_pReal)
if (any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' f_trans' if (any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' f_trans'
allocate(dst%Lambda_sl (prm%sum_N_sl,NipcMyPhase),source=0.0_pReal) allocate(dst%Lambda_sl (prm%sum_N_sl,Nconstituents),source=0.0_pReal)
allocate(dst%tau_pass (prm%sum_N_sl,NipcMyPhase),source=0.0_pReal) allocate(dst%tau_pass (prm%sum_N_sl,Nconstituents),source=0.0_pReal)
allocate(dst%Lambda_tw (prm%sum_N_tw,NipcMyPhase),source=0.0_pReal) allocate(dst%Lambda_tw (prm%sum_N_tw,Nconstituents),source=0.0_pReal)
allocate(dst%tau_hat_tw (prm%sum_N_tw,NipcMyPhase),source=0.0_pReal) allocate(dst%tau_hat_tw (prm%sum_N_tw,Nconstituents),source=0.0_pReal)
allocate(dst%tau_r_tw (prm%sum_N_tw,NipcMyPhase),source=0.0_pReal) allocate(dst%tau_r_tw (prm%sum_N_tw,Nconstituents),source=0.0_pReal)
allocate(dst%V_tw (prm%sum_N_tw,NipcMyPhase),source=0.0_pReal) allocate(dst%V_tw (prm%sum_N_tw,Nconstituents),source=0.0_pReal)
allocate(dst%Lambda_tr (prm%sum_N_tr,NipcMyPhase),source=0.0_pReal) allocate(dst%Lambda_tr (prm%sum_N_tr,Nconstituents),source=0.0_pReal)
allocate(dst%tau_hat_tr (prm%sum_N_tr,NipcMyPhase),source=0.0_pReal) allocate(dst%tau_hat_tr (prm%sum_N_tr,Nconstituents),source=0.0_pReal)
allocate(dst%tau_r_tr (prm%sum_N_tr,NipcMyPhase),source=0.0_pReal) allocate(dst%tau_r_tr (prm%sum_N_tr,Nconstituents),source=0.0_pReal)
allocate(dst%V_tr (prm%sum_N_tr,NipcMyPhase),source=0.0_pReal) allocate(dst%V_tr (prm%sum_N_tr,Nconstituents),source=0.0_pReal)
plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally plasticState(p)%state0 = plasticState(p)%state ! ToDo: this could be done centrally

20
src/constitutive_plastic_isotropic.f90
View File

@ -53,10 +53,10 @@ module function plastic_isotropic_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity logical, dimension(:), allocatable :: myPlasticity
integer :: & integer :: &
Ninstance, & Ninstances, &
p, & p, &
i, & i, &
NipcMyPhase, & Nconstituents, &
sizeState, sizeDotState sizeState, sizeDotState
real(pReal) :: & real(pReal) :: &
xi_0 !< initial critical stress xi_0 !< initial critical stress
@ -70,16 +70,16 @@ module function plastic_isotropic_init() result(myPlasticity)
print'(/,a)', ' <<<+- plastic_isotropic init -+>>>' print'(/,a)', ' <<<+- plastic_isotropic init -+>>>'
myPlasticity = plastic_active('isotropic') myPlasticity = plastic_active('isotropic')
Ninstance = count(myPlasticity) Ninstances = count(myPlasticity)
print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstances == 0) return
print*, 'Maiti and Eisenlohr, Scripta Materialia 145:37–40, 2018' print*, 'Maiti and Eisenlohr, Scripta Materialia 145:37–40, 2018'
print*, 'https://doi.org/10.1016/j.scriptamat.2017.09.047' print*, 'https://doi.org/10.1016/j.scriptamat.2017.09.047'
allocate(param(Ninstance)) allocate(param(Ninstances))
allocate(state(Ninstance)) allocate(state(Ninstances))
allocate(dotState(Ninstance)) allocate(dotState(Ninstances))
phases => config_material%get('phase') phases => config_material%get('phase')
i = 0 i = 0
@ -130,11 +130,11 @@ module function plastic_isotropic_init() result(myPlasticity)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
NipcMyPhase = count(material_phaseAt == p) * discretization_nIP Nconstituents = count(material_phaseAt == p) * discretization_nIPs
sizeDotState = size(['xi ','gamma']) sizeDotState = size(['xi ','gamma'])
sizeState = sizeDotState sizeState = sizeDotState
call constitutive_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,0) call constitutive_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! state aliases and initialization ! state aliases and initialization

24
src/constitutive_plastic_kinehardening.f90
View File

@ -62,9 +62,9 @@ module function plastic_kinehardening_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity logical, dimension(:), allocatable :: myPlasticity
integer :: & integer :: &
Ninstance, & Ninstances, &
p, i, o, & p, i, o, &
NipcMyPhase, & Nconstituents, &
sizeState, sizeDeltaState, sizeDotState, & sizeState, sizeDeltaState, sizeDotState, &
startIndex, endIndex startIndex, endIndex
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
@ -82,14 +82,14 @@ module function plastic_kinehardening_init() result(myPlasticity)
print'(/,a)', ' <<<+- plastic_kinehardening init -+>>>' print'(/,a)', ' <<<+- plastic_kinehardening init -+>>>'
myPlasticity = plastic_active('kinehardening') myPlasticity = plastic_active('kinehardening')
Ninstance = count(myPlasticity) Ninstances = count(myPlasticity)
print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstances == 0) return
allocate(param(Ninstance)) allocate(param(Ninstances))
allocate(state(Ninstance)) allocate(state(Ninstances))
allocate(dotState(Ninstance)) allocate(dotState(Ninstances))
allocate(deltaState(Ninstance)) allocate(deltaState(Ninstances))
phases => config_material%get('phase') phases => config_material%get('phase')
i = 0 i = 0
@ -174,19 +174,19 @@ module function plastic_kinehardening_init() result(myPlasticity)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
NipcMyPhase = count(material_phaseAt == p) * discretization_nIP Nconstituents = count(material_phaseAt == p) * discretization_nIPs
sizeDotState = size(['crss ','crss_back', 'accshear ']) * prm%sum_N_sl!ToDo: adjust names, ask Philip sizeDotState = size(['crss ','crss_back', 'accshear ']) * prm%sum_N_sl!ToDo: adjust names, ask Philip
sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%sum_N_sl !ToDo: adjust names sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%sum_N_sl !ToDo: adjust names
sizeState = sizeDotState + sizeDeltaState sizeState = sizeDotState + sizeDeltaState
call constitutive_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,sizeDeltaState) call constitutive_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,sizeDeltaState)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! state aliases and initialization ! state aliases and initialization
startIndex = 1 startIndex = 1
endIndex = prm%sum_N_sl endIndex = prm%sum_N_sl
stt%crss => plasticState(p)%state (startIndex:endIndex,:) stt%crss => plasticState(p)%state (startIndex:endIndex,:)
stt%crss = spread(xi_0, 2, NipcMyPhase) stt%crss = spread(xi_0, 2, Nconstituents)
dot%crss => plasticState(p)%dotState(startIndex:endIndex,:) dot%crss => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi' if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'

14
src/constitutive_plastic_none.f90
View File

@ -16,9 +16,9 @@ module function plastic_none_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity logical, dimension(:), allocatable :: myPlasticity
integer :: & integer :: &
Ninstance, & Ninstances, &
p, & p, &
NipcMyPhase Nconstituents
class(tNode), pointer :: & class(tNode), pointer :: &
phases, & phases, &
phase, & phase, &
@ -34,15 +34,15 @@ module function plastic_none_init() result(myPlasticity)
if(pl%get_asString('type') == 'none') myPlasticity(p) = .true. if(pl%get_asString('type') == 'none') myPlasticity(p) = .true.
enddo enddo
Ninstance = count(myPlasticity) Ninstances = count(myPlasticity)
print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstances == 0) return
do p = 1, phases%length do p = 1, phases%length
phase => phases%get(p) phase => phases%get(p)
if(.not. myPlasticity(p)) cycle if(.not. myPlasticity(p)) cycle
NipcMyPhase = count(material_phaseAt == p) * discretization_nIP Nconstituents = count(material_phaseAt == p) * discretization_nIPs
call constitutive_allocateState(plasticState(p),NipcMyPhase,0,0,0) call constitutive_allocateState(plasticState(p),Nconstituents,0,0,0)
enddo enddo
end function plastic_none_init end function plastic_none_init

84
src/constitutive_plastic_nonlocal.f90
View File

@ -153,7 +153,7 @@ submodule(constitutive:constitutive_plastic) plastic_nonlocal
state, & state, &
state0 state0
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances)
type(tNonlocalMicrostructure), dimension(:), allocatable :: microstructure type(tNonlocalMicrostructure), dimension(:), allocatable :: microstructure
@ -168,9 +168,9 @@ module function plastic_nonlocal_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity logical, dimension(:), allocatable :: myPlasticity
integer :: & integer :: &
Ninstance, & Ninstances, &
p, i, & p, i, &
NipcMyPhase, & Nconstituents, &
sizeState, sizeDotState, sizeDependentState, sizeDeltaState, & sizeState, sizeDotState, sizeDependentState, sizeDeltaState, &
s1, s2, & s1, s2, &
s, t, l s, t, l
@ -188,9 +188,9 @@ module function plastic_nonlocal_init() result(myPlasticity)
print'(/,a)', ' <<<+- plastic_nonlocal init -+>>>' print'(/,a)', ' <<<+- plastic_nonlocal init -+>>>'
myPlasticity = plastic_active('nonlocal') myPlasticity = plastic_active('nonlocal')
Ninstance = count(myPlasticity) Ninstances = count(myPlasticity)
print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstance == 0) then if(Ninstances == 0) then
call geometry_plastic_nonlocal_disable call geometry_plastic_nonlocal_disable
return return
endif endif
@ -201,12 +201,12 @@ module function plastic_nonlocal_init() result(myPlasticity)
print*, 'Kords, Dissertation RWTH Aachen, 2014' print*, 'Kords, Dissertation RWTH Aachen, 2014'
print*, 'http://publications.rwth-aachen.de/record/229993' print*, 'http://publications.rwth-aachen.de/record/229993'
allocate(param(Ninstance)) allocate(param(Ninstances))
allocate(state(Ninstance)) allocate(state(Ninstances))
allocate(state0(Ninstance)) allocate(state0(Ninstances))
allocate(dotState(Ninstance)) allocate(dotState(Ninstances))
allocate(deltaState(Ninstance)) allocate(deltaState(Ninstances))
allocate(microstructure(Ninstance)) allocate(microstructure(Ninstances))
phases => config_material%get('phase') phases => config_material%get('phase')
i = 0 i = 0
@ -391,7 +391,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
NipcMyPhase = count(material_phaseAt==p) * discretization_nIP Nconstituents = count(material_phaseAt==p) * discretization_nIPs
sizeDotState = size([ 'rhoSglEdgePosMobile ','rhoSglEdgeNegMobile ', & sizeDotState = size([ 'rhoSglEdgePosMobile ','rhoSglEdgeNegMobile ', &
'rhoSglScrewPosMobile ','rhoSglScrewNegMobile ', & 'rhoSglScrewPosMobile ','rhoSglScrewNegMobile ', &
'rhoSglEdgePosImmobile ','rhoSglEdgeNegImmobile ', & 'rhoSglEdgePosImmobile ','rhoSglEdgeNegImmobile ', &
@ -405,7 +405,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
'maxDipoleHeightEdge ','maxDipoleHeightScrew' ]) * prm%sum_N_sl !< other dependent state variables that are not updated by microstructure 'maxDipoleHeightEdge ','maxDipoleHeightScrew' ]) * prm%sum_N_sl !< other dependent state variables that are not updated by microstructure
sizeDeltaState = sizeDotState sizeDeltaState = sizeDotState
call constitutive_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,sizeDeltaState) call constitutive_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,sizeDeltaState)
plasticState(p)%nonlocal = pl%get_asBool('nonlocal') plasticState(p)%nonlocal = pl%get_asBool('nonlocal')
if(plasticState(p)%nonlocal .and. .not. allocated(IPneighborhood)) & if(plasticState(p)%nonlocal .and. .not. allocated(IPneighborhood)) &
@ -476,26 +476,26 @@ module function plastic_nonlocal_init() result(myPlasticity)
dot%rho_dip_scr => plasticState(p)%dotState (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:) dot%rho_dip_scr => plasticState(p)%dotState (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
del%rho_dip_scr => plasticState(p)%deltaState (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:) del%rho_dip_scr => plasticState(p)%deltaState (9*prm%sum_N_sl+1:10*prm%sum_N_sl,:)
stt%gamma => plasticState(p)%state (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:NipcMyPhase) stt%gamma => plasticState(p)%state (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents)
dot%gamma => plasticState(p)%dotState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:NipcMyPhase) dot%gamma => plasticState(p)%dotState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents)
del%gamma => plasticState(p)%deltaState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:NipcMyPhase) del%gamma => plasticState(p)%deltaState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents)
plasticState(p)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl ) = pl%get_asFloat('atol_gamma', defaultVal = 1.0e-2_pReal) plasticState(p)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl ) = pl%get_asFloat('atol_gamma', defaultVal = 1.0e-2_pReal)
if(any(plasticState(p)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl) < 0.0_pReal)) & if(any(plasticState(p)%atol(10*prm%sum_N_sl+1:11*prm%sum_N_sl) < 0.0_pReal)) &
extmsg = trim(extmsg)//' atol_gamma' extmsg = trim(extmsg)//' atol_gamma'
plasticState(p)%slipRate => plasticState(p)%dotState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:NipcMyPhase) plasticState(p)%slipRate => plasticState(p)%dotState (10*prm%sum_N_sl + 1:11*prm%sum_N_sl,1:Nconstituents)
stt%rho_forest => plasticState(p)%state (11*prm%sum_N_sl + 1:12*prm%sum_N_sl,1:NipcMyPhase) stt%rho_forest => plasticState(p)%state (11*prm%sum_N_sl + 1:12*prm%sum_N_sl,1:Nconstituents)
stt%v => plasticState(p)%state (12*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:NipcMyPhase) stt%v => plasticState(p)%state (12*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nconstituents)
stt%v_edg_pos => plasticState(p)%state (12*prm%sum_N_sl + 1:13*prm%sum_N_sl,1:NipcMyPhase) stt%v_edg_pos => plasticState(p)%state (12*prm%sum_N_sl + 1:13*prm%sum_N_sl,1:Nconstituents)
stt%v_edg_neg => plasticState(p)%state (13*prm%sum_N_sl + 1:14*prm%sum_N_sl,1:NipcMyPhase) stt%v_edg_neg => plasticState(p)%state (13*prm%sum_N_sl + 1:14*prm%sum_N_sl,1:Nconstituents)
stt%v_scr_pos => plasticState(p)%state (14*prm%sum_N_sl + 1:15*prm%sum_N_sl,1:NipcMyPhase) stt%v_scr_pos => plasticState(p)%state (14*prm%sum_N_sl + 1:15*prm%sum_N_sl,1:Nconstituents)
stt%v_scr_neg => plasticState(p)%state (15*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:NipcMyPhase) stt%v_scr_neg => plasticState(p)%state (15*prm%sum_N_sl + 1:16*prm%sum_N_sl,1:Nconstituents)
allocate(dst%tau_pass(prm%sum_N_sl,NipcMyPhase),source=0.0_pReal) allocate(dst%tau_pass(prm%sum_N_sl,Nconstituents),source=0.0_pReal)
allocate(dst%tau_back(prm%sum_N_sl,NipcMyPhase),source=0.0_pReal) allocate(dst%tau_back(prm%sum_N_sl,Nconstituents),source=0.0_pReal)
end associate end associate
if (NipcMyPhase > 0) call stateInit(ini,p,NipcMyPhase,i) if (Nconstituents > 0) call stateInit(ini,p,Nconstituents,i)
plasticState(p)%state0 = plasticState(p)%state plasticState(p)%state0 = plasticState(p)%state
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -505,12 +505,12 @@ module function plastic_nonlocal_init() result(myPlasticity)
enddo enddo
allocate(compatibility(2,maxval(param%sum_N_sl),maxval(param%sum_N_sl),nIPneighbors,& allocate(compatibility(2,maxval(param%sum_N_sl),maxval(param%sum_N_sl),nIPneighbors,&
discretization_nIP,discretization_nElem), source=0.0_pReal) discretization_nIPs,discretization_Nelems), source=0.0_pReal)
! BEGIN DEPRECATED---------------------------------------------------------------------------------- ! BEGIN DEPRECATED----------------------------------------------------------------------------------
allocate(iRhoU(maxval(param%sum_N_sl),4,Ninstance), source=0) allocate(iRhoU(maxval(param%sum_N_sl),4,Ninstances), source=0)
allocate(iV(maxval(param%sum_N_sl),4,Ninstance), source=0) allocate(iV(maxval(param%sum_N_sl),4,Ninstances), source=0)
allocate(iD(maxval(param%sum_N_sl),2,Ninstance), source=0) allocate(iD(maxval(param%sum_N_sl),2,Ninstances), source=0)
i = 0 i = 0
do p = 1, phases%length do p = 1, phases%length
@ -519,7 +519,7 @@ module function plastic_nonlocal_init() result(myPlasticity)
if(.not. myPlasticity(p)) cycle if(.not. myPlasticity(p)) cycle
i = i + 1 i = i + 1
NipcMyPhase = count(material_phaseAt==p) * discretization_nIP Nconstituents = count(material_phaseAt==p) * discretization_nIPs
l = 0 l = 0
do t = 1,4 do t = 1,4
do s = 1,param(i)%sum_N_sl do s = 1,param(i)%sum_N_sl
@ -976,7 +976,7 @@ module subroutine plastic_nonlocal_dotState(Mp, F, Fp, Temperature,timestep, &
real(pReal), dimension(3,3), intent(in) :: & real(pReal), dimension(3,3), intent(in) :: &
Mp !< MandelStress Mp !< MandelStress
real(pReal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: & real(pReal), dimension(3,3,homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems), intent(in) :: &
F, & !< elastic deformation gradient F, & !< elastic deformation gradient
Fp !< plastic deformation gradient Fp !< plastic deformation gradient
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
@ -1176,7 +1176,7 @@ end subroutine plastic_nonlocal_dotState
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
function rhoDotFlux(F,Fp,timestep, instance,of,ip,el) function rhoDotFlux(F,Fp,timestep, instance,of,ip,el)
real(pReal), dimension(3,3,homogenization_maxNgrains,discretization_nIP,discretization_nElem), intent(in) :: & real(pReal), dimension(3,3,homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems), intent(in) :: &
F, & !< elastic deformation gradient F, & !< elastic deformation gradient
Fp !< plastic deformation gradient Fp !< plastic deformation gradient
real(pReal), intent(in) :: & real(pReal), intent(in) :: &
@ -1416,7 +1416,7 @@ end function rhoDotFlux
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e) module subroutine plastic_nonlocal_updateCompatibility(orientation,instance,i,e)
type(rotation), dimension(1,discretization_nIP,discretization_nElem), intent(in) :: & type(rotation), dimension(1,discretization_nIPs,discretization_Nelems), intent(in) :: &
orientation ! crystal orientation orientation ! crystal orientation
integer, intent(in) :: & integer, intent(in) :: &
instance, & instance, &
@ -1601,13 +1601,13 @@ end subroutine plastic_nonlocal_results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief populates the initial dislocation density !> @brief populates the initial dislocation density
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine stateInit(ini,phase,NipcMyPhase,instance) subroutine stateInit(ini,phase,Nconstituents,instance)
type(tInitialParameters) :: & type(tInitialParameters) :: &
ini ini
integer,intent(in) :: & integer,intent(in) :: &
phase, & phase, &
NipcMyPhase, & Nconstituents, &
instance instance
integer :: & integer :: &
e, & e, &
@ -1625,15 +1625,15 @@ subroutine stateInit(ini,phase,NipcMyPhase,instance)
totalVolume, & totalVolume, &
densityBinning, & densityBinning, &
minimumIpVolume minimumIpVolume
real(pReal), dimension(NipcMyPhase) :: & real(pReal), dimension(Nconstituents) :: &
volume volume
associate(stt => state(instance)) associate(stt => state(instance))
if (ini%random_rho_u > 0.0_pReal) then ! randomly distribute dislocation segments on random slip system and of random type in the volume if (ini%random_rho_u > 0.0_pReal) then ! randomly distribute dislocation segments on random slip system and of random type in the volume
do e = 1,discretization_nElem do e = 1,discretization_Nelems
do i = 1,discretization_nIP do i = 1,discretization_nIPs
if (material_phaseAt(1,e) == phase) volume(material_phasememberAt(1,i,e)) = IPvolume(i,e) if (material_phaseAt(1,e) == phase) volume(material_phasememberAt(1,i,e)) = IPvolume(i,e)
enddo enddo
enddo enddo
@ -1645,13 +1645,13 @@ subroutine stateInit(ini,phase,NipcMyPhase,instance)
meanDensity = 0.0_pReal meanDensity = 0.0_pReal
do while(meanDensity < ini%random_rho_u) do while(meanDensity < ini%random_rho_u)
call random_number(rnd) call random_number(rnd)
phasemember = nint(rnd(1)*real(NipcMyPhase,pReal) + 0.5_pReal) phasemember = nint(rnd(1)*real(Nconstituents,pReal) + 0.5_pReal)
s = nint(rnd(2)*real(sum(ini%N_sl),pReal)*4.0_pReal + 0.5_pReal) s = nint(rnd(2)*real(sum(ini%N_sl),pReal)*4.0_pReal + 0.5_pReal)
meanDensity = meanDensity + densityBinning * volume(phasemember) / totalVolume meanDensity = meanDensity + densityBinning * volume(phasemember) / totalVolume
stt%rhoSglMobile(s,phasemember) = densityBinning stt%rhoSglMobile(s,phasemember) = densityBinning
enddo enddo
else ! homogeneous distribution with noise else ! homogeneous distribution with noise
do e = 1, NipcMyPhase do e = 1, Nconstituents
do f = 1,size(ini%N_sl,1) do f = 1,size(ini%N_sl,1)
from = 1 + sum(ini%N_sl(1:f-1)) from = 1 + sum(ini%N_sl(1:f-1))
upto = sum(ini%N_sl(1:f)) upto = sum(ini%N_sl(1:f))

24
src/constitutive_plastic_phenopowerlaw.f90
View File

@ -70,9 +70,9 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
logical, dimension(:), allocatable :: myPlasticity logical, dimension(:), allocatable :: myPlasticity
integer :: & integer :: &
Ninstance, & Ninstances, &
p, i, & p, i, &
NipcMyPhase, & Nconstituents, &
sizeState, sizeDotState, & sizeState, sizeDotState, &
startIndex, endIndex startIndex, endIndex
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
@ -91,13 +91,13 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
print'(/,a)', ' <<<+- plastic_phenopowerlaw init -+>>>' print'(/,a)', ' <<<+- plastic_phenopowerlaw init -+>>>'
myPlasticity = plastic_active('phenopowerlaw') myPlasticity = plastic_active('phenopowerlaw')
Ninstance = count(myPlasticity) Ninstances = count(myPlasticity)
print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstances == 0) return
allocate(param(Ninstance)) allocate(param(Ninstances))
allocate(state(Ninstance)) allocate(state(Ninstances))
allocate(dotState(Ninstance)) allocate(dotState(Ninstances))
phases => config_material%get('phase') phases => config_material%get('phase')
i = 0 i = 0
@ -224,20 +224,20 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate state arrays ! allocate state arrays
NipcMyPhase = count(material_phaseAt == p) * discretization_nIP Nconstituents = count(material_phaseAt == p) * discretization_nIPs
sizeDotState = size(['xi_sl ','gamma_sl']) * prm%sum_N_sl & sizeDotState = size(['xi_sl ','gamma_sl']) * prm%sum_N_sl &
+ size(['xi_tw ','gamma_tw']) * prm%sum_N_tw + size(['xi_tw ','gamma_tw']) * prm%sum_N_tw
sizeState = sizeDotState sizeState = sizeDotState
call constitutive_allocateState(plasticState(p),NipcMyPhase,sizeState,sizeDotState,0) call constitutive_allocateState(plasticState(p),Nconstituents,sizeState,sizeDotState,0)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! state aliases and initialization ! state aliases and initialization
startIndex = 1 startIndex = 1
endIndex = prm%sum_N_sl endIndex = prm%sum_N_sl
stt%xi_slip => plasticState(p)%state (startIndex:endIndex,:) stt%xi_slip => plasticState(p)%state (startIndex:endIndex,:)
stt%xi_slip = spread(xi_0_sl, 2, NipcMyPhase) stt%xi_slip = spread(xi_0_sl, 2, Nconstituents)
dot%xi_slip => plasticState(p)%dotState(startIndex:endIndex,:) dot%xi_slip => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi' if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'
@ -245,7 +245,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity)
startIndex = endIndex + 1 startIndex = endIndex + 1
endIndex = endIndex + prm%sum_N_tw endIndex = endIndex + prm%sum_N_tw
stt%xi_twin => plasticState(p)%state (startIndex:endIndex,:) stt%xi_twin => plasticState(p)%state (startIndex:endIndex,:)
stt%xi_twin = spread(xi_0_tw, 2, NipcMyPhase) stt%xi_twin = spread(xi_0_tw, 2, Nconstituents)
dot%xi_twin => plasticState(p)%dotState(startIndex:endIndex,:) dot%xi_twin => plasticState(p)%dotState(startIndex:endIndex,:)
plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal) plasticState(p)%atol(startIndex:endIndex) = pl%get_asFloat('atol_xi',defaultVal=1.0_pReal)
if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi' if(any(plasticState(p)%atol(startIndex:endIndex) < 0.0_pReal)) extmsg = trim(extmsg)//' atol_xi'

2
src/constitutive_thermal.f90
View File

@ -95,7 +95,7 @@ module subroutine constitutive_thermal_getRateAndItsTangents(TDot, dTDot_dT, T,
homog = material_homogenizationAt(el) homog = material_homogenizationAt(el)
instance = thermal_typeInstance(homog) instance = thermal_typeInstance(homog)
do grain = 1, homogenization_Ngrains(homog) do grain = 1, homogenization_Nconstituents(homog)
phase = material_phaseAt(grain,el) phase = material_phaseAt(grain,el)
constituent = material_phasememberAt(grain,ip,el) constituent = material_phasememberAt(grain,ip,el)
do source = 1, phase_Nsources(phase) do source = 1, phase_Nsources(phase)

44
src/crystallite.f90
View File

@ -135,7 +135,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine crystallite_init subroutine crystallite_init
logical, dimension(discretization_nIP,discretization_nElem) :: devNull logical, dimension(discretization_nIPs,discretization_Nelems) :: devNull
integer :: & integer :: &
c, & !< counter in integration point component loop c, & !< counter in integration point component loop
i, & !< counter in integration point loop i, & !< counter in integration point loop
@ -162,9 +162,9 @@ subroutine crystallite_init
debugCrystallite%ip = config_debug%get_asInt('integrationpoint', defaultVal=1) debugCrystallite%ip = config_debug%get_asInt('integrationpoint', defaultVal=1)
debugCrystallite%grain = config_debug%get_asInt('grain', defaultVal=1) debugCrystallite%grain = config_debug%get_asInt('grain', defaultVal=1)
cMax = homogenization_maxNgrains cMax = homogenization_maxNconstituents
iMax = discretization_nIP iMax = discretization_nIPs
eMax = discretization_nElem eMax = discretization_Nelems
allocate(crystallite_partitionedF(3,3,cMax,iMax,eMax),source=0.0_pReal) allocate(crystallite_partitionedF(3,3,cMax,iMax,eMax),source=0.0_pReal)
@ -253,7 +253,7 @@ subroutine crystallite_init
! initialize ! initialize
!$OMP PARALLEL DO PRIVATE(i,c) !$OMP PARALLEL DO PRIVATE(i,c)
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1), FEsolving_execIP(2); do c = 1, homogenization_Ngrains(material_homogenizationAt(e)) do i = FEsolving_execIP(1), FEsolving_execIP(2); do c = 1, homogenization_Nconstituents(material_homogenizationAt(e))
crystallite_Fp0(1:3,1:3,c,i,e) = material_orientation0(c,i,e)%asMatrix() ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005) crystallite_Fp0(1:3,1:3,c,i,e) = material_orientation0(c,i,e)%asMatrix() ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005)
crystallite_Fp0(1:3,1:3,c,i,e) = crystallite_Fp0(1:3,1:3,c,i,e) & crystallite_Fp0(1:3,1:3,c,i,e) = crystallite_Fp0(1:3,1:3,c,i,e) &
/ math_det33(crystallite_Fp0(1:3,1:3,c,i,e))**(1.0_pReal/3.0_pReal) / math_det33(crystallite_Fp0(1:3,1:3,c,i,e))**(1.0_pReal/3.0_pReal)
@ -279,7 +279,7 @@ subroutine crystallite_init
!$OMP PARALLEL DO !$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2) do i = FEsolving_execIP(1),FEsolving_execIP(2)
do c = 1,homogenization_Ngrains(material_homogenizationAt(e)) do c = 1,homogenization_Nconstituents(material_homogenizationAt(e))
call constitutive_dependentState(crystallite_partitionedF0(1:3,1:3,c,i,e), & call constitutive_dependentState(crystallite_partitionedF0(1:3,1:3,c,i,e), &
crystallite_partitionedFp0(1:3,1:3,c,i,e), & crystallite_partitionedFp0(1:3,1:3,c,i,e), &
c,i,e) ! update dependent state variables to be consistent with basic states c,i,e) ! update dependent state variables to be consistent with basic states
@ -307,7 +307,7 @@ end subroutine crystallite_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
function crystallite_stress() function crystallite_stress()
logical, dimension(discretization_nIP,discretization_nElem) :: crystallite_stress logical, dimension(discretization_nIPs,discretization_Nelems) :: crystallite_stress
real(pReal) :: & real(pReal) :: &
formerSubStep formerSubStep
integer :: & integer :: &
@ -317,7 +317,7 @@ function crystallite_stress()
e, & !< counter in element loop e, & !< counter in element loop
startIP, endIP, & startIP, endIP, &
s s
logical, dimension(homogenization_maxNgrains,discretization_nIP,discretization_nElem) :: todo !ToDo: need to set some values to false for different Ngrains logical, dimension(homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems) :: todo !ToDo: need to set some values to false for different Ngrains
real(pReal), dimension(:,:,:,:,:), allocatable :: & real(pReal), dimension(:,:,:,:,:), allocatable :: &
subLp0,& !< plastic velocity grad at start of crystallite inc subLp0,& !< plastic velocity grad at start of crystallite inc
subLi0 !< intermediate velocity grad at start of crystallite inc subLi0 !< intermediate velocity grad at start of crystallite inc
@ -335,7 +335,7 @@ function crystallite_stress()
crystallite_subStep = 0.0_pReal crystallite_subStep = 0.0_pReal
!$OMP PARALLEL DO !$OMP PARALLEL DO
elementLooping1: do e = FEsolving_execElem(1),FEsolving_execElem(2) elementLooping1: do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2); do c = 1,homogenization_Ngrains(material_homogenizationAt(e)) do i = FEsolving_execIP(1),FEsolving_execIP(2); do c = 1,homogenization_Nconstituents(material_homogenizationAt(e))
homogenizationRequestsCalculation: if (crystallite_requested(c,i,e)) then homogenizationRequestsCalculation: if (crystallite_requested(c,i,e)) then
plasticState (material_phaseAt(c,e))%subState0( :,material_phaseMemberAt(c,i,e)) = & plasticState (material_phaseAt(c,e))%subState0( :,material_phaseMemberAt(c,i,e)) = &
plasticState (material_phaseAt(c,e))%partitionedState0(:,material_phaseMemberAt(c,i,e)) plasticState (material_phaseAt(c,e))%partitionedState0(:,material_phaseMemberAt(c,i,e))
@ -362,7 +362,7 @@ function crystallite_stress()
endIP = startIP endIP = startIP
else singleRun else singleRun
startIP = 1 startIP = 1
endIP = discretization_nIP endIP = discretization_nIPs
endif singleRun endif singleRun
NiterationCrystallite = 0 NiterationCrystallite = 0
@ -376,7 +376,7 @@ function crystallite_stress()
!$OMP PARALLEL DO PRIVATE(formerSubStep) !$OMP PARALLEL DO PRIVATE(formerSubStep)
elementLooping3: do e = FEsolving_execElem(1),FEsolving_execElem(2) elementLooping3: do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2) do i = FEsolving_execIP(1),FEsolving_execIP(2)
do c = 1,homogenization_Ngrains(material_homogenizationAt(e)) do c = 1,homogenization_Nconstituents(material_homogenizationAt(e))
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! wind forward ! wind forward
if (crystallite_converged(c,i,e)) then if (crystallite_converged(c,i,e)) then
@ -472,7 +472,7 @@ subroutine crystallite_initializeRestorationPoints(i,e)
c, & !< constituent number c, & !< constituent number
s s
do c = 1,homogenization_Ngrains(material_homogenizationAt(e)) do c = 1,homogenization_Nconstituents(material_homogenizationAt(e))
crystallite_partitionedFp0(1:3,1:3,c,i,e) = crystallite_Fp0(1:3,1:3,c,i,e) crystallite_partitionedFp0(1:3,1:3,c,i,e) = crystallite_Fp0(1:3,1:3,c,i,e)
crystallite_partitionedLp0(1:3,1:3,c,i,e) = crystallite_Lp0(1:3,1:3,c,i,e) crystallite_partitionedLp0(1:3,1:3,c,i,e) = crystallite_Lp0(1:3,1:3,c,i,e)
crystallite_partitionedFi0(1:3,1:3,c,i,e) = crystallite_Fi0(1:3,1:3,c,i,e) crystallite_partitionedFi0(1:3,1:3,c,i,e) = crystallite_Fi0(1:3,1:3,c,i,e)
@ -503,7 +503,7 @@ subroutine crystallite_windForward(i,e)
c, & !< constituent number c, & !< constituent number
s s
do c = 1,homogenization_Ngrains(material_homogenizationAt(e)) do c = 1,homogenization_Nconstituents(material_homogenizationAt(e))
crystallite_partitionedF0 (1:3,1:3,c,i,e) = crystallite_partitionedF(1:3,1:3,c,i,e) crystallite_partitionedF0 (1:3,1:3,c,i,e) = crystallite_partitionedF(1:3,1:3,c,i,e)
crystallite_partitionedFp0(1:3,1:3,c,i,e) = crystallite_Fp (1:3,1:3,c,i,e) crystallite_partitionedFp0(1:3,1:3,c,i,e) = crystallite_Fp (1:3,1:3,c,i,e)
crystallite_partitionedLp0(1:3,1:3,c,i,e) = crystallite_Lp (1:3,1:3,c,i,e) crystallite_partitionedLp0(1:3,1:3,c,i,e) = crystallite_Lp (1:3,1:3,c,i,e)
@ -536,7 +536,7 @@ subroutine crystallite_restore(i,e,includeL)
c, & !< constituent number c, & !< constituent number
s s
do c = 1,homogenization_Ngrains(material_homogenizationAt(e)) do c = 1,homogenization_Nconstituents(material_homogenizationAt(e))
if (includeL) then if (includeL) then
crystallite_Lp(1:3,1:3,c,i,e) = crystallite_partitionedLp0(1:3,1:3,c,i,e) crystallite_Lp(1:3,1:3,c,i,e) = crystallite_partitionedLp0(1:3,1:3,c,i,e)
crystallite_Li(1:3,1:3,c,i,e) = crystallite_partitionedLi0(1:3,1:3,c,i,e) crystallite_Li(1:3,1:3,c,i,e) = crystallite_partitionedLi0(1:3,1:3,c,i,e)
@ -697,7 +697,7 @@ subroutine crystallite_orientations
!$OMP PARALLEL DO !$OMP PARALLEL DO
do e = FEsolving_execElem(1),FEsolving_execElem(2) do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1),FEsolving_execIP(2) do i = FEsolving_execIP(1),FEsolving_execIP(2)
do c = 1,homogenization_Ngrains(material_homogenizationAt(e)) do c = 1,homogenization_Nconstituents(material_homogenizationAt(e))
call crystallite_orientation(c,i,e)%fromMatrix(transpose(math_rotationalPart(crystallite_Fe(1:3,1:3,c,i,e)))) call crystallite_orientation(c,i,e)%fromMatrix(transpose(math_rotationalPart(crystallite_Fe(1:3,1:3,c,i,e))))
enddo; enddo; enddo enddo; enddo; enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
@ -821,11 +821,11 @@ subroutine crystallite_results
real(pReal), allocatable, dimension(:,:,:) :: select_tensors real(pReal), allocatable, dimension(:,:,:) :: select_tensors
integer :: e,i,c,j integer :: e,i,c,j
allocate(select_tensors(3,3,count(material_phaseAt==instance)*discretization_nIP)) allocate(select_tensors(3,3,count(material_phaseAt==instance)*discretization_nIPs))
j=0 j=0
do e = 1, size(material_phaseAt,2) do e = 1, size(material_phaseAt,2)
do i = 1, discretization_nIP do i = 1, discretization_nIPs
do c = 1, size(material_phaseAt,1) !ToDo: this needs to be changed for varying Ngrains do c = 1, size(material_phaseAt,1) !ToDo: this needs to be changed for varying Ngrains
if (material_phaseAt(c,e) == instance) then if (material_phaseAt(c,e) == instance) then
j = j + 1 j = j + 1
@ -848,11 +848,11 @@ subroutine crystallite_results
type(rotation), allocatable, dimension(:) :: select_rotations type(rotation), allocatable, dimension(:) :: select_rotations
integer :: e,i,c,j integer :: e,i,c,j
allocate(select_rotations(count(material_phaseAt==instance)*homogenization_maxNgrains*discretization_nIP)) allocate(select_rotations(count(material_phaseAt==instance)*homogenization_maxNconstituents*discretization_nIPs))
j=0 j=0
do e = 1, size(material_phaseAt,2) do e = 1, size(material_phaseAt,2)
do i = 1, discretization_nIP do i = 1, discretization_nIPs
do c = 1, size(material_phaseAt,1) !ToDo: this needs to be changed for varying Ngrains do c = 1, size(material_phaseAt,1) !ToDo: this needs to be changed for varying Ngrains
if (material_phaseAt(c,e) == instance) then if (material_phaseAt(c,e) == instance) then
j = j + 1 j = j + 1
@ -1571,7 +1571,7 @@ subroutine crystallite_restartWrite
call HDF5_closeGroup(groupHandle) call HDF5_closeGroup(groupHandle)
groupHandle = HDF5_addGroup(fileHandle,'materialpoint') groupHandle = HDF5_addGroup(fileHandle,'materialpoint')
do i = 1, material_Nhomogenization do i = 1, size(material_name_homogenization)
write(datasetName,'(i0,a)') i,'_omega_homogenization' write(datasetName,'(i0,a)') i,'_omega_homogenization'
call HDF5_write(groupHandle,homogState(i)%state,datasetName) call HDF5_write(groupHandle,homogState(i)%state,datasetName)
enddo enddo
@ -1612,7 +1612,7 @@ subroutine crystallite_restartRead
call HDF5_closeGroup(groupHandle) call HDF5_closeGroup(groupHandle)
groupHandle = HDF5_openGroup(fileHandle,'materialpoint') groupHandle = HDF5_openGroup(fileHandle,'materialpoint')
do i = 1, material_Nhomogenization do i = 1,size(material_name_homogenization)
write(datasetName,'(i0,a)') i,'_omega_homogenization' write(datasetName,'(i0,a)') i,'_omega_homogenization'
call HDF5_read(groupHandle,homogState(i)%state0,datasetName) call HDF5_read(groupHandle,homogState(i)%state0,datasetName)
enddo enddo
@ -1645,7 +1645,7 @@ subroutine crystallite_forward
do j = 1,phase_Nsources(i) do j = 1,phase_Nsources(i)
sourceState(i)%p(j)%state0 = sourceState(i)%p(j)%state sourceState(i)%p(j)%state0 = sourceState(i)%p(j)%state
enddo; enddo enddo; enddo
do i = 1, material_Nhomogenization do i = 1,size(material_name_homogenization)
homogState (i)%state0 = homogState (i)%state homogState (i)%state0 = homogState (i)%state
thermalState(i)%state0 = thermalState(i)%state thermalState(i)%state0 = thermalState(i)%state
damageState (i)%state0 = damageState (i)%state damageState (i)%state0 = damageState (i)%state

18
src/damage_local.f90
View File

@ -42,7 +42,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine damage_local_init subroutine damage_local_init
integer :: Ninstance,NofMyHomog,h integer :: Ninstances,Nmaterialpoints,h
class(tNode), pointer :: & class(tNode), pointer :: &
num_generic, & num_generic, &
material_homogenization, & material_homogenization, &
@ -57,8 +57,8 @@ subroutine damage_local_init
num%residualStiffness = num_generic%get_asFloat('residualStiffness', defaultVal=1.0e-6_pReal) num%residualStiffness = num_generic%get_asFloat('residualStiffness', defaultVal=1.0e-6_pReal)
if (num%residualStiffness < 0.0_pReal) call IO_error(301,ext_msg='residualStiffness') if (num%residualStiffness < 0.0_pReal) call IO_error(301,ext_msg='residualStiffness')
Ninstance = count(damage_type == DAMAGE_local_ID) Ninstances = count(damage_type == DAMAGE_local_ID)
allocate(param(Ninstance)) allocate(param(Ninstances))
material_homogenization => config_material%get('homogenization') material_homogenization => config_material%get('homogenization')
do h = 1, material_homogenization%length do h = 1, material_homogenization%length
@ -73,11 +73,11 @@ subroutine damage_local_init
prm%output = homogDamage%get_asStrings('output',defaultVal=emptyStringArray) prm%output = homogDamage%get_asStrings('output',defaultVal=emptyStringArray)
#endif #endif
NofMyHomog = count(material_homogenizationAt == h) Nmaterialpoints = count(material_homogenizationAt == h)
damageState(h)%sizeState = 1 damageState(h)%sizeState = 1
allocate(damageState(h)%state0 (1,NofMyHomog), source=damage_initialPhi(h)) allocate(damageState(h)%state0 (1,Nmaterialpoints), source=damage_initialPhi(h))
allocate(damageState(h)%subState0(1,NofMyHomog), source=damage_initialPhi(h)) allocate(damageState(h)%subState0(1,Nmaterialpoints), source=damage_initialPhi(h))
allocate(damageState(h)%state (1,NofMyHomog), source=damage_initialPhi(h)) allocate(damageState(h)%state (1,Nmaterialpoints), source=damage_initialPhi(h))
nullify(damageMapping(h)%p) nullify(damageMapping(h)%p)
damageMapping(h)%p => material_homogenizationMemberAt damageMapping(h)%p => material_homogenizationMemberAt
@ -143,8 +143,8 @@ subroutine damage_local_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip, el
call constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) call constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el)
phiDot = phiDot/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) phiDot = phiDot/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)
dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)
end subroutine damage_local_getSourceAndItsTangent end subroutine damage_local_getSourceAndItsTangent

12
src/damage_none.f90
View File

@ -16,18 +16,18 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine damage_none_init subroutine damage_none_init
integer :: h,NofMyHomog integer :: h,Nmaterialpoints
print'(/,a)', ' <<<+- damage_none init -+>>>'; flush(6) print'(/,a)', ' <<<+- damage_none init -+>>>'; flush(6)
do h = 1, material_Nhomogenization do h = 1, size(material_name_homogenization)
if (damage_type(h) /= DAMAGE_NONE_ID) cycle if (damage_type(h) /= DAMAGE_NONE_ID) cycle
NofMyHomog = count(material_homogenizationAt == h) Nmaterialpoints = count(material_homogenizationAt == h)
damageState(h)%sizeState = 0 damageState(h)%sizeState = 0
allocate(damageState(h)%state0 (0,NofMyHomog)) allocate(damageState(h)%state0 (0,Nmaterialpoints))
allocate(damageState(h)%subState0(0,NofMyHomog)) allocate(damageState(h)%subState0(0,Nmaterialpoints))
allocate(damageState(h)%state (0,NofMyHomog)) allocate(damageState(h)%state (0,Nmaterialpoints))
deallocate(damage(h)%p) deallocate(damage(h)%p)
allocate (damage(h)%p(1), source=damage_initialPhi(h)) allocate (damage(h)%p(1), source=damage_initialPhi(h))

26
src/damage_nonlocal.f90
View File

@ -46,7 +46,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine damage_nonlocal_init subroutine damage_nonlocal_init
integer :: Ninstance,NofMyHomog,h integer :: Ninstances,Nmaterialpoints,h
class(tNode), pointer :: & class(tNode), pointer :: &
num_generic, & num_generic, &
material_homogenization, & material_homogenization, &
@ -60,8 +60,8 @@ subroutine damage_nonlocal_init
num_generic => config_numerics%get('generic',defaultVal= emptyDict) num_generic => config_numerics%get('generic',defaultVal= emptyDict)
num%charLength = num_generic%get_asFloat('charLength',defaultVal=1.0_pReal) num%charLength = num_generic%get_asFloat('charLength',defaultVal=1.0_pReal)
Ninstance = count(damage_type == DAMAGE_nonlocal_ID) Ninstances = count(damage_type == DAMAGE_nonlocal_ID)
allocate(param(Ninstance)) allocate(param(Ninstances))
material_homogenization => config_material%get('homogenization') material_homogenization => config_material%get('homogenization')
do h = 1, material_homogenization%length do h = 1, material_homogenization%length
@ -76,11 +76,11 @@ subroutine damage_nonlocal_init
prm%output = homogDamage%get_asStrings('output',defaultVal=emptyStringArray) prm%output = homogDamage%get_asStrings('output',defaultVal=emptyStringArray)
#endif #endif
NofMyHomog = count(material_homogenizationAt == h) Nmaterialpoints = count(material_homogenizationAt == h)
damageState(h)%sizeState = 1 damageState(h)%sizeState = 1
allocate(damageState(h)%state0 (1,NofMyHomog), source=damage_initialPhi(h)) allocate(damageState(h)%state0 (1,Nmaterialpoints), source=damage_initialPhi(h))
allocate(damageState(h)%subState0(1,NofMyHomog), source=damage_initialPhi(h)) allocate(damageState(h)%subState0(1,Nmaterialpoints), source=damage_initialPhi(h))
allocate(damageState(h)%state (1,NofMyHomog), source=damage_initialPhi(h)) allocate(damageState(h)%state (1,Nmaterialpoints), source=damage_initialPhi(h))
nullify(damageMapping(h)%p) nullify(damageMapping(h)%p)
damageMapping(h)%p => material_homogenizationMemberAt damageMapping(h)%p => material_homogenizationMemberAt
@ -110,8 +110,8 @@ subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, phi, ip,
dPhiDot_dPhi = 0.0_pReal dPhiDot_dPhi = 0.0_pReal
call constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el) call constitutive_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ip, el)
phiDot = phiDot/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) phiDot = phiDot/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)
dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)
end subroutine damage_nonlocal_getSourceAndItsTangent end subroutine damage_nonlocal_getSourceAndItsTangent
@ -132,13 +132,13 @@ function damage_nonlocal_getDiffusion(ip,el)
homog = material_homogenizationAt(el) homog = material_homogenizationAt(el)
damage_nonlocal_getDiffusion = 0.0_pReal damage_nonlocal_getDiffusion = 0.0_pReal
do grain = 1, homogenization_Ngrains(homog) do grain = 1, homogenization_Nconstituents(homog)
damage_nonlocal_getDiffusion = damage_nonlocal_getDiffusion + & damage_nonlocal_getDiffusion = damage_nonlocal_getDiffusion + &
crystallite_push33ToRef(grain,ip,el,lattice_D(1:3,1:3,material_phaseAt(grain,el))) crystallite_push33ToRef(grain,ip,el,lattice_D(1:3,1:3,material_phaseAt(grain,el)))
enddo enddo
damage_nonlocal_getDiffusion = & damage_nonlocal_getDiffusion = &
num%charLength**2*damage_nonlocal_getDiffusion/real(homogenization_Ngrains(homog),pReal) num%charLength**2*damage_nonlocal_getDiffusion/real(homogenization_Nconstituents(homog),pReal)
end function damage_nonlocal_getDiffusion end function damage_nonlocal_getDiffusion
@ -156,12 +156,12 @@ real(pReal) function damage_nonlocal_getMobility(ip,el)
damage_nonlocal_getMobility = 0.0_pReal damage_nonlocal_getMobility = 0.0_pReal
do ipc = 1, homogenization_Ngrains(material_homogenizationAt(el)) do ipc = 1, homogenization_Nconstituents(material_homogenizationAt(el))
damage_nonlocal_getMobility = damage_nonlocal_getMobility + lattice_M(material_phaseAt(ipc,el)) damage_nonlocal_getMobility = damage_nonlocal_getMobility + lattice_M(material_phaseAt(ipc,el))
enddo enddo
damage_nonlocal_getMobility = damage_nonlocal_getMobility/& damage_nonlocal_getMobility = damage_nonlocal_getMobility/&
real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)
end function damage_nonlocal_getMobility end function damage_nonlocal_getMobility

8
src/discretization.f90
View File

@ -11,8 +11,8 @@ module discretization
private private
integer, public, protected :: & integer, public, protected :: &
discretization_nIP, & discretization_nIPs, &
discretization_nElem discretization_Nelems
integer, public, protected, dimension(:), allocatable :: & integer, public, protected, dimension(:), allocatable :: &
discretization_materialAt discretization_materialAt
@ -51,8 +51,8 @@ subroutine discretization_init(materialAt,&
print'(/,a)', ' <<<+- discretization init -+>>>'; flush(6) print'(/,a)', ' <<<+- discretization init -+>>>'; flush(6)
discretization_nElem = size(materialAt,1) discretization_Nelems = size(materialAt,1)
discretization_nIP = size(IPcoords0,2)/discretization_nElem discretization_nIPs = size(IPcoords0,2)/discretization_Nelems
discretization_materialAt = materialAt discretization_materialAt = materialAt

2
src/grid/DAMASK_grid.f90
View File

@ -210,7 +210,7 @@ program DAMASK_grid
if(.not. step_discretization%contains('t')) call IO_error(error_ID=837,ext_msg = 't missing') if(.not. step_discretization%contains('t')) call IO_error(error_ID=837,ext_msg = 't missing')
if(.not. step_discretization%contains('N')) call IO_error(error_ID=837,ext_msg = 'N missing') if(.not. step_discretization%contains('N')) call IO_error(error_ID=837,ext_msg = 'N missing')
loadCases(l)%time = step_discretization%get_asFloat('t') loadCases(l)%time = step_discretization%get_asFloat('t')
loadCases(l)%incs = step_discretization%get_asFloat('N') loadCases(l)%incs = step_discretization%get_asInt ('N')
loadCases(l)%logscale = step_discretization%get_asBool ('log_timestep', defaultVal= .false.) loadCases(l)%logscale = step_discretization%get_asBool ('log_timestep', defaultVal= .false.)
loadCases(l)%outputfrequency = step_discretization%get_asInt ('f_out', defaultVal=1) loadCases(l)%outputfrequency = step_discretization%get_asInt ('f_out', defaultVal=1)
loadCases(l)%restartfrequency = step_discretization%get_asInt ('f_restart', defaultVal=huge(0)) loadCases(l)%restartfrequency = step_discretization%get_asInt ('f_restart', defaultVal=huge(0))

30
src/grid/grid_mech_FEM.f90
View File

@ -236,7 +236,7 @@ subroutine grid_mech_FEM_init
F = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3) F = spread(spread(spread(math_I3,3,grid(1)),4,grid(2)),5,grid3)
endif restartRead endif restartRead
materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent homogenization_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
call utilities_updateCoords(F) call utilities_updateCoords(F)
call utilities_constitutiveResponse(P_current,temp33_Real,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2 call utilities_constitutiveResponse(P_current,temp33_Real,C_volAvg,devNull, & ! stress field, stress avg, global average of stiffness and (min+max)/2
F, & ! target F F, & ! target F
@ -364,7 +364,7 @@ subroutine grid_mech_FEM_forward(cutBack,guess,timeinc,timeinc_old,loadCaseTime,
F_lastInc = F F_lastInc = F
materialpoint_F0 = reshape(F, [3,3,1,product(grid(1:2))*grid3]) homogenization_F0 = reshape(F, [3,3,1,product(grid(1:2))*grid3])
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -557,9 +557,9 @@ subroutine formResidual(da_local,x_local, &
ii = i-xstart+1; jj = j-ystart+1; kk = k-zstart+1 ii = i-xstart+1; jj = j-ystart+1; kk = k-zstart+1
ele = ele + 1 ele = ele + 1
f_elem = matmul(transpose(BMat),transpose(P_current(1:3,1:3,ii,jj,kk)))*detJ + & f_elem = matmul(transpose(BMat),transpose(P_current(1:3,1:3,ii,jj,kk)))*detJ + &
matmul(HGMat,x_elem)*(materialpoint_dPdF(1,1,1,1,1,ele) + & matmul(HGMat,x_elem)*(homogenization_dPdF(1,1,1,1,1,ele) + &
materialpoint_dPdF(2,2,2,2,1,ele) + & homogenization_dPdF(2,2,2,2,1,ele) + &
materialpoint_dPdF(3,3,3,3,1,ele))/3.0_pReal homogenization_dPdF(3,3,3,3,1,ele))/3.0_pReal
ctr = 0 ctr = 0
do kk = 0, 1; do jj = 0, 1; do ii = 0, 1 do kk = 0, 1; do jj = 0, 1; do ii = 0, 1
ctr = ctr + 1 ctr = ctr + 1
@ -636,18 +636,18 @@ subroutine formJacobian(da_local,x_local,Jac_pre,Jac,dummy,ierr)
row = col row = col
ele = ele + 1 ele = ele + 1
K_ele = 0.0 K_ele = 0.0
K_ele(1 :8 ,1 :8 ) = HGMat*(materialpoint_dPdF(1,1,1,1,1,ele) + & K_ele(1 :8 ,1 :8 ) = HGMat*(homogenization_dPdF(1,1,1,1,1,ele) + &
materialpoint_dPdF(2,2,2,2,1,ele) + & homogenization_dPdF(2,2,2,2,1,ele) + &
materialpoint_dPdF(3,3,3,3,1,ele))/3.0_pReal homogenization_dPdF(3,3,3,3,1,ele))/3.0_pReal
K_ele(9 :16,9 :16) = HGMat*(materialpoint_dPdF(1,1,1,1,1,ele) + & K_ele(9 :16,9 :16) = HGMat*(homogenization_dPdF(1,1,1,1,1,ele) + &
materialpoint_dPdF(2,2,2,2,1,ele) + & homogenization_dPdF(2,2,2,2,1,ele) + &
materialpoint_dPdF(3,3,3,3,1,ele))/3.0_pReal homogenization_dPdF(3,3,3,3,1,ele))/3.0_pReal
K_ele(17:24,17:24) = HGMat*(materialpoint_dPdF(1,1,1,1,1,ele) + & K_ele(17:24,17:24) = HGMat*(homogenization_dPdF(1,1,1,1,1,ele) + &
materialpoint_dPdF(2,2,2,2,1,ele) + & homogenization_dPdF(2,2,2,2,1,ele) + &
materialpoint_dPdF(3,3,3,3,1,ele))/3.0_pReal homogenization_dPdF(3,3,3,3,1,ele))/3.0_pReal
K_ele = K_ele + & K_ele = K_ele + &
matmul(transpose(BMatFull), & matmul(transpose(BMatFull), &
matmul(reshape(reshape(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,ele), & matmul(reshape(reshape(homogenization_dPdF(1:3,1:3,1:3,1:3,1,ele), &
shape=[3,3,3,3], order=[2,1,4,3]),shape=[9,9]),BMatFull))*detJ shape=[3,3,3,3], order=[2,1,4,3]),shape=[9,9]),BMatFull))*detJ
call MatSetValuesStencil(Jac,24,row,24,col,K_ele,ADD_VALUES,ierr) call MatSetValuesStencil(Jac,24,row,24,col,K_ele,ADD_VALUES,ierr)
CHKERRQ(ierr) CHKERRQ(ierr)

4
src/grid/grid_mech_spectral_basic.f90
View File

@ -198,7 +198,7 @@ subroutine grid_mech_spectral_basic_init
F = reshape(F_lastInc,[9,grid(1),grid(2),grid3]) F = reshape(F_lastInc,[9,grid(1),grid(2),grid3])
endif restartRead endif restartRead
materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent homogenization_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
call utilities_updateCoords(reshape(F,shape(F_lastInc))) call utilities_updateCoords(reshape(F,shape(F_lastInc)))
call utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 call utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2
reshape(F,shape(F_lastInc)), & ! target F reshape(F,shape(F_lastInc)), & ! target F
@ -324,7 +324,7 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,timeinc,timeinc_old,lo
rotation_BC%rotate(F_aimDot,active=.true.)) rotation_BC%rotate(F_aimDot,active=.true.))
F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid3]) F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid3])
materialpoint_F0 = reshape(F, [3,3,1,product(grid(1:2))*grid3]) homogenization_F0 = reshape(F, [3,3,1,product(grid(1:2))*grid3])
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

6
src/grid/grid_mech_spectral_polarisation.f90
View File

@ -224,7 +224,7 @@ subroutine grid_mech_spectral_polarisation_init
F_tau_lastInc = 2.0_pReal*F_lastInc F_tau_lastInc = 2.0_pReal*F_lastInc
endif restartRead endif restartRead
materialpoint_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent homogenization_F0 = reshape(F_lastInc, [3,3,1,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
call utilities_updateCoords(reshape(F,shape(F_lastInc))) call utilities_updateCoords(reshape(F,shape(F_lastInc)))
call utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2 call utilities_constitutiveResponse(P,temp33_Real,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2
reshape(F,shape(F_lastInc)), & ! target F reshape(F,shape(F_lastInc)), & ! target F
@ -364,7 +364,7 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,timeinc,timeinc
F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3]) F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3])
F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid3]) F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid3])
materialpoint_F0 = reshape(F,[3,3,1,product(grid(1:2))*grid3]) homogenization_F0 = reshape(F,[3,3,1,product(grid(1:2))*grid3])
endif endif
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
@ -606,7 +606,7 @@ subroutine formResidual(in, FandF_tau, &
do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1) do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1)
e = e + 1 e = e + 1
residual_F(1:3,1:3,i,j,k) = & residual_F(1:3,1:3,i,j,k) = &
math_mul3333xx33(math_invSym3333(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,e) + C_scale), & math_mul3333xx33(math_invSym3333(homogenization_dPdF(1:3,1:3,1:3,1:3,1,e) + C_scale), &
residual_F(1:3,1:3,i,j,k) - matmul(F(1:3,1:3,i,j,k), & residual_F(1:3,1:3,i,j,k) - matmul(F(1:3,1:3,i,j,k), &
math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))) & math_mul3333xx33(C_scale,F_tau(1:3,1:3,i,j,k) - F(1:3,1:3,i,j,k) - math_I3))) &
+ residual_F_tau(1:3,1:3,i,j,k) + residual_F_tau(1:3,1:3,i,j,k)

20
src/grid/spectral_utilities.f90
View File

@ -802,7 +802,7 @@ end subroutine utilities_fourierTensorDivergence
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief calculate constitutive response from materialpoint_F0 to F during timeinc !> @brief calculate constitutive response from homogenization_F0 to F during timeinc
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,& subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
F,timeinc,rotation_BC) F,timeinc,rotation_BC)
@ -824,11 +824,11 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
print'(/,a)', ' ... evaluating constitutive response ......................................' print'(/,a)', ' ... evaluating constitutive response ......................................'
flush(IO_STDOUT) flush(IO_STDOUT)
materialpoint_F = reshape(F,[3,3,1,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field homogenization_F = reshape(F,[3,3,1,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field
call materialpoint_stressAndItsTangent(timeinc) ! calculate P field call materialpoint_stressAndItsTangent(timeinc) ! calculate P field
P = reshape(materialpoint_P, [3,3,grid(1),grid(2),grid3]) P = reshape(homogenization_P, [3,3,grid(1),grid(2),grid3])
P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt ! average of P P_av = sum(sum(sum(P,dim=5),dim=4),dim=3) * wgt ! average of P
call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
if (debugRotation) & if (debugRotation) &
@ -845,13 +845,13 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
dPdF_min = huge(1.0_pReal) dPdF_min = huge(1.0_pReal)
dPdF_norm_min = huge(1.0_pReal) dPdF_norm_min = huge(1.0_pReal)
do i = 1, product(grid(1:2))*grid3 do i = 1, product(grid(1:2))*grid3
if (dPdF_norm_max < sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal)) then if (dPdF_norm_max < sum(homogenization_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal)) then
dPdF_max = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,i) dPdF_max = homogenization_dPdF(1:3,1:3,1:3,1:3,1,i)
dPdF_norm_max = sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal) dPdF_norm_max = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal)
endif endif
if (dPdF_norm_min > sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal)) then if (dPdF_norm_min > sum(homogenization_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal)) then
dPdF_min = materialpoint_dPdF(1:3,1:3,1:3,1:3,1,i) dPdF_min = homogenization_dPdF(1:3,1:3,1:3,1:3,1,i)
dPdF_norm_min = sum(materialpoint_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal) dPdF_norm_min = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal)
endif endif
end do end do
@ -869,7 +869,7 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
C_minmaxAvg = 0.5_pReal*(dPdF_max + dPdF_min) C_minmaxAvg = 0.5_pReal*(dPdF_max + dPdF_min)
C_volAvg = sum(sum(materialpoint_dPdF,dim=6),dim=5) C_volAvg = sum(sum(homogenization_dPdF,dim=6),dim=5)
call MPI_Allreduce(MPI_IN_PLACE,C_volAvg,81,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,C_volAvg,81,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
C_volAvg = C_volAvg * wgt C_volAvg = C_volAvg * wgt

76
src/homogenization.f90
View File

@ -31,12 +31,12 @@ module homogenization
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! General variables for the homogenization at a material point ! General variables for the homogenization at a material point
real(pReal), dimension(:,:,:,:), allocatable, public :: & real(pReal), dimension(:,:,:,:), allocatable, public :: &
materialpoint_F0, & !< def grad of IP at start of FE increment homogenization_F0, & !< def grad of IP at start of FE increment
materialpoint_F !< def grad of IP to be reached at end of FE increment homogenization_F !< def grad of IP to be reached at end of FE increment
real(pReal), dimension(:,:,:,:), allocatable, public, protected :: & real(pReal), dimension(:,:,:,:), allocatable, public, protected :: &
materialpoint_P !< first P--K stress of IP homogenization_P !< first P--K stress of IP
real(pReal), dimension(:,:,:,:,:,:), allocatable, public, protected :: & real(pReal), dimension(:,:,:,:,:,:), allocatable, public, protected :: &
materialpoint_dPdF !< tangent of first P--K stress at IP homogenization_dPdF !< tangent of first P--K stress at IP
type :: tNumerics type :: tNumerics
integer :: & integer :: &
@ -158,7 +158,7 @@ subroutine homogenization_init
debugHomog%grain = config_debug%get_asInt('grain',defaultVal = 1) debugHomog%grain = config_debug%get_asInt('grain',defaultVal = 1)
if (debugHomog%grain < 1 & if (debugHomog%grain < 1 &
.or. debugHomog%grain > homogenization_Ngrains(material_homogenizationAt(debugHomog%element))) & .or. debugHomog%grain > homogenization_Nconstituents(material_homogenizationAt(debugHomog%element))) &
call IO_error(602,ext_msg='constituent', el=debugHomog%element, g=debugHomog%grain) call IO_error(602,ext_msg='constituent', el=debugHomog%element, g=debugHomog%grain)
@ -181,10 +181,10 @@ subroutine homogenization_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! allocate and initialize global variables ! allocate and initialize global variables
allocate(materialpoint_dPdF(3,3,3,3,discretization_nIP,discretization_nElem), source=0.0_pReal) allocate(homogenization_dPdF(3,3,3,3,discretization_nIPs,discretization_Nelems), source=0.0_pReal)
materialpoint_F0 = spread(spread(math_I3,3,discretization_nIP),4,discretization_nElem) ! initialize to identity homogenization_F0 = spread(spread(math_I3,3,discretization_nIPs),4,discretization_Nelems) ! initialize to identity
materialpoint_F = materialpoint_F0 ! initialize to identity homogenization_F = homogenization_F0 ! initialize to identity
allocate(materialpoint_P(3,3,discretization_nIP,discretization_nElem), source=0.0_pReal) allocate(homogenization_P(3,3,discretization_nIPs,discretization_Nelems), source=0.0_pReal)
print'(/,a)', ' <<<+- homogenization init -+>>>'; flush(IO_STDOUT) print'(/,a)', ' <<<+- homogenization init -+>>>'; flush(IO_STDOUT)
@ -213,13 +213,13 @@ subroutine materialpoint_stressAndItsTangent(dt)
i, & !< integration point number i, & !< integration point number
e, & !< element number e, & !< element number
myNgrains myNgrains
real(pReal), dimension(discretization_nIP,discretization_nElem) :: & real(pReal), dimension(discretization_nIPs,discretization_Nelems) :: &
subFrac, & subFrac, &
subStep subStep
logical, dimension(discretization_nIP,discretization_nElem) :: & logical, dimension(discretization_nIPs,discretization_Nelems) :: &
requested, & requested, &
converged converged
logical, dimension(2,discretization_nIP,discretization_nElem) :: & logical, dimension(2,discretization_nIPs,discretization_Nelems) :: &
doneAndHappy doneAndHappy
@ -257,7 +257,7 @@ subroutine materialpoint_stressAndItsTangent(dt)
!$OMP PARALLEL DO !$OMP PARALLEL DO
elementLooping1: do e = FEsolving_execElem(1),FEsolving_execElem(2) elementLooping1: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(material_homogenizationAt(e)) myNgrains = homogenization_Nconstituents(material_homogenizationAt(e))
IpLooping1: do i = FEsolving_execIP(1),FEsolving_execIP(2) IpLooping1: do i = FEsolving_execIP(1),FEsolving_execIP(2)
if (converged(i,e)) then if (converged(i,e)) then
@ -327,11 +327,11 @@ subroutine materialpoint_stressAndItsTangent(dt)
! deformation partitioning ! deformation partitioning
!$OMP PARALLEL DO PRIVATE(myNgrains) !$OMP PARALLEL DO PRIVATE(myNgrains)
elementLooping2: do e = FEsolving_execElem(1),FEsolving_execElem(2) elementLooping2: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Ngrains(material_homogenizationAt(e)) myNgrains = homogenization_Nconstituents(material_homogenizationAt(e))
IpLooping2: do i = FEsolving_execIP(1),FEsolving_execIP(2) IpLooping2: do i = FEsolving_execIP(1),FEsolving_execIP(2)
if(requested(i,e) .and. .not. doneAndHappy(1,i,e)) then ! requested but not yet done if(requested(i,e) .and. .not. doneAndHappy(1,i,e)) then ! requested but not yet done
call partitionDeformation(materialpoint_F0(1:3,1:3,i,e) & call partitionDeformation(homogenization_F0(1:3,1:3,i,e) &
+ (materialpoint_F(1:3,1:3,i,e)-materialpoint_F0(1:3,1:3,i,e))& + (homogenization_F(1:3,1:3,i,e)-homogenization_F0(1:3,1:3,i,e))&
*(subStep(i,e)+subFrac(i,e)), & *(subStep(i,e)+subFrac(i,e)), &
i,e) i,e)
crystallite_dt(1:myNgrains,i,e) = dt*subStep(i,e) ! propagate materialpoint dt to grains crystallite_dt(1:myNgrains,i,e) = dt*subStep(i,e) ! propagate materialpoint dt to grains
@ -357,8 +357,8 @@ subroutine materialpoint_stressAndItsTangent(dt)
doneAndHappy(1:2,i,e) = [.true.,.false.] doneAndHappy(1:2,i,e) = [.true.,.false.]
else else
doneAndHappy(1:2,i,e) = updateState(dt*subStep(i,e), & doneAndHappy(1:2,i,e) = updateState(dt*subStep(i,e), &
materialpoint_F0(1:3,1:3,i,e) & homogenization_F0(1:3,1:3,i,e) &
+ (materialpoint_F(1:3,1:3,i,e)-materialpoint_F0(1:3,1:3,i,e)) & + (homogenization_F(1:3,1:3,i,e)-homogenization_F0(1:3,1:3,i,e)) &
*(subStep(i,e)+subFrac(i,e)), & *(subStep(i,e)+subFrac(i,e)), &
i,e) i,e)
converged(i,e) = all(doneAndHappy(1:2,i,e)) ! converged if done and happy converged(i,e) = all(doneAndHappy(1:2,i,e)) ! converged if done and happy
@ -408,12 +408,12 @@ subroutine partitionDeformation(subF,ip,el)
case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
call mech_isostrain_partitionDeformation(& call mech_isostrain_partitionDeformation(&
crystallite_partitionedF(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), & crystallite_partitionedF(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
subF) subF)
case (HOMOGENIZATION_RGC_ID) chosenHomogenization case (HOMOGENIZATION_RGC_ID) chosenHomogenization
call mech_RGC_partitionDeformation(& call mech_RGC_partitionDeformation(&
crystallite_partitionedF(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), & crystallite_partitionedF(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
subF,& subF,&
ip, & ip, &
el) el)
@ -437,19 +437,19 @@ function updateState(subdt,subF,ip,el)
el !< element number el !< element number
integer :: c integer :: c
logical, dimension(2) :: updateState logical, dimension(2) :: updateState
real(pReal) :: dPdFs(3,3,3,3,homogenization_Ngrains(material_homogenizationAt(el))) real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
updateState = .true. updateState = .true.
chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el))) chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
case (HOMOGENIZATION_RGC_ID) chosenHomogenization case (HOMOGENIZATION_RGC_ID) chosenHomogenization
do c=1,homogenization_Ngrains(material_homogenizationAt(el)) do c=1,homogenization_Nconstituents(material_homogenizationAt(el))
dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el) dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el)
enddo enddo
updateState = & updateState = &
updateState .and. & updateState .and. &
mech_RGC_updateState(crystallite_P(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), & 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_Ngrains(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_Ngrains(material_homogenizationAt(el)),ip,el),& crystallite_partitionedF0(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el),&
subF,& subF,&
subdt, & subdt, &
dPdFs, & dPdFs, &
@ -487,33 +487,33 @@ subroutine averageStressAndItsTangent(ip,el)
ip, & !< integration point ip, & !< integration point
el !< element number el !< element number
integer :: c integer :: c
real(pReal) :: dPdFs(3,3,3,3,homogenization_Ngrains(material_homogenizationAt(el))) real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el))) chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
case (HOMOGENIZATION_NONE_ID) chosenHomogenization case (HOMOGENIZATION_NONE_ID) chosenHomogenization
materialpoint_P(1:3,1:3,ip,el) = crystallite_P(1:3,1:3,1,ip,el) homogenization_P(1:3,1:3,ip,el) = crystallite_P(1:3,1:3,1,ip,el)
materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el) = crystallite_stressTangent(1,ip,el) homogenization_dPdF(1:3,1:3,1:3,1:3,ip,el) = crystallite_stressTangent(1,ip,el)
case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
do c = 1, homogenization_Ngrains(material_homogenizationAt(el)) do c = 1, homogenization_Nconstituents(material_homogenizationAt(el))
dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el) dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el)
enddo enddo
call mech_isostrain_averageStressAndItsTangent(& call mech_isostrain_averageStressAndItsTangent(&
materialpoint_P(1:3,1:3,ip,el), & homogenization_P(1:3,1:3,ip,el), &
materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el),& homogenization_dPdF(1:3,1:3,1:3,1:3,ip,el),&
crystallite_P(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), & crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
dPdFs, & dPdFs, &
homogenization_typeInstance(material_homogenizationAt(el))) homogenization_typeInstance(material_homogenizationAt(el)))
case (HOMOGENIZATION_RGC_ID) chosenHomogenization case (HOMOGENIZATION_RGC_ID) chosenHomogenization
do c = 1, homogenization_Ngrains(material_homogenizationAt(el)) do c = 1, homogenization_Nconstituents(material_homogenizationAt(el))
dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el) dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el)
enddo enddo
call mech_RGC_averageStressAndItsTangent(& call mech_RGC_averageStressAndItsTangent(&
materialpoint_P(1:3,1:3,ip,el), & homogenization_P(1:3,1:3,ip,el), &
materialpoint_dPdF(1:3,1:3,1:3,1:3,ip,el),& homogenization_dPdF(1:3,1:3,1:3,1:3,ip,el),&
crystallite_P(1:3,1:3,1:homogenization_Ngrains(material_homogenizationAt(el)),ip,el), & crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
dPdFs, & dPdFs, &
homogenization_typeInstance(material_homogenizationAt(el))) homogenization_typeInstance(material_homogenizationAt(el)))
end select chosenHomogenization end select chosenHomogenization
@ -539,10 +539,10 @@ subroutine homogenization_results
group = trim(group_base)//'/generic' group = trim(group_base)//'/generic'
call results_closeGroup(results_addGroup(group)) call results_closeGroup(results_addGroup(group))
!temp = reshape(materialpoint_F,[3,3,discretization_nIP*discretization_nElem]) !temp = reshape(homogenization_F,[3,3,discretization_nIPs*discretization_Nelems])
!call results_writeDataset(group,temp,'F',& !call results_writeDataset(group,temp,'F',&
! 'deformation gradient','1') ! 'deformation gradient','1')
!temp = reshape(materialpoint_P,[3,3,discretization_nIP*discretization_nElem]) !temp = reshape(homogenization_P,[3,3,discretization_nIPs*discretization_Nelems])
!call results_writeDataset(group,temp,'P',& !call results_writeDataset(group,temp,'P',&
! '1st Piola-Kirchhoff stress','Pa') ! '1st Piola-Kirchhoff stress','Pa')

38
src/homogenization_mech_RGC.f90
View File

@ -81,9 +81,9 @@ module subroutine mech_RGC_init(num_homogMech)
num_homogMech !< pointer to mechanical homogenization numerics data num_homogMech !< pointer to mechanical homogenization numerics data
integer :: & integer :: &
Ninstance, & Ninstances, &
h, & h, &
NofMyHomog, & Nmaterialpoints, &
sizeState, nIntFaceTot sizeState, nIntFaceTot
class (tNode), pointer :: & class (tNode), pointer :: &
@ -94,8 +94,8 @@ module subroutine mech_RGC_init(num_homogMech)
print'(/,a)', ' <<<+- homogenization_mech_rgc init -+>>>' print'(/,a)', ' <<<+- homogenization_mech_rgc init -+>>>'
Ninstance = count(homogenization_type == HOMOGENIZATION_RGC_ID) Ninstances = count(homogenization_type == HOMOGENIZATION_RGC_ID)
print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
print*, 'Tjahjanto et al., International Journal of Material Forming 2(1):939–942, 2009' print*, 'Tjahjanto et al., International Journal of Material Forming 2(1):939–942, 2009'
print*, 'https://doi.org/10.1007/s12289-009-0619-1'//IO_EOL print*, 'https://doi.org/10.1007/s12289-009-0619-1'//IO_EOL
@ -105,10 +105,10 @@ module subroutine mech_RGC_init(num_homogMech)
allocate(param(Ninstance)) allocate(param(Ninstances))
allocate(state(Ninstance)) allocate(state(Ninstances))
allocate(state0(Ninstance)) allocate(state0(Ninstances))
allocate(dependentState(Ninstance)) allocate(dependentState(Ninstances))
num_RGC => num_homogMech%get('RGC',defaultVal=emptyDict) num_RGC => num_homogMech%get('RGC',defaultVal=emptyDict)
@ -164,7 +164,7 @@ module subroutine mech_RGC_init(num_homogMech)
#endif #endif
prm%N_constituents = homogMech%get_asInts('cluster_size',requiredSize=3) prm%N_constituents = homogMech%get_asInts('cluster_size',requiredSize=3)
if (homogenization_Ngrains(h) /= product(prm%N_constituents)) & if (homogenization_Nconstituents(h) /= product(prm%N_constituents)) &
call IO_error(211,ext_msg='N_constituents (mech_RGC)') call IO_error(211,ext_msg='N_constituents (mech_RGC)')
prm%xi_alpha = homogMech%get_asFloat('xi_alpha') prm%xi_alpha = homogMech%get_asFloat('xi_alpha')
@ -173,7 +173,7 @@ module subroutine mech_RGC_init(num_homogMech)
prm%D_alpha = homogMech%get_asFloats('D_alpha', requiredSize=3) prm%D_alpha = homogMech%get_asFloats('D_alpha', requiredSize=3)
prm%a_g = homogMech%get_asFloats('a_g', requiredSize=3) prm%a_g = homogMech%get_asFloats('a_g', requiredSize=3)
NofMyHomog = count(material_homogenizationAt == h) Nmaterialpoints = count(material_homogenizationAt == h)
nIntFaceTot = 3*( (prm%N_constituents(1)-1)*prm%N_constituents(2)*prm%N_constituents(3) & 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)-1)*prm%N_constituents(3) &
+ prm%N_constituents(1)*prm%N_constituents(2)*(prm%N_constituents(3)-1)) + prm%N_constituents(1)*prm%N_constituents(2)*(prm%N_constituents(3)-1))
@ -181,24 +181,24 @@ module subroutine mech_RGC_init(num_homogMech)
+ size(['avg constitutive work ','average penalty energy']) + size(['avg constitutive work ','average penalty energy'])
homogState(h)%sizeState = sizeState homogState(h)%sizeState = sizeState
allocate(homogState(h)%state0 (sizeState,NofMyHomog), source=0.0_pReal) allocate(homogState(h)%state0 (sizeState,Nmaterialpoints), source=0.0_pReal)
allocate(homogState(h)%subState0(sizeState,NofMyHomog), source=0.0_pReal) allocate(homogState(h)%subState0(sizeState,Nmaterialpoints), source=0.0_pReal)
allocate(homogState(h)%state (sizeState,NofMyHomog), source=0.0_pReal) allocate(homogState(h)%state (sizeState,Nmaterialpoints), source=0.0_pReal)
stt%relaxationVector => homogState(h)%state(1:nIntFaceTot,:) stt%relaxationVector => homogState(h)%state(1:nIntFaceTot,:)
st0%relaxationVector => homogState(h)%state0(1:nIntFaceTot,:) st0%relaxationVector => homogState(h)%state0(1:nIntFaceTot,:)
stt%work => homogState(h)%state(nIntFaceTot+1,:) stt%work => homogState(h)%state(nIntFaceTot+1,:)
stt%penaltyEnergy => homogState(h)%state(nIntFaceTot+2,:) stt%penaltyEnergy => homogState(h)%state(nIntFaceTot+2,:)
allocate(dst%volumeDiscrepancy( NofMyHomog), source=0.0_pReal) allocate(dst%volumeDiscrepancy( Nmaterialpoints), source=0.0_pReal)
allocate(dst%relaxationRate_avg( NofMyHomog), source=0.0_pReal) allocate(dst%relaxationRate_avg( Nmaterialpoints), source=0.0_pReal)
allocate(dst%relaxationRate_max( NofMyHomog), source=0.0_pReal) allocate(dst%relaxationRate_max( Nmaterialpoints), source=0.0_pReal)
allocate(dst%mismatch( 3,NofMyHomog), source=0.0_pReal) allocate(dst%mismatch( 3,Nmaterialpoints), source=0.0_pReal)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! assigning cluster orientations ! assigning cluster orientations
dependentState(homogenization_typeInstance(h))%orientation = spread(eu2om(prm%a_g*inRad),3,NofMyHomog) dependentState(homogenization_typeInstance(h))%orientation = spread(eu2om(prm%a_g*inRad),3,Nmaterialpoints)
!dst%orientation = spread(eu2om(prm%a_g*inRad),3,NofMyHomog) ifort version 18.0.1 crashes (for whatever reason) !dst%orientation = spread(eu2om(prm%a_g*inRad),3,Nmaterialpoints) ifort version 18.0.1 crashes (for whatever reason)
end associate end associate

22
src/homogenization_mech_isostrain.f90
View File

@ -18,7 +18,7 @@ submodule(homogenization) homogenization_mech_isostrain
mapping mapping
end type end type
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances)
contains contains
@ -29,9 +29,9 @@ contains
module subroutine mech_isostrain_init module subroutine mech_isostrain_init
integer :: & integer :: &
Ninstance, & Ninstances, &
h, & h, &
NofMyHomog Nmaterialpoints
class(tNode), pointer :: & class(tNode), pointer :: &
material_homogenization, & material_homogenization, &
homog, & homog, &
@ -39,10 +39,10 @@ module subroutine mech_isostrain_init
print'(/,a)', ' <<<+- homogenization_mech_isostrain init -+>>>' print'(/,a)', ' <<<+- homogenization_mech_isostrain init -+>>>'
Ninstance = count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID) Ninstances = count(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)
print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
allocate(param(Ninstance)) ! one container of parameters per instance allocate(param(Ninstances)) ! one container of parameters per instance
material_homogenization => config_material%get('homogenization') material_homogenization => config_material%get('homogenization')
do h = 1, size(homogenization_type) do h = 1, size(homogenization_type)
@ -51,7 +51,7 @@ module subroutine mech_isostrain_init
homogMech => homog%get('mech') homogMech => homog%get('mech')
associate(prm => param(homogenization_typeInstance(h))) associate(prm => param(homogenization_typeInstance(h)))
prm%N_constituents = homogenization_Ngrains(h) prm%N_constituents = homogenization_Nconstituents(h)
select case(homogMech%get_asString('mapping',defaultVal = 'sum')) select case(homogMech%get_asString('mapping',defaultVal = 'sum'))
case ('sum') case ('sum')
prm%mapping = parallel_ID prm%mapping = parallel_ID
@ -61,11 +61,11 @@ module subroutine mech_isostrain_init
call IO_error(211,ext_msg='sum'//' (mech_isostrain)') call IO_error(211,ext_msg='sum'//' (mech_isostrain)')
end select end select
NofMyHomog = count(material_homogenizationAt == h) Nmaterialpoints = count(material_homogenizationAt == h)
homogState(h)%sizeState = 0 homogState(h)%sizeState = 0
allocate(homogState(h)%state0 (0,NofMyHomog)) allocate(homogState(h)%state0 (0,Nmaterialpoints))
allocate(homogState(h)%subState0(0,NofMyHomog)) allocate(homogState(h)%subState0(0,Nmaterialpoints))
allocate(homogState(h)%state (0,NofMyHomog)) allocate(homogState(h)%state (0,Nmaterialpoints))
end associate end associate

18
src/homogenization_mech_none.f90
View File

@ -14,26 +14,26 @@ contains
module subroutine mech_none_init module subroutine mech_none_init
integer :: & integer :: &
Ninstance, & Ninstances, &
h, & h, &
NofMyHomog Nmaterialpoints
print'(/,a)', ' <<<+- homogenization_mech_none init -+>>>' print'(/,a)', ' <<<+- homogenization_mech_none init -+>>>'
Ninstance = count(homogenization_type == HOMOGENIZATION_NONE_ID) Ninstances = count(homogenization_type == HOMOGENIZATION_NONE_ID)
print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
do h = 1, size(homogenization_type) do h = 1, size(homogenization_type)
if(homogenization_type(h) /= HOMOGENIZATION_NONE_ID) cycle if(homogenization_type(h) /= HOMOGENIZATION_NONE_ID) cycle
if(homogenization_Ngrains(h) /= 1) & if(homogenization_Nconstituents(h) /= 1) &
call IO_error(211,ext_msg='N_constituents (mech_none)') call IO_error(211,ext_msg='N_constituents (mech_none)')
NofMyHomog = count(material_homogenizationAt == h) Nmaterialpoints = count(material_homogenizationAt == h)
homogState(h)%sizeState = 0 homogState(h)%sizeState = 0
allocate(homogState(h)%state0 (0,NofMyHomog)) allocate(homogState(h)%state0 (0,Nmaterialpoints))
allocate(homogState(h)%subState0(0,NofMyHomog)) allocate(homogState(h)%subState0(0,Nmaterialpoints))
allocate(homogState(h)%state (0,NofMyHomog)) allocate(homogState(h)%state (0,Nmaterialpoints))
enddo enddo

12
src/kinematics_cleavage_opening.f90
View File

@ -20,7 +20,7 @@ submodule(constitutive:constitutive_damage) kinematics_cleavage_opening
cleavage_systems cleavage_systems
end type tParameters end type tParameters
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances)
contains contains
@ -35,7 +35,7 @@ module function kinematics_cleavage_opening_init(kinematics_length) result(myKin
integer, intent(in) :: kinematics_length integer, intent(in) :: kinematics_length
logical, dimension(:,:), allocatable :: myKinematics logical, dimension(:,:), allocatable :: myKinematics
integer :: Ninstance,p,k integer :: Ninstances,p,k
integer, dimension(:), allocatable :: N_cl !< active number of cleavage systems per family integer, dimension(:), allocatable :: N_cl !< active number of cleavage systems per family
character(len=pStringLen) :: extmsg = '' character(len=pStringLen) :: extmsg = ''
class(tNode), pointer :: & class(tNode), pointer :: &
@ -48,12 +48,12 @@ module function kinematics_cleavage_opening_init(kinematics_length) result(myKin
print'(/,a)', ' <<<+- kinematics_cleavage_opening init -+>>>' print'(/,a)', ' <<<+- kinematics_cleavage_opening init -+>>>'
myKinematics = kinematics_active('cleavage_opening',kinematics_length) myKinematics = kinematics_active('cleavage_opening',kinematics_length)
Ninstance = count(myKinematics) Ninstances = count(myKinematics)
print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstances == 0) return
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(param(Ninstance)) allocate(param(Ninstances))
allocate(kinematics_cleavage_opening_instance(phases%length), source=0) allocate(kinematics_cleavage_opening_instance(phases%length), source=0)
do p = 1, phases%length do p = 1, phases%length

12
src/kinematics_slipplane_opening.f90
View File

@ -22,7 +22,7 @@ submodule(constitutive:constitutive_damage) kinematics_slipplane_opening
P_n P_n
end type tParameters end type tParameters
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances)
contains contains
@ -37,7 +37,7 @@ module function kinematics_slipplane_opening_init(kinematics_length) result(myKi
integer, intent(in) :: kinematics_length integer, intent(in) :: kinematics_length
logical, dimension(:,:), allocatable :: myKinematics logical, dimension(:,:), allocatable :: myKinematics
integer :: Ninstance,p,i,k integer :: Ninstances,p,i,k
character(len=pStringLen) :: extmsg = '' character(len=pStringLen) :: extmsg = ''
integer, dimension(:), allocatable :: N_sl integer, dimension(:), allocatable :: N_sl
real(pReal), dimension(:,:), allocatable :: d,n,t real(pReal), dimension(:,:), allocatable :: d,n,t
@ -51,13 +51,13 @@ module function kinematics_slipplane_opening_init(kinematics_length) result(myKi
print'(/,a)', ' <<<+- kinematics_slipplane init -+>>>' print'(/,a)', ' <<<+- kinematics_slipplane init -+>>>'
myKinematics = kinematics_active('slipplane_opening',kinematics_length) myKinematics = kinematics_active('slipplane_opening',kinematics_length)
Ninstance = count(myKinematics) Ninstances = count(myKinematics)
print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstances == 0) return
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(kinematics_slipplane_opening_instance(phases%length), source=0) allocate(kinematics_slipplane_opening_instance(phases%length), source=0)
allocate(param(Ninstance)) allocate(param(Ninstances))
do p = 1, phases%length do p = 1, phases%length
if(any(myKinematics(:,p))) kinematics_slipplane_opening_instance(p) = count(myKinematics(:,1:p)) if(any(myKinematics(:,p))) kinematics_slipplane_opening_instance(p) = count(myKinematics(:,1:p))

10
src/kinematics_thermal_expansion.f90
View File

@ -29,7 +29,7 @@ module function kinematics_thermal_expansion_init(kinematics_length) result(myKi
integer, intent(in) :: kinematics_length integer, intent(in) :: kinematics_length
logical, dimension(:,:), allocatable :: myKinematics logical, dimension(:,:), allocatable :: myKinematics
integer :: Ninstance,p,i,k integer :: Ninstances,p,i,k
real(pReal), dimension(:), allocatable :: temp real(pReal), dimension(:), allocatable :: temp
class(tNode), pointer :: & class(tNode), pointer :: &
phases, & phases, &
@ -41,12 +41,12 @@ module function kinematics_thermal_expansion_init(kinematics_length) result(myKi
print'(/,a)', ' <<<+- kinematics_thermal_expansion init -+>>>' print'(/,a)', ' <<<+- kinematics_thermal_expansion init -+>>>'
myKinematics = kinematics_active('thermal_expansion',kinematics_length) myKinematics = kinematics_active('thermal_expansion',kinematics_length)
Ninstance = count(myKinematics) Ninstances = count(myKinematics)
print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstances == 0) return
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(param(Ninstance)) allocate(param(Ninstances))
allocate(kinematics_thermal_expansion_instance(phases%length), source=0) allocate(kinematics_thermal_expansion_instance(phases%length), source=0)
do p = 1, phases%length do p = 1, phases%length

195
src/material.f90
View File

@ -52,7 +52,7 @@ module material
HOMOGENIZATION_RGC_ID HOMOGENIZATION_RGC_ID
end enum end enum
character(len=pStringLen), public, protected, allocatable, dimension(:) :: & character(len=pStringLen), public, protected, allocatable, dimension(:) :: &
material_name_phase, & !< name of each phase material_name_phase, & !< name of each phase
material_name_homogenization !< name of each homogenization material_name_homogenization !< name of each homogenization
@ -64,13 +64,10 @@ module material
homogenization_type !< type of each homogenization homogenization_type !< type of each homogenization
integer, public, protected :: & integer, public, protected :: &
material_Nhomogenization !< number of homogenizations homogenization_maxNconstituents !< max number of grains in any USED homogenization
integer, public, protected :: &
homogenization_maxNgrains !< max number of grains in any USED homogenization
integer, dimension(:), allocatable, public, protected :: & integer, dimension(:), allocatable, public, protected :: &
homogenization_Ngrains, & !< number of grains in each homogenization homogenization_Nconstituents, & !< number of grains in each homogenization
homogenization_typeInstance, & !< instance of particular type of each homogenization homogenization_typeInstance, & !< instance of particular type of each homogenization
thermal_typeInstance, & !< instance of particular type of each thermal transport thermal_typeInstance, & !< instance of particular type of each thermal transport
damage_typeInstance !< instance of particular type of each nonlocal damage damage_typeInstance !< instance of particular type of each nonlocal damage
@ -83,9 +80,9 @@ module material
material_homogenizationAt !< homogenization ID of each element material_homogenizationAt !< homogenization ID of each element
integer, dimension(:,:), allocatable, public, target :: & ! (ip,elem) ToDo: ugly target for mapping hack integer, dimension(:,:), allocatable, public, target :: & ! (ip,elem) ToDo: ugly target for mapping hack
material_homogenizationMemberAt !< position of the element within its homogenization instance material_homogenizationMemberAt !< position of the element within its homogenization instance
integer, dimension(:,:), allocatable, public, protected :: & ! (constituent,elem) integer, dimension(:,:), allocatable, public, protected :: & ! (constituent,elem)
material_phaseAt !< phase ID of each element material_phaseAt !< phase ID of each element
integer, dimension(:,:,:), allocatable, public, protected :: & ! (constituent,elem) integer, dimension(:,:,:), allocatable, public, protected :: & ! (constituent,IP,elem)
material_phaseMemberAt !< position of the element within its phase instance material_phaseMemberAt !< position of the element within its phase instance
type(tState), allocatable, dimension(:), public :: & type(tState), allocatable, dimension(:), public :: &
@ -96,11 +93,6 @@ module material
type(Rotation), dimension(:,:,:), allocatable, public, protected :: & type(Rotation), dimension(:,:,:), allocatable, public, protected :: &
material_orientation0 !< initial orientation of each grain,IP,element material_orientation0 !< initial orientation of each grain,IP,element
integer, dimension(:), allocatable, private :: &
material_Nconstituents !< number of constituents in each material
! BEGIN DEPRECATED ! BEGIN DEPRECATED
integer, dimension(:,:), allocatable, private, target :: mappingHomogenizationConst !< mapping from material points to offset in constant state/field integer, dimension(:,:), allocatable, private, target :: mappingHomogenizationConst !< mapping from material points to offset in constant state/field
! END DEPRECATED ! END DEPRECATED
@ -157,28 +149,10 @@ contains
subroutine material_init(restart) subroutine material_init(restart)
logical, intent(in) :: restart logical, intent(in) :: restart
integer :: myHomog
integer :: ph, myHomog
class(tNode), pointer :: &
phases, &
material_homogenization
character(len=pStringLen) :: sectionName
print'(/,a)', ' <<<+- material init -+>>>'; flush(IO_STDOUT) print'(/,a)', ' <<<+- material init -+>>>'; flush(IO_STDOUT)
phases => config_material%get('phase')
allocate(material_name_phase(phases%length))
do ph = 1, phases%length
write(sectionName,'(i0,a)') ph,'_'
material_name_phase(ph) = trim(adjustl(sectionName))//phases%getKey(ph) !ToDO: No reason to do. Update damage tests
enddo
material_homogenization => config_material%get('homogenization')
allocate(material_name_homogenization(material_homogenization%length))
do myHomog = 1, material_homogenization%length
write(sectionName,'(i0,a)') myHomog,'_'
material_name_homogenization(myHomog) = trim(adjustl(sectionName))//material_homogenization%getKey(myHomog)
enddo
call material_parseMaterial call material_parseMaterial
print*, 'Material parsed' print*, 'Material parsed'
@ -187,34 +161,32 @@ subroutine material_init(restart)
print*, 'Homogenization parsed' print*, 'Homogenization parsed'
if(homogenization_maxNgrains > size(material_phaseAt,1)) call IO_error(148) allocate(homogState (size(material_name_homogenization)))
allocate(thermalState (size(material_name_homogenization)))
allocate(damageState (size(material_name_homogenization)))
allocate(homogState (material_Nhomogenization)) allocate(thermalMapping (size(material_name_homogenization)))
allocate(thermalState (material_Nhomogenization)) allocate(damageMapping (size(material_name_homogenization)))
allocate(damageState (material_Nhomogenization))
allocate(thermalMapping (material_Nhomogenization)) allocate(temperature (size(material_name_homogenization)))
allocate(damageMapping (material_Nhomogenization)) allocate(damage (size(material_name_homogenization)))
allocate(temperature (material_Nhomogenization)) allocate(temperatureRate (size(material_name_homogenization)))
allocate(damage (material_Nhomogenization))
allocate(temperatureRate (material_Nhomogenization))
if (.not. restart) then if (.not. restart) then
call results_openJobFile call results_openJobFile
call results_mapping_constituent(material_phaseAt,material_phaseMemberAt,material_name_phase) call results_mapping_constituent(material_phaseAt,material_phaseMemberAt,material_name_phase)
call results_mapping_materialpoint(material_homogenizationAt,material_homogenizationMemberAt,material_name_homogenization) call results_mapping_homogenization(material_homogenizationAt,material_homogenizationMemberAt,material_name_homogenization)
call results_closeJobFile call results_closeJobFile
endif endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! BEGIN DEPRECATED ! BEGIN DEPRECATED
allocate(mappingHomogenizationConst( discretization_nIP,discretization_nElem),source=1) allocate(mappingHomogenizationConst( discretization_nIPs,discretization_Nelems),source=1)
! hack needed to initialize field values used during constitutive initialization ! hack needed to initialize field values used during constitutive initialization
do myHomog = 1,material_Nhomogenization do myHomog = 1, size(material_name_homogenization)
thermalMapping (myHomog)%p => mappingHomogenizationConst thermalMapping (myHomog)%p => mappingHomogenizationConst
damageMapping (myHomog)%p => mappingHomogenizationConst damageMapping (myHomog)%p => mappingHomogenizationConst
allocate(temperature (myHomog)%p(1), source=thermal_initialT(myHomog)) allocate(temperature (myHomog)%p(1), source=thermal_initialT(myHomog))
@ -225,6 +197,7 @@ subroutine material_init(restart)
end subroutine material_init end subroutine material_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief parses the homogenization part from the material configuration !> @brief parses the homogenization part from the material configuration
! ToDo: This should be done in homogenization ! ToDo: This should be done in homogenization
@ -241,22 +214,19 @@ subroutine material_parseHomogenization
integer :: h integer :: h
material_homogenization => config_material%get('homogenization') material_homogenization => config_material%get('homogenization')
material_Nhomogenization = material_homogenization%length
allocate(homogenization_type(material_Nhomogenization), source=HOMOGENIZATION_undefined_ID) allocate(homogenization_type(size(material_name_homogenization)), source=HOMOGENIZATION_undefined_ID)
allocate(thermal_type(material_Nhomogenization), source=THERMAL_isothermal_ID) allocate(thermal_type(size(material_name_homogenization)), source=THERMAL_isothermal_ID)
allocate(damage_type (material_Nhomogenization), source=DAMAGE_none_ID) allocate(damage_type (size(material_name_homogenization)), source=DAMAGE_none_ID)
allocate(homogenization_typeInstance(material_Nhomogenization), source=0) allocate(homogenization_typeInstance(size(material_name_homogenization)), source=0)
allocate(thermal_typeInstance(material_Nhomogenization), source=0) allocate(thermal_typeInstance(size(material_name_homogenization)), source=0)
allocate(damage_typeInstance(material_Nhomogenization), source=0) allocate(damage_typeInstance(size(material_name_homogenization)), source=0)
allocate(homogenization_Ngrains(material_Nhomogenization), source=0) allocate(thermal_initialT(size(material_name_homogenization)), source=300.0_pReal)
allocate(thermal_initialT(material_Nhomogenization), source=300.0_pReal) allocate(damage_initialPhi(size(material_name_homogenization)), source=1.0_pReal)
allocate(damage_initialPhi(material_Nhomogenization), source=1.0_pReal)
do h=1, material_Nhomogenization do h=1, size(material_name_homogenization)
homog => material_homogenization%get(h) homog => material_homogenization%get(h)
homogMech => homog%get('mech') homogMech => homog%get('mech')
homogenization_Ngrains(h) = homog%get_asInt('N_constituents')
select case (homogMech%get_asString('type')) select case (homogMech%get_asString('type'))
case('none') case('none')
homogenization_type(h) = HOMOGENIZATION_NONE_ID homogenization_type(h) = HOMOGENIZATION_NONE_ID
@ -302,15 +272,12 @@ subroutine material_parseHomogenization
endif endif
enddo enddo
do h=1, material_Nhomogenization do h=1, size(material_name_homogenization)
homogenization_typeInstance(h) = count(homogenization_type(1:h) == homogenization_type(h)) homogenization_typeInstance(h) = count(homogenization_type(1:h) == homogenization_type(h))
thermal_typeInstance(h) = count(thermal_type (1:h) == thermal_type (h)) thermal_typeInstance(h) = count(thermal_type (1:h) == thermal_type (h))
damage_typeInstance(h) = count(damage_type (1:h) == damage_type (h)) damage_typeInstance(h) = count(damage_type (1:h) == damage_type (h))
enddo enddo
homogenization_maxNgrains = maxval(homogenization_Ngrains)
end subroutine material_parseHomogenization end subroutine material_parseHomogenization
@ -324,7 +291,8 @@ subroutine material_parseMaterial
constituents, & !> list of constituents constituents, & !> list of constituents
constituent, & !> constituent definition constituent, & !> constituent definition
phases, & phases, &
homogenizations homogenizations, &
homogenization
integer, dimension(:), allocatable :: & integer, dimension(:), allocatable :: &
counterPhase, & counterPhase, &
@ -333,42 +301,40 @@ subroutine material_parseMaterial
real(pReal) :: & real(pReal) :: &
frac frac
integer :: & integer :: &
e, & e, i, c, &
i, & h
m, &
c, &
maxNconstituents
materials => config_material%get('material') materials => config_material%get('material')
if(any(discretization_materialAt > materials%length)) & phases => config_material%get('phase')
call IO_error(155,ext_msg='More materials requested than found in material.yaml')
allocate(material_Nconstituents(materials%length),source=0)
do m = 1, materials%length
material => materials%get(m)
constituents => material%get('constituents')
material_Nconstituents(m) = constituents%length
enddo
maxNconstituents = maxval(material_Nconstituents)
allocate(material_homogenizationAt(discretization_nElem),source=0)
allocate(material_homogenizationMemberAt(discretization_nIP,discretization_nElem),source=0)
allocate(material_phaseAt(maxNconstituents,discretization_nElem),source=0)
allocate(material_phaseMemberAt(maxNconstituents,discretization_nIP,discretization_nElem),source=0)
allocate(material_orientation0(maxNconstituents,discretization_nIP,discretization_nElem))
phases => config_material%get('phase')
allocate(counterPhase(phases%length),source=0)
homogenizations => config_material%get('homogenization') homogenizations => config_material%get('homogenization')
call sanityCheck(materials, homogenizations)
material_name_phase = getKeys(phases)
material_name_homogenization = getKeys(homogenizations)
allocate(homogenization_Nconstituents(homogenizations%length))
do h=1, homogenizations%length
homogenization => homogenizations%get(h)
homogenization_Nconstituents(h) = homogenization%get_asInt('N_constituents')
enddo
homogenization_maxNconstituents = maxval(homogenization_Nconstituents)
allocate(counterPhase(phases%length),source=0)
allocate(counterHomogenization(homogenizations%length),source=0) allocate(counterHomogenization(homogenizations%length),source=0)
do e = 1, discretization_nElem allocate(material_homogenizationAt(discretization_Nelems),source=0)
allocate(material_homogenizationMemberAt(discretization_nIPs,discretization_Nelems),source=0)
allocate(material_phaseAt(homogenization_maxNconstituents,discretization_Nelems),source=0)
allocate(material_phaseMemberAt(homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems),source=0)
allocate(material_orientation0(homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems))
do e = 1, discretization_Nelems
material => materials%get(discretization_materialAt(e)) material => materials%get(discretization_materialAt(e))
constituents => material%get('constituents') constituents => material%get('constituents')
material_homogenizationAt(e) = homogenizations%getIndex(material%get_asString('homogenization')) material_homogenizationAt(e) = homogenizations%getIndex(material%get_asString('homogenization'))
do i = 1, discretization_nIP do i = 1, discretization_nIPs
counterHomogenization(material_homogenizationAt(e)) = counterHomogenization(material_homogenizationAt(e)) + 1 counterHomogenization(material_homogenizationAt(e)) = counterHomogenization(material_homogenizationAt(e)) + 1
material_homogenizationMemberAt(i,e) = counterHomogenization(material_homogenizationAt(e)) material_homogenizationMemberAt(i,e) = counterHomogenization(material_homogenizationAt(e))
enddo enddo
@ -379,11 +345,11 @@ subroutine material_parseMaterial
frac = frac + constituent%get_asFloat('fraction') frac = frac + constituent%get_asFloat('fraction')
material_phaseAt(c,e) = phases%getIndex(constituent%get_asString('phase')) material_phaseAt(c,e) = phases%getIndex(constituent%get_asString('phase'))
do i = 1, discretization_nIP do i = 1, discretization_nIPs
counterPhase(material_phaseAt(c,e)) = counterPhase(material_phaseAt(c,e)) + 1 counterPhase(material_phaseAt(c,e)) = counterPhase(material_phaseAt(c,e)) + 1
material_phaseMemberAt(c,i,e) = counterPhase(material_phaseAt(c,e)) material_phaseMemberAt(c,i,e) = counterPhase(material_phaseAt(c,e))
call material_orientation0(c,i,e)%fromQuaternion(constituent%get_asFloats('O',requiredSize=4)) call material_orientation0(c,i,e)%fromQuaternion(constituent%get_asFloats('O',requiredSize=4)) ! should be done in crystallite
enddo enddo
enddo enddo
@ -394,4 +360,49 @@ subroutine material_parseMaterial
end subroutine material_parseMaterial end subroutine material_parseMaterial
!--------------------------------------------------------------------------------------------------
!> @brief Check if material.yaml is consistent and contains sufficient # of materials
!--------------------------------------------------------------------------------------------------
subroutine sanityCheck(materials,homogenizations)
class(tNode), intent(in) :: materials, &
homogenizations
class(tNode), pointer :: material, &
homogenization, &
constituents
integer :: m
if(maxval(discretization_materialAt) > materials%length) &
call IO_error(155,ext_msg='More materials requested than found in material.yaml')
do m = 1, materials%length
material => materials%get(m)
constituents => material%get('constituents')
homogenization => homogenizations%get(material%get_asString('homogenization'))
if(constituents%length /= homogenization%get_asInt('N_constituents')) call IO_error(148)
enddo
end subroutine sanityCheck
!--------------------------------------------------------------------------------------------------
!> @brief Get all keys from a dictionary (currently with #_ prefix)
!--------------------------------------------------------------------------------------------------
function getKeys(dict)
class(tNode), intent(in) :: dict
character(len=pStringLen), dimension(:), allocatable :: getKeys
integer :: i
character(len=pStringLen) :: sectionName
allocate(getKeys(dict%length))
do i=1, dict%length
write(sectionName,'(i0,a)') i,'_'
getKeys(i) = trim(adjustl(sectionName))//dict%getKey(i) !ToDo: remove prefix
enddo
end function getKeys
end module material end module material

2
src/mesh/FEM_utilities.f90
View File

@ -164,7 +164,7 @@ subroutine utilities_constitutiveResponse(timeinc,P_av,forwardData)
cutBack = .false. ! reset cutBack status cutBack = .false. ! reset cutBack status
P_av = sum(sum(materialpoint_P,dim=4),dim=3) * wgt ! average of P P_av = sum(sum(homogenization_P,dim=4),dim=3) * wgt ! average of P
call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr) call MPI_Allreduce(MPI_IN_PLACE,P_av,9,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
end subroutine utilities_constitutiveResponse end subroutine utilities_constitutiveResponse

20
src/mesh/mesh_mech_FEM.f90
View File

@ -400,15 +400,15 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
(((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp+1)*dimPlex)) (((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp+1)*dimPlex))
enddo enddo
enddo enddo
materialpoint_F(1:dimPlex,1:dimPlex,qPt+1,cell+1) = & homogenization_F(1:dimPlex,1:dimPlex,qPt+1,cell+1) = &
reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex], order=[2,1]) reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex], order=[2,1])
enddo enddo
if (num%BBarStabilisation) then if (num%BBarStabilisation) then
detFAvg = math_det33(sum(materialpoint_F(1:3,1:3,1:nQuadrature,cell+1),dim=3)/real(nQuadrature)) detFAvg = math_det33(sum(homogenization_F(1:3,1:3,1:nQuadrature,cell+1),dim=3)/real(nQuadrature))
do qPt = 1, nQuadrature do qPt = 1, nQuadrature
materialpoint_F(1:dimPlex,1:dimPlex,qPt,cell+1) = & homogenization_F(1:dimPlex,1:dimPlex,qPt,cell+1) = &
materialpoint_F(1:dimPlex,1:dimPlex,qPt,cell+1)* & homogenization_F(1:dimPlex,1:dimPlex,qPt,cell+1)* &
(detFAvg/math_det33(materialpoint_F(1:3,1:3,qPt,cell+1)))**(1.0/real(dimPlex)) (detFAvg/math_det33(homogenization_F(1:3,1:3,qPt,cell+1)))**(1.0/real(dimPlex))
enddo enddo
endif endif
@ -443,7 +443,7 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
enddo enddo
f_scal = f_scal + & f_scal = f_scal + &
matmul(transpose(BMat), & matmul(transpose(BMat), &
reshape(transpose(materialpoint_P(1:dimPlex,1:dimPlex,qPt+1,cell+1)), & reshape(transpose(homogenization_P(1:dimPlex,1:dimPlex,qPt+1,cell+1)), &
shape=[dimPlex*dimPlex]))*qWeights(qPt+1) shape=[dimPlex*dimPlex]))*qWeights(qPt+1)
enddo enddo
f_scal = f_scal*abs(detJ) f_scal = f_scal*abs(detJ)
@ -545,7 +545,7 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
(((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp+1)*dimPlex)) (((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp+1)*dimPlex))
enddo enddo
enddo enddo
MatA = matmul(reshape(reshape(materialpoint_dPdF(1:dimPlex,1:dimPlex,1:dimPlex,1:dimPlex,qPt+1,cell+1), & MatA = matmul(reshape(reshape(homogenization_dPdF(1:dimPlex,1:dimPlex,1:dimPlex,1:dimPlex,qPt+1,cell+1), &
shape=[dimPlex,dimPlex,dimPlex,dimPlex], order=[2,1,4,3]), & shape=[dimPlex,dimPlex,dimPlex,dimPlex], order=[2,1,4,3]), &
shape=[dimPlex*dimPlex,dimPlex*dimPlex]),BMat)*qWeights(qPt+1) shape=[dimPlex*dimPlex,dimPlex*dimPlex]),BMat)*qWeights(qPt+1)
if (num%BBarStabilisation) then if (num%BBarStabilisation) then
@ -553,12 +553,12 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
FInv = math_inv33(F) FInv = math_inv33(F)
K_eA = K_eA + matmul(transpose(BMat),MatA)*math_det33(FInv)**(1.0/real(dimPlex)) K_eA = K_eA + matmul(transpose(BMat),MatA)*math_det33(FInv)**(1.0/real(dimPlex))
K_eB = K_eB - & K_eB = K_eB - &
matmul(transpose(matmul(reshape(materialpoint_F(1:dimPlex,1:dimPlex,qPt+1,cell+1), & matmul(transpose(matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,qPt+1,cell+1), &
shape=[dimPlex*dimPlex,1]), & shape=[dimPlex*dimPlex,1]), &
matmul(reshape(FInv(1:dimPlex,1:dimPlex), & matmul(reshape(FInv(1:dimPlex,1:dimPlex), &
shape=[1,dimPlex*dimPlex],order=[2,1]),BMat))),MatA) shape=[1,dimPlex*dimPlex],order=[2,1]),BMat))),MatA)
MatB = MatB + & MatB = MatB + &
matmul(reshape(materialpoint_F(1:dimPlex,1:dimPlex,qPt+1,cell+1),shape=[1,dimPlex*dimPlex]),MatA) matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,qPt+1,cell+1),shape=[1,dimPlex*dimPlex]),MatA)
FAvg = FAvg + F FAvg = FAvg + F
BMatAvg = BMatAvg + BMat BMatAvg = BMatAvg + BMat
else else
@ -630,7 +630,7 @@ subroutine FEM_mech_forward(guess,timeinc,timeinc_old,fieldBC)
! forward last inc ! forward last inc
if (guess .and. .not. cutBack) then if (guess .and. .not. cutBack) then
ForwardData = .True. ForwardData = .True.
materialpoint_F0 = materialpoint_F homogenization_F0 = homogenization_F
call SNESGetDM(mech_snes,dm_local,ierr); CHKERRQ(ierr) !< retrieve mesh info from mech_snes into dm_local call SNESGetDM(mech_snes,dm_local,ierr); CHKERRQ(ierr) !< retrieve mesh info from mech_snes into dm_local
call DMGetSection(dm_local,section,ierr); CHKERRQ(ierr) call DMGetSection(dm_local,section,ierr); CHKERRQ(ierr)
call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr) call DMGetLocalVector(dm_local,x_local,ierr); CHKERRQ(ierr)

24
src/results.f90
View File

@ -56,7 +56,7 @@ module results
results_addAttribute, & results_addAttribute, &
results_removeLink, & results_removeLink, &
results_mapping_constituent, & results_mapping_constituent, &
results_mapping_materialpoint results_mapping_homogenization
contains contains
subroutine results_init(restart) subroutine results_init(restart)
@ -644,7 +644,7 @@ end subroutine results_mapping_constituent
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief adds the unique mapping from spatial position and constituent ID to results !> @brief adds the unique mapping from spatial position and constituent ID to results
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine results_mapping_materialpoint(homogenizationAt,memberAtLocal,label) subroutine results_mapping_homogenization(homogenizationAt,memberAtLocal,label)
integer, dimension(:), intent(in) :: homogenizationAt !< homogenization section at (element) integer, dimension(:), intent(in) :: homogenizationAt !< homogenization section at (element)
integer, dimension(:,:), intent(in) :: memberAtLocal !< homogenization member at (IP,element) integer, dimension(:,:), intent(in) :: memberAtLocal !< homogenization member at (IP,element)
@ -711,13 +711,13 @@ subroutine results_mapping_materialpoint(homogenizationAt,memberAtLocal,label)
! MPI settings and communication ! MPI settings and communication
#ifdef PETSc #ifdef PETSc
call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, ierr) call h5pset_dxpl_mpio_f(plist_id, H5FD_MPIO_COLLECTIVE_F, ierr)
if (ierr < 0) call IO_error(1,ext_msg='results_mapping_materialpoint: h5pset_dxpl_mpio_f') if (ierr < 0) call IO_error(1,ext_msg='results_mapping_homogenization: h5pset_dxpl_mpio_f')
call MPI_allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get output at each process call MPI_allreduce(MPI_IN_PLACE,writeSize,worldsize,MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr) ! get output at each process
if (ierr /= 0) call IO_error(894,ext_msg='results_mapping_materialpoint: MPI_allreduce/writeSize') if (ierr /= 0) call IO_error(894,ext_msg='results_mapping_homogenization: MPI_allreduce/writeSize')
call MPI_allreduce(MPI_IN_PLACE,memberOffset,size(memberOffset),MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)! get offset at each process call MPI_allreduce(MPI_IN_PLACE,memberOffset,size(memberOffset),MPI_INT,MPI_SUM,PETSC_COMM_WORLD,ierr)! get offset at each process
if (ierr /= 0) call IO_error(894,ext_msg='results_mapping_materialpoint: MPI_allreduce/memberOffset') if (ierr /= 0) call IO_error(894,ext_msg='results_mapping_homogenization: MPI_allreduce/memberOffset')
#endif #endif
myShape = int([writeSize(worldrank)], HSIZE_T) myShape = int([writeSize(worldrank)], HSIZE_T)
@ -727,13 +727,13 @@ subroutine results_mapping_materialpoint(homogenizationAt,memberAtLocal,label)
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! create dataspace in memory (local shape = hyperslab) and in file (global shape) ! create dataspace in memory (local shape = hyperslab) and in file (global shape)
call h5screate_simple_f(1,myShape,memspace_id,ierr,myShape) call h5screate_simple_f(1,myShape,memspace_id,ierr,myShape)
if (ierr < 0) call IO_error(1,ext_msg='results_mapping_materialpoint: h5screate_simple_f/memspace_id') if (ierr < 0) call IO_error(1,ext_msg='results_mapping_homogenization: h5screate_simple_f/memspace_id')
call h5screate_simple_f(1,totalShape,filespace_id,ierr,totalShape) call h5screate_simple_f(1,totalShape,filespace_id,ierr,totalShape)
if (ierr < 0) call IO_error(1,ext_msg='results_mapping_materialpoint: h5screate_simple_f/filespace_id') if (ierr < 0) call IO_error(1,ext_msg='results_mapping_homogenization: h5screate_simple_f/filespace_id')
call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myOffset, myShape, ierr) call h5sselect_hyperslab_f(filespace_id, H5S_SELECT_SET_F, myOffset, myShape, ierr)
if (ierr < 0) call IO_error(1,ext_msg='results_mapping_materialpoint: h5sselect_hyperslab_f') if (ierr < 0) call IO_error(1,ext_msg='results_mapping_homogenization: h5sselect_hyperslab_f')
!--------------------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------------------
! expand phaseAt to consider IPs (is not stored per IP) ! expand phaseAt to consider IPs (is not stored per IP)
@ -753,14 +753,14 @@ subroutine results_mapping_materialpoint(homogenizationAt,memberAtLocal,label)
loc_id = results_openGroup('/mapping') loc_id = results_openGroup('/mapping')
call h5dcreate_f(loc_id, 'homogenization', dtype_id, filespace_id, dset_id, ierr) call h5dcreate_f(loc_id, 'homogenization', dtype_id, filespace_id, dset_id, ierr)
if (ierr < 0) call IO_error(1,ext_msg='results_mapping_materialpoint: h5dcreate_f') if (ierr < 0) call IO_error(1,ext_msg='results_mapping_homogenization: h5dcreate_f')
call h5dwrite_f(dset_id, name_id, reshape(label(pack(homogenizationAtMaterialpoint,.true.)),myShape), & call h5dwrite_f(dset_id, name_id, reshape(label(pack(homogenizationAtMaterialpoint,.true.)),myShape), &
myShape, ierr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) myShape, ierr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
if (ierr < 0) call IO_error(1,ext_msg='results_mapping_materialpoint: h5dwrite_f/name_id') if (ierr < 0) call IO_error(1,ext_msg='results_mapping_homogenization: h5dwrite_f/name_id')
call h5dwrite_f(dset_id, position_id, reshape(pack(memberAtGlobal,.true.),myShape), & call h5dwrite_f(dset_id, position_id, reshape(pack(memberAtGlobal,.true.),myShape), &
myShape, ierr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id) myShape, ierr, file_space_id = filespace_id, mem_space_id = memspace_id, xfer_prp = plist_id)
if (ierr < 0) call IO_error(1,ext_msg='results_mapping_materialpoint: h5dwrite_f/position_id') if (ierr < 0) call IO_error(1,ext_msg='results_mapping_homogenization: h5dwrite_f/position_id')
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
! close all ! close all
@ -776,7 +776,7 @@ subroutine results_mapping_materialpoint(homogenizationAt,memberAtLocal,label)
! for backward compatibility ! for backward compatibility
call results_setLink('/mapping/homogenization','/mapping/cellResults/materialpoint') call results_setLink('/mapping/homogenization','/mapping/cellResults/materialpoint')
end subroutine results_mapping_materialpoint end subroutine results_mapping_homogenization
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------

16
src/source_damage_anisoBrittle.f90
View File

@ -25,7 +25,7 @@ submodule (constitutive:constitutive_damage) source_damage_anisoBrittle
output output
end type tParameters end type tParameters
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances)
contains contains
@ -45,19 +45,19 @@ module function source_damage_anisoBrittle_init(source_length) result(mySources)
phase, & phase, &
sources, & sources, &
src src
integer :: Ninstance,sourceOffset,NipcMyPhase,p integer :: Ninstances,sourceOffset,Nconstituents,p
integer, dimension(:), allocatable :: N_cl integer, dimension(:), allocatable :: N_cl
character(len=pStringLen) :: extmsg = '' character(len=pStringLen) :: extmsg = ''
print'(/,a)', ' <<<+- source_damage_anisoBrittle init -+>>>' print'(/,a)', ' <<<+- source_damage_anisoBrittle init -+>>>'
mySources = source_active('damage_anisoBrittle',source_length) mySources = source_active('damage_anisoBrittle',source_length)
Ninstance = count(mySources) Ninstances = count(mySources)
print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstances == 0) return
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(param(Ninstance)) allocate(param(Ninstances))
allocate(source_damage_anisoBrittle_offset (phases%length), source=0) allocate(source_damage_anisoBrittle_offset (phases%length), source=0)
allocate(source_damage_anisoBrittle_instance(phases%length), source=0) allocate(source_damage_anisoBrittle_instance(phases%length), source=0)
@ -100,8 +100,8 @@ module function source_damage_anisoBrittle_init(source_length) result(mySources)
if (any(prm%g_crit < 0.0_pReal)) extmsg = trim(extmsg)//' g_crit' if (any(prm%g_crit < 0.0_pReal)) extmsg = trim(extmsg)//' g_crit'
if (any(prm%s_crit < 0.0_pReal)) extmsg = trim(extmsg)//' s_crit' if (any(prm%s_crit < 0.0_pReal)) extmsg = trim(extmsg)//' s_crit'
NipcMyPhase = count(material_phaseAt==p) * discretization_nIP Nconstituents = count(material_phaseAt==p) * discretization_nIPs
call constitutive_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0) call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,1,1,0)
sourceState(p)%p(sourceOffset)%atol = src%get_asFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal) sourceState(p)%p(sourceOffset)%atol = src%get_asFloat('anisobrittle_atol',defaultVal=1.0e-3_pReal)
if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol' if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisobrittle_atol'

16
src/source_damage_anisoDuctile.f90
View File

@ -19,7 +19,7 @@ submodule(constitutive:constitutive_damage) source_damage_anisoDuctile
output output
end type tParameters end type tParameters
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances)
contains contains
@ -39,19 +39,19 @@ module function source_damage_anisoDuctile_init(source_length) result(mySources)
pl, & pl, &
sources, & sources, &
src src
integer :: Ninstance,sourceOffset,NipcMyPhase,p integer :: Ninstances,sourceOffset,Nconstituents,p
integer, dimension(:), allocatable :: N_sl integer, dimension(:), allocatable :: N_sl
character(len=pStringLen) :: extmsg = '' character(len=pStringLen) :: extmsg = ''
print'(/,a)', ' <<<+- source_damage_anisoDuctile init -+>>>' print'(/,a)', ' <<<+- source_damage_anisoDuctile init -+>>>'
mySources = source_active('damage_anisoDuctile',source_length) mySources = source_active('damage_anisoDuctile',source_length)
Ninstance = count(mySources) Ninstances = count(mySources)
print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstances == 0) return
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(param(Ninstance)) allocate(param(Ninstances))
allocate(source_damage_anisoDuctile_offset (phases%length), source=0) allocate(source_damage_anisoDuctile_offset (phases%length), source=0)
allocate(source_damage_anisoDuctile_instance(phases%length), source=0) allocate(source_damage_anisoDuctile_instance(phases%length), source=0)
@ -84,8 +84,8 @@ module function source_damage_anisoDuctile_init(source_length) result(mySources)
if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q' if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q'
if (any(prm%gamma_crit < 0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit' if (any(prm%gamma_crit < 0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit'
NipcMyPhase=count(material_phaseAt==p) * discretization_nIP Nconstituents=count(material_phaseAt==p) * discretization_nIPs
call constitutive_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0) call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,1,1,0)
sourceState(p)%p(sourceOffset)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal) sourceState(p)%p(sourceOffset)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal)
if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol' if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol'

16
src/source_damage_isoBrittle.f90
View File

@ -17,7 +17,7 @@ submodule(constitutive:constitutive_damage) source_damage_isoBrittle
output output
end type tParameters end type tParameters
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances)
contains contains
@ -36,18 +36,18 @@ module function source_damage_isoBrittle_init(source_length) result(mySources)
phase, & phase, &
sources, & sources, &
src src
integer :: Ninstance,sourceOffset,NipcMyPhase,p integer :: Ninstances,sourceOffset,Nconstituents,p
character(len=pStringLen) :: extmsg = '' character(len=pStringLen) :: extmsg = ''
print'(/,a)', ' <<<+- source_damage_isoBrittle init -+>>>' print'(/,a)', ' <<<+- source_damage_isoBrittle init -+>>>'
mySources = source_active('damage_isoBrittle',source_length) mySources = source_active('damage_isoBrittle',source_length)
Ninstance = count(mySources) Ninstances = count(mySources)
print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstances == 0) return
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(param(Ninstance)) allocate(param(Ninstances))
allocate(source_damage_isoBrittle_offset (phases%length), source=0) allocate(source_damage_isoBrittle_offset (phases%length), source=0)
allocate(source_damage_isoBrittle_instance(phases%length), source=0) allocate(source_damage_isoBrittle_instance(phases%length), source=0)
@ -73,8 +73,8 @@ module function source_damage_isoBrittle_init(source_length) result(mySources)
! sanity checks ! sanity checks
if (prm%W_crit <= 0.0_pReal) extmsg = trim(extmsg)//' W_crit' if (prm%W_crit <= 0.0_pReal) extmsg = trim(extmsg)//' W_crit'
NipcMyPhase = count(material_phaseAt==p) * discretization_nIP Nconstituents = count(material_phaseAt==p) * discretization_nIPs
call constitutive_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,1) call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,1,1,1)
sourceState(p)%p(sourceOffset)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal) sourceState(p)%p(sourceOffset)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal)
if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol' if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol'

16
src/source_damage_isoDuctile.f90
View File

@ -18,7 +18,7 @@ submodule (constitutive:constitutive_damage) source_damage_isoDuctile
output output
end type tParameters end type tParameters
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances)
contains contains
@ -38,18 +38,18 @@ module function source_damage_isoDuctile_init(source_length) result(mySources)
phase, & phase, &
sources, & sources, &
src src
integer :: Ninstance,sourceOffset,NipcMyPhase,p integer :: Ninstances,sourceOffset,Nconstituents,p
character(len=pStringLen) :: extmsg = '' character(len=pStringLen) :: extmsg = ''
print'(/,a)', ' <<<+- source_damage_isoDuctile init -+>>>' print'(/,a)', ' <<<+- source_damage_isoDuctile init -+>>>'
mySources = source_active('damage_isoDuctile',source_length) mySources = source_active('damage_isoDuctile',source_length)
Ninstance = count(mySources) Ninstances = count(mySources)
print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstances == 0) return
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(param(Ninstance)) allocate(param(Ninstances))
allocate(source_damage_isoDuctile_offset (phases%length), source=0) allocate(source_damage_isoDuctile_offset (phases%length), source=0)
allocate(source_damage_isoDuctile_instance(phases%length), source=0) allocate(source_damage_isoDuctile_instance(phases%length), source=0)
@ -77,8 +77,8 @@ module function source_damage_isoDuctile_init(source_length) result(mySources)
if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q' if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q'
if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit' if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit'
NipcMyPhase=count(material_phaseAt==p) * discretization_nIP Nconstituents=count(material_phaseAt==p) * discretization_nIPs
call constitutive_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0) call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,1,1,0)
sourceState(p)%p(sourceOffset)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal) sourceState(p)%p(sourceOffset)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal)
if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol' if(any(sourceState(p)%p(sourceOffset)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol'

18
src/source_thermal_dissipation.f90
View File

@ -15,7 +15,7 @@ submodule(constitutive:constitutive_thermal) source_thermal_dissipation
kappa !< TAYLOR-QUINNEY factor kappa !< TAYLOR-QUINNEY factor
end type tParameters end type tParameters
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances)
contains contains
@ -35,17 +35,17 @@ module function source_thermal_dissipation_init(source_length) result(mySources)
phase, & phase, &
sources, & sources, &
src src
integer :: Ninstance,sourceOffset,NipcMyPhase,p integer :: Ninstances,sourceOffset,Nconstituents,p
print'(/,a)', ' <<<+- source_thermal_dissipation init -+>>>' print'(/,a)', ' <<<+- source_thermal_dissipation init -+>>>'
mySources = source_active('thermal_dissipation',source_length) mySources = source_active('thermal_dissipation',source_length)
Ninstance = count(mySources) Ninstances = count(mySources)
print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstances == 0) return
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(param(Ninstance)) allocate(param(Ninstances))
allocate(source_thermal_dissipation_offset (phases%length), source=0) allocate(source_thermal_dissipation_offset (phases%length), source=0)
allocate(source_thermal_dissipation_instance(phases%length), source=0) allocate(source_thermal_dissipation_instance(phases%length), source=0)
@ -61,8 +61,8 @@ module function source_thermal_dissipation_init(source_length) result(mySources)
src => sources%get(sourceOffset) src => sources%get(sourceOffset)
prm%kappa = src%get_asFloat('kappa') prm%kappa = src%get_asFloat('kappa')
NipcMyPhase = count(material_phaseAt==p) * discretization_nIP Nconstituents = count(material_phaseAt==p) * discretization_nIPs
call constitutive_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,0,0,0) call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,0,0,0)
end associate end associate
endif endif
@ -74,7 +74,7 @@ end function source_thermal_dissipation_init
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
!> @brief Ninstances dissipation rate !> @brief Ninstancess dissipation rate
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
module subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT, Tstar, Lp, phase) module subroutine source_thermal_dissipation_getRateAndItsTangent(TDot, dTDot_dT, Tstar, Lp, phase)

16
src/source_thermal_externalheat.f90
View File

@ -19,7 +19,7 @@ submodule(constitutive:constitutive_thermal) source_thermal_externalheat
nIntervals nIntervals
end type tParameters end type tParameters
type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstance) type(tParameters), dimension(:), allocatable :: param !< containers of constitutive parameters (len Ninstances)
contains contains
@ -39,17 +39,17 @@ module function source_thermal_externalheat_init(source_length) result(mySources
phase, & phase, &
sources, & sources, &
src src
integer :: Ninstance,sourceOffset,NipcMyPhase,p integer :: Ninstances,sourceOffset,Nconstituents,p
print'(/,a)', ' <<<+- source_thermal_externalHeat init -+>>>' print'(/,a)', ' <<<+- source_thermal_externalHeat init -+>>>'
mySources = source_active('thermal_externalheat',source_length) mySources = source_active('thermal_externalheat',source_length)
Ninstance = count(mySources) Ninstances = count(mySources)
print'(a,i2)', ' # instances: ',Ninstance; flush(IO_STDOUT) print'(a,i2)', ' # instances: ',Ninstances; flush(IO_STDOUT)
if(Ninstance == 0) return if(Ninstances == 0) return
phases => config_material%get('phase') phases => config_material%get('phase')
allocate(param(Ninstance)) allocate(param(Ninstances))
allocate(source_thermal_externalheat_offset (phases%length), source=0) allocate(source_thermal_externalheat_offset (phases%length), source=0)
allocate(source_thermal_externalheat_instance(phases%length), source=0) allocate(source_thermal_externalheat_instance(phases%length), source=0)
@ -69,8 +69,8 @@ module function source_thermal_externalheat_init(source_length) result(mySources
prm%f_T = src%get_asFloats('f_T',requiredSize = size(prm%t_n)) prm%f_T = src%get_asFloats('f_T',requiredSize = size(prm%t_n))
NipcMyPhase = count(material_phaseAt==p) * discretization_nIP Nconstituents = count(material_phaseAt==p) * discretization_nIPs
call constitutive_allocateState(sourceState(p)%p(sourceOffset),NipcMyPhase,1,1,0) call constitutive_allocateState(sourceState(p)%p(sourceOffset),Nconstituents,1,1,0)
end associate end associate
endif endif

32
src/thermal_adiabatic.f90
View File

@ -40,7 +40,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine thermal_adiabatic_init subroutine thermal_adiabatic_init
integer :: maxNinstance,h,NofMyHomog integer :: maxNinstances,h,Nmaterialpoints
class(tNode), pointer :: & class(tNode), pointer :: &
material_homogenization, & material_homogenization, &
homog, & homog, &
@ -48,13 +48,13 @@ subroutine thermal_adiabatic_init
print'(/,a)', ' <<<+- thermal_adiabatic init -+>>>'; flush(6) print'(/,a)', ' <<<+- thermal_adiabatic init -+>>>'; flush(6)
maxNinstance = count(thermal_type == THERMAL_adiabatic_ID) maxNinstances = count(thermal_type == THERMAL_adiabatic_ID)
if (maxNinstance == 0) return if (maxNinstances == 0) return
allocate(param(maxNinstance)) allocate(param(maxNinstances))
material_homogenization => config_material%get('homogenization') material_homogenization => config_material%get('homogenization')
do h = 1, material_Nhomogenization do h = 1, size(material_name_homogenization)
if (thermal_type(h) /= THERMAL_adiabatic_ID) cycle if (thermal_type(h) /= THERMAL_adiabatic_ID) cycle
homog => material_homogenization%get(h) homog => material_homogenization%get(h)
homogThermal => homog%get('thermal') homogThermal => homog%get('thermal')
@ -67,17 +67,17 @@ subroutine thermal_adiabatic_init
prm%output = homogThermal%get_asStrings('output',defaultVal=emptyStringArray) prm%output = homogThermal%get_asStrings('output',defaultVal=emptyStringArray)
#endif #endif
NofMyHomog=count(material_homogenizationAt==h) Nmaterialpoints=count(material_homogenizationAt==h)
thermalState(h)%sizeState = 1 thermalState(h)%sizeState = 1
allocate(thermalState(h)%state0 (1,NofMyHomog), source=thermal_initialT(h)) allocate(thermalState(h)%state0 (1,Nmaterialpoints), source=thermal_initialT(h))
allocate(thermalState(h)%subState0(1,NofMyHomog), source=thermal_initialT(h)) allocate(thermalState(h)%subState0(1,Nmaterialpoints), source=thermal_initialT(h))
allocate(thermalState(h)%state (1,NofMyHomog), source=thermal_initialT(h)) allocate(thermalState(h)%state (1,Nmaterialpoints), source=thermal_initialT(h))
thermalMapping(h)%p => material_homogenizationMemberAt thermalMapping(h)%p => material_homogenizationMemberAt
deallocate(temperature(h)%p) deallocate(temperature(h)%p)
temperature(h)%p => thermalState(h)%state(1,:) temperature(h)%p => thermalState(h)%state(1,:)
deallocate(temperatureRate(h)%p) deallocate(temperatureRate(h)%p)
allocate (temperatureRate(h)%p(NofMyHomog), source=0.0_pReal) allocate (temperatureRate(h)%p(Nmaterialpoints), source=0.0_pReal)
end associate end associate
enddo enddo
@ -145,8 +145,8 @@ subroutine thermal_adiabatic_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
homog = material_homogenizationAt(el) 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, crystallite_S, crystallite_Lp, ip, el)
Tdot = Tdot/real(homogenization_Ngrains(homog),pReal) Tdot = Tdot/real(homogenization_Nconstituents(homog),pReal)
dTdot_dT = dTdot_dT/real(homogenization_Ngrains(homog),pReal) dTdot_dT = dTdot_dT/real(homogenization_Nconstituents(homog),pReal)
end subroutine thermal_adiabatic_getSourceAndItsTangent end subroutine thermal_adiabatic_getSourceAndItsTangent
@ -167,13 +167,13 @@ function thermal_adiabatic_getSpecificHeat(ip,el)
thermal_adiabatic_getSpecificHeat = 0.0_pReal thermal_adiabatic_getSpecificHeat = 0.0_pReal
do grain = 1, homogenization_Ngrains(material_homogenizationAt(el)) do grain = 1, homogenization_Nconstituents(material_homogenizationAt(el))
thermal_adiabatic_getSpecificHeat = thermal_adiabatic_getSpecificHeat & thermal_adiabatic_getSpecificHeat = thermal_adiabatic_getSpecificHeat &
+ lattice_c_p(material_phaseAt(grain,el)) + lattice_c_p(material_phaseAt(grain,el))
enddo enddo
thermal_adiabatic_getSpecificHeat = thermal_adiabatic_getSpecificHeat & thermal_adiabatic_getSpecificHeat = thermal_adiabatic_getSpecificHeat &
/ real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) / real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)
end function thermal_adiabatic_getSpecificHeat end function thermal_adiabatic_getSpecificHeat
@ -193,13 +193,13 @@ function thermal_adiabatic_getMassDensity(ip,el)
thermal_adiabatic_getMassDensity = 0.0_pReal thermal_adiabatic_getMassDensity = 0.0_pReal
do grain = 1, homogenization_Ngrains(material_homogenizationAt(el)) do grain = 1, homogenization_Nconstituents(material_homogenizationAt(el))
thermal_adiabatic_getMassDensity = thermal_adiabatic_getMassDensity & thermal_adiabatic_getMassDensity = thermal_adiabatic_getMassDensity &
+ lattice_rho(material_phaseAt(grain,el)) + lattice_rho(material_phaseAt(grain,el))
enddo enddo
thermal_adiabatic_getMassDensity = thermal_adiabatic_getMassDensity & thermal_adiabatic_getMassDensity = thermal_adiabatic_getMassDensity &
/ real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) / real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)
end function thermal_adiabatic_getMassDensity end function thermal_adiabatic_getMassDensity

36
src/thermal_conduction.f90
View File

@ -41,7 +41,7 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine thermal_conduction_init subroutine thermal_conduction_init
integer :: Ninstance,NofMyHomog,h integer :: Ninstances,Nmaterialpoints,h
class(tNode), pointer :: & class(tNode), pointer :: &
material_homogenization, & material_homogenization, &
homog, & homog, &
@ -49,11 +49,11 @@ subroutine thermal_conduction_init
print'(/,a)', ' <<<+- thermal_conduction init -+>>>'; flush(6) print'(/,a)', ' <<<+- thermal_conduction init -+>>>'; flush(6)
Ninstance = count(thermal_type == THERMAL_conduction_ID) Ninstances = count(thermal_type == THERMAL_conduction_ID)
allocate(param(Ninstance)) allocate(param(Ninstances))
material_homogenization => config_material%get('homogenization') material_homogenization => config_material%get('homogenization')
do h = 1, material_Nhomogenization do h = 1, size(material_name_homogenization)
if (thermal_type(h) /= THERMAL_conduction_ID) cycle if (thermal_type(h) /= THERMAL_conduction_ID) cycle
homog => material_homogenization%get(h) homog => material_homogenization%get(h)
homogThermal => homog%get('thermal') homogThermal => homog%get('thermal')
@ -65,17 +65,17 @@ subroutine thermal_conduction_init
prm%output = homogThermal%get_asStrings('output',defaultVal=emptyStringArray) prm%output = homogThermal%get_asStrings('output',defaultVal=emptyStringArray)
#endif #endif
NofMyHomog=count(material_homogenizationAt==h) Nmaterialpoints=count(material_homogenizationAt==h)
thermalState(h)%sizeState = 0 thermalState(h)%sizeState = 0
allocate(thermalState(h)%state0 (0,NofMyHomog)) allocate(thermalState(h)%state0 (0,Nmaterialpoints))
allocate(thermalState(h)%subState0(0,NofMyHomog)) allocate(thermalState(h)%subState0(0,Nmaterialpoints))
allocate(thermalState(h)%state (0,NofMyHomog)) allocate(thermalState(h)%state (0,Nmaterialpoints))
thermalMapping(h)%p => material_homogenizationMemberAt thermalMapping(h)%p => material_homogenizationMemberAt
deallocate(temperature (h)%p) deallocate(temperature (h)%p)
allocate (temperature (h)%p(NofMyHomog), source=thermal_initialT(h)) allocate (temperature (h)%p(Nmaterialpoints), source=thermal_initialT(h))
deallocate(temperatureRate(h)%p) deallocate(temperatureRate(h)%p)
allocate (temperatureRate(h)%p(NofMyHomog), source=0.0_pReal) allocate (temperatureRate(h)%p(Nmaterialpoints), source=0.0_pReal)
end associate end associate
enddo enddo
@ -104,8 +104,8 @@ subroutine thermal_conduction_getSourceAndItsTangent(Tdot, dTdot_dT, T, ip, el)
homog = material_homogenizationAt(el) 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, crystallite_S,crystallite_Lp ,ip, el)
Tdot = Tdot/real(homogenization_Ngrains(homog),pReal) Tdot = Tdot/real(homogenization_Nconstituents(homog),pReal)
dTdot_dT = dTdot_dT/real(homogenization_Ngrains(homog),pReal) dTdot_dT = dTdot_dT/real(homogenization_Nconstituents(homog),pReal)
end subroutine thermal_conduction_getSourceAndItsTangent end subroutine thermal_conduction_getSourceAndItsTangent
@ -125,13 +125,13 @@ function thermal_conduction_getConductivity(ip,el)
thermal_conduction_getConductivity = 0.0_pReal thermal_conduction_getConductivity = 0.0_pReal
do grain = 1, homogenization_Ngrains(material_homogenizationAt(el)) do grain = 1, homogenization_Nconstituents(material_homogenizationAt(el))
thermal_conduction_getConductivity = thermal_conduction_getConductivity + & thermal_conduction_getConductivity = thermal_conduction_getConductivity + &
crystallite_push33ToRef(grain,ip,el,lattice_K(:,:,material_phaseAt(grain,el))) crystallite_push33ToRef(grain,ip,el,lattice_K(:,:,material_phaseAt(grain,el)))
enddo enddo
thermal_conduction_getConductivity = thermal_conduction_getConductivity & thermal_conduction_getConductivity = thermal_conduction_getConductivity &
/ real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) / real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)
end function thermal_conduction_getConductivity end function thermal_conduction_getConductivity
@ -151,13 +151,13 @@ function thermal_conduction_getSpecificHeat(ip,el)
thermal_conduction_getSpecificHeat = 0.0_pReal thermal_conduction_getSpecificHeat = 0.0_pReal
do grain = 1, homogenization_Ngrains(material_homogenizationAt(el)) do grain = 1, homogenization_Nconstituents(material_homogenizationAt(el))
thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat & thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat &
+ lattice_c_p(material_phaseAt(grain,el)) + lattice_c_p(material_phaseAt(grain,el))
enddo enddo
thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat & thermal_conduction_getSpecificHeat = thermal_conduction_getSpecificHeat &
/ real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) / real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)
end function thermal_conduction_getSpecificHeat end function thermal_conduction_getSpecificHeat
@ -178,13 +178,13 @@ function thermal_conduction_getMassDensity(ip,el)
thermal_conduction_getMassDensity = 0.0_pReal thermal_conduction_getMassDensity = 0.0_pReal
do grain = 1, homogenization_Ngrains(material_homogenizationAt(el)) do grain = 1, homogenization_Nconstituents(material_homogenizationAt(el))
thermal_conduction_getMassDensity = thermal_conduction_getMassDensity & thermal_conduction_getMassDensity = thermal_conduction_getMassDensity &
+ lattice_rho(material_phaseAt(grain,el)) + lattice_rho(material_phaseAt(grain,el))
enddo enddo
thermal_conduction_getMassDensity = thermal_conduction_getMassDensity & thermal_conduction_getMassDensity = thermal_conduction_getMassDensity &
/ real(homogenization_Ngrains(material_homogenizationAt(el)),pReal) / real(homogenization_Nconstituents(material_homogenizationAt(el)),pReal)
end function thermal_conduction_getMassDensity end function thermal_conduction_getMassDensity

12
src/thermal_isothermal.f90
View File

@ -16,18 +16,18 @@ contains
!-------------------------------------------------------------------------------------------------- !--------------------------------------------------------------------------------------------------
subroutine thermal_isothermal_init subroutine thermal_isothermal_init
integer :: h,NofMyHomog integer :: h,Nmaterialpoints
print'(/,a)', ' <<<+- thermal_isothermal init -+>>>'; flush(6) print'(/,a)', ' <<<+- thermal_isothermal init -+>>>'; flush(6)
do h = 1, material_Nhomogenization do h = 1, size(material_name_homogenization)
if (thermal_type(h) /= THERMAL_isothermal_ID) cycle if (thermal_type(h) /= THERMAL_isothermal_ID) cycle
NofMyHomog = count(material_homogenizationAt == h) Nmaterialpoints = count(material_homogenizationAt == h)
thermalState(h)%sizeState = 0 thermalState(h)%sizeState = 0
allocate(thermalState(h)%state0 (0,NofMyHomog)) allocate(thermalState(h)%state0 (0,Nmaterialpoints))
allocate(thermalState(h)%subState0(0,NofMyHomog)) allocate(thermalState(h)%subState0(0,Nmaterialpoints))
allocate(thermalState(h)%state (0,NofMyHomog)) allocate(thermalState(h)%state (0,Nmaterialpoints))
deallocate(temperature (h)%p) deallocate(temperature (h)%p)
allocate (temperature (h)%p(1), source=thermal_initialT(h)) allocate (temperature (h)%p(1), source=thermal_initialT(h))