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

Merge branch 'plain-field-variables' into 'development' 

multiphysics restructure

See merge request damask/DAMASK!304
This commit is contained in:
Sharan Roongta 2020-12-21 16:24:18 +01:00
parent 35b833e2ad 3884549e19
commit d1c7ec0680
32 changed files with 393 additions and 359 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

22
src/CPFEM.f90
View File

@ -89,8 +89,8 @@ subroutine CPFEM_initAll
call lattice_init
call material_init(.false.)
call constitutive_init
call crystallite_init
call homogenization_init
call crystallite_init
call CPFEM_init
call config_deallocate
@ -153,7 +153,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
H
integer(pInt) elCP, & ! crystal plasticity element number
i, j, k, l, m, n, ph, homog, mySource
i, j, k, l, m, n, ph, homog, mySource,ma
real(pReal), parameter :: ODD_STRESS = 1e15_pReal, & !< return value for stress if terminallyIll
ODD_JACOBIAN = 1e50_pReal !< return value for jacobian if terminallyIll
@ -161,6 +161,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
elCP = mesh_FEM2DAMASK_elem(elFE)
ma = (elCP-1) * discretization_nIPs + ip
if (debugCPFEM%basic .and. elCP == debugCPFEM%element .and. ip == debugCPFEM%ip) then
print'(/,a)', '#############################################'
print'(a1,a22,1x,i8,a13)', '#','element', elCP, '#'
@ -181,11 +183,11 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
chosenThermal1: select case (thermal_type(material_homogenizationAt(elCP)))
case (THERMAL_conduction_ID) chosenThermal1
temperature(material_homogenizationAt(elCP))%p(thermalMapping(material_homogenizationAt(elCP))%p(ip,elCP)) = &
temperature(material_homogenizationAt(elCP))%p(material_homogenizationMemberAt(ip,elCP)) = &
temperature_inp
end select chosenThermal1
homogenization_F0(1:3,1:3,ip,elCP) = ffn
homogenization_F(1:3,1:3,ip,elCP) = ffn1
homogenization_F0(1:3,1:3,ma) = ffn
homogenization_F(1:3,1:3,ma) = ffn1
if (iand(mode, CPFEM_CALCRESULTS) /= 0_pInt) then
@ -212,17 +214,17 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
else terminalIllness
! translate from P to sigma
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(homogenization_F(1:3,1:3,ip,elCP))
Kirchhoff = matmul(homogenization_P(1:3,1:3,ma), transpose(homogenization_F(1:3,1:3,ma)))
J_inverse = 1.0_pReal / math_det33(homogenization_F(1:3,1:3,ma))
CPFEM_cs(1:6,ip,elCP) = math_sym33to6(J_inverse * Kirchhoff,weighted=.false.)
! translate from dP/dF to dCS/dE
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
H(i,j,k,l) = H(i,j,k,l) &
+ homogenization_F(j,m,ip,elCP) * homogenization_F(l,n,ip,elCP) &
* homogenization_dPdF(i,m,k,n,ip,elCP) &
- math_delta(j,l) * homogenization_F(i,m,ip,elCP) * homogenization_P(k,m,ip,elCP) &
+ homogenization_F(j,m,ma) * homogenization_F(l,n,ma) &
* homogenization_dPdF(i,m,k,n,ma) &
- math_delta(j,l) * homogenization_F(i,m,ma) * homogenization_P(k,m,ma) &
+ 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))
enddo; enddo; enddo; enddo; enddo; enddo

2
src/CPFEM2.f90
View File

@ -63,8 +63,8 @@ subroutine CPFEM_initAll
#endif
call material_init(restart=interface_restartInc>0)
call constitutive_init
call crystallite_init
call homogenization_init
call crystallite_init
call CPFEM_init
call config_deallocate

1
src/commercialFEM_fileList.f90
View File

@ -52,6 +52,7 @@
#include "damage_local.f90"
#include "damage_nonlocal.f90"
#include "homogenization.f90"
#include "homogenization_mech.f90"
#include "homogenization_mech_none.f90"
#include "homogenization_mech_isostrain.f90"
#include "homogenization_mech_RGC.f90"

24
src/constitutive.f90
View File

@ -19,25 +19,21 @@ module constitutive
implicit none
private
integer(kind(ELASTICITY_undefined_ID)), dimension(:), allocatable :: & !ToDo: old intel compiler complains about protected
phase_elasticity !< elasticity of each phase
integer(kind(PLASTICITY_undefined_ID)), dimension(:), allocatable :: & !ToDo: old intel compiler complains about protected
integer(kind(PLASTICITY_undefined_ID)), dimension(:), allocatable :: &
phase_plasticity !< plasticity of each phase
integer(kind(SOURCE_undefined_ID)), dimension(:,:), allocatable :: & ! ToDo: old intel compiler complains about protected
integer(kind(SOURCE_undefined_ID)), dimension(:,:), allocatable :: &
phase_source, & !< active sources mechanisms of each phase
phase_kinematics, & !< active kinematic mechanisms of each phase
phase_stiffnessDegradation !< active stiffness degradation mechanisms of each phase
phase_kinematics !< active kinematic mechanisms of each phase
integer, dimension(:), allocatable, public :: & ! ToDo: old intel compiler complains about protected
integer, dimension(:), allocatable, public :: & !< ToDo: should be protected (bug in Intel compiler)
phase_Nsources, & !< number of source mechanisms active in each phase
phase_Nkinematics, & !< number of kinematic mechanisms active in each phase
phase_NstiffnessDegradations, & !< number of stiffness degradation mechanisms active in each phase
phase_plasticityInstance, & !< instance of particular plasticity of each phase
phase_elasticityInstance !< instance of particular elasticity of each phase
logical, dimension(:), allocatable, public :: & ! ToDo: old intel compiler complains about protected
logical, dimension(:), allocatable, public :: & ! ToDo: should be protected (bug in Intel Compiler)
phase_localPlasticity !< flags phases with local constitutive law
type(tPlasticState), allocatable, dimension(:), public :: &
@ -634,10 +630,10 @@ pure function constitutive_initialFi(ipc, ip, el)
KinematicsLoop: do k = 1, phase_Nkinematics(phase) !< Warning: small initial strain assumption
kinematicsType: select case (phase_kinematics(k,phase))
case (KINEMATICS_thermal_expansion_ID) kinematicsType
homog = material_homogenizationAt(el)
offset = thermalMapping(homog)%p(ip,el)
constitutive_initialFi = &
constitutive_initialFi + kinematics_thermal_expansion_initialStrain(homog,phase,offset)
homog = material_homogenizationAt(el)
offset = material_homogenizationMemberAt(ip,el)
constitutive_initialFi = constitutive_initialFi &
+ kinematics_thermal_expansion_initialStrain(homog,phase,offset)
end select kinematicsType
enddo KinematicsLoop
@ -674,7 +670,7 @@ function constitutive_collectDotState(S, FArray, Fi, FpArray, subdt, ipc, ip, el
logical :: broken
ho = material_homogenizationAt(el)
tme = thermalMapping(ho)%p(ip,el)
tme = material_homogenizationMemberAt(ip,el)
instance = phase_plasticityInstance(phase)
Mp = matmul(matmul(transpose(Fi),Fi),S)

12
src/constitutive_mech.f90
View File

@ -3,6 +3,12 @@
!----------------------------------------------------------------------------------------------------
submodule(constitutive) constitutive_mech
integer(kind(ELASTICITY_undefined_ID)), dimension(:), allocatable :: &
phase_elasticity !< elasticity of each phase
integer(kind(SOURCE_undefined_ID)), dimension(:,:), allocatable :: &
phase_stiffnessDegradation !< active stiffness degradation mechanisms of each phase
interface
module function plastic_none_init() result(myPlasticity)
@ -326,7 +332,7 @@ module subroutine constitutive_hooke_SandItsTangents(S, dS_dFe, dS_dFi, &
DegradationLoop: do d = 1, phase_NstiffnessDegradations(material_phaseAt(ipc,el))
degradationType: select case(phase_stiffnessDegradation(d,material_phaseAt(ipc,el)))
case (STIFFNESS_DEGRADATION_damage_ID) degradationType
C = C * damage(ho)%p(damageMapping(ho)%p(ip,el))**2
C = C * damage(ho)%p(material_homogenizationMemberAt(ip,el))**2
end select degradationType
enddo DegradationLoop
@ -360,7 +366,7 @@ module subroutine constitutive_plastic_dependentState(F, Fp, ipc, ip, el)
instance, of
ho = material_homogenizationAt(el)
tme = thermalMapping(ho)%p(ip,el)
tme = material_homogenizationMemberAt(ip,el)
of = material_phasememberAt(ipc,ip,el)
instance = phase_plasticityInstance(material_phaseAt(ipc,el))
@ -407,7 +413,7 @@ module subroutine constitutive_plastic_LpAndItsTangents(Lp, dLp_dS, dLp_dFi, &
i, j, instance, of
ho = material_homogenizationAt(el)
tme = thermalMapping(ho)%p(ip,el)
tme = material_homogenizationMemberAt(ip,el)
Mp = matmul(matmul(transpose(Fi),Fi),S)
of = material_phasememberAt(ipc,ip,el)

43
src/crystallite.f90
View File

@ -69,9 +69,9 @@ module crystallite
real(pReal), dimension(:,:,:,:,:), allocatable, public :: &
crystallite_partitionedF !< def grad to be reached at end of homog inc
logical, dimension(:,:,:), allocatable, public :: &
logical, dimension(:,:,:), allocatable, public :: &
crystallite_requested !< used by upper level (homogenization) to request crystallite calculation
logical, dimension(:,:,:), allocatable :: &
logical, dimension(:,:,:), allocatable :: &
crystallite_converged !< convergence flag
type :: tOutput !< new requested output (per phase)
@ -135,8 +135,8 @@ contains
!--------------------------------------------------------------------------------------------------
subroutine crystallite_init
logical, dimension(discretization_nIPs,discretization_Nelems) :: devNull
integer :: &
p, &
c, & !< counter in integration point component loop
i, & !< counter in integration point loop
e, & !< counter in element loop
@ -238,19 +238,25 @@ subroutine crystallite_init
phases => config_material%get('phase')
allocate(output_constituent(phases%length))
do c = 1, phases%length
phase => phases%get(c)
do p = 1, phases%length
phase => phases%get(p)
mech => phase%get('mechanics',defaultVal = emptyDict)
#if defined(__GFORTRAN__)
output_constituent(c)%label = output_asStrings(mech)
output_constituent(p)%label = output_asStrings(mech)
#else
output_constituent(c)%label = mech%get_asStrings('output',defaultVal=emptyStringArray)
output_constituent(p)%label = mech%get_asStrings('output',defaultVal=emptyStringArray)
#endif
enddo
#ifdef DEBUG
if (debugCrystallite%basic) then
print'(a42,1x,i10)', ' # of elements: ', eMax
print'(a42,1x,i10)', ' # of integration points/element: ', iMax
print'(a42,1x,i10)', 'max # of constituents/integration point: ', cMax
flush(IO_STDOUT)
endif
#endif
!--------------------------------------------------------------------------------------------------
! initialize
!$OMP PARALLEL DO PRIVATE(i,c)
do e = FEsolving_execElem(1),FEsolving_execElem(2)
do i = FEsolving_execIP(1), FEsolving_execIP(2); do c = 1, homogenization_Nconstituents(material_homogenizationAt(e))
@ -288,16 +294,6 @@ subroutine crystallite_init
enddo
!$OMP END PARALLEL DO
devNull = crystallite_stress()
#ifdef DEBUG
if (debugCrystallite%basic) then
print'(a42,1x,i10)', ' # of elements: ', eMax
print'(a42,1x,i10)', ' # of integration points/element: ', iMax
print'(a42,1x,i10)', 'max # of constituents/integration point: ', cMax
flush(IO_STDOUT)
endif
#endif
end subroutine crystallite_init
@ -321,14 +317,11 @@ function crystallite_stress()
subLp0,& !< plastic velocity grad at start of crystallite inc
subLi0 !< intermediate velocity grad at start of crystallite inc
todo = .false.
subLp0 = crystallite_partitionedLp0
subLi0 = crystallite_partitionedLi0
!--------------------------------------------------------------------------------------------------
! initialize to starting condition
crystallite_subStep = 0.0_pReal
@ -435,8 +428,6 @@ function crystallite_stress()
! integrate --- requires fully defined state array (basic + dependent state)
where(.not. crystallite_converged .and. crystallite_subStep > num%subStepMinCryst) & ! do not try non-converged but fully cutbacked any further
todo = .true. ! TODO: again unroll this into proper elementloop to avoid N^2 for single point evaluation
enddo cutbackLooping
! return whether converged or not
@ -471,10 +462,10 @@ subroutine crystallite_initializeRestorationPoints(i,e)
crystallite_partitionedS0(1:3,1:3,c,i,e) = crystallite_S0(1:3,1:3,c,i,e)
plasticState(material_phaseAt(c,e))%partitionedState0(:,material_phasememberAt(c,i,e)) = &
plasticState(material_phaseAt(c,e))%state0( :,material_phasememberAt(c,i,e))
plasticState(material_phaseAt(c,e))%state0( :,material_phasememberAt(c,i,e))
do s = 1, phase_Nsources(material_phaseAt(c,e))
sourceState(material_phaseAt(c,e))%p(s)%partitionedState0(:,material_phasememberAt(c,i,e)) = &
sourceState(material_phaseAt(c,e))%p(s)%state0( :,material_phasememberAt(c,i,e))
sourceState(material_phaseAt(c,e))%p(s)%state0( :,material_phasememberAt(c,i,e))
enddo
enddo

9
src/damage_local.f90
View File

@ -75,13 +75,10 @@ subroutine damage_local_init
Nmaterialpoints = count(material_homogenizationAt == h)
damageState(h)%sizeState = 1
allocate(damageState(h)%state0 (1,Nmaterialpoints), source=damage_initialPhi(h))
allocate(damageState(h)%subState0(1,Nmaterialpoints), source=damage_initialPhi(h))
allocate(damageState(h)%state (1,Nmaterialpoints), source=damage_initialPhi(h))
allocate(damageState(h)%state0 (1,Nmaterialpoints), source=1.0_pReal)
allocate(damageState(h)%subState0(1,Nmaterialpoints), source=1.0_pReal)
allocate(damageState(h)%state (1,Nmaterialpoints), source=1.0_pReal)
nullify(damageMapping(h)%p)
damageMapping(h)%p => material_homogenizationMemberAt
deallocate(damage(h)%p)
damage(h)%p => damageState(h)%state(1,:)
end associate

4
src/damage_none.f90
View File

@ -3,6 +3,7 @@
!> @brief material subroutine for constant damage field
!--------------------------------------------------------------------------------------------------
module damage_none
use prec
use config
use material
@ -29,8 +30,7 @@ subroutine damage_none_init
allocate(damageState(h)%subState0(0,Nmaterialpoints))
allocate(damageState(h)%state (0,Nmaterialpoints))
deallocate(damage(h)%p)
allocate (damage(h)%p(1), source=damage_initialPhi(h))
allocate (damage(h)%p(Nmaterialpoints), source=1.0_pReal)
enddo

11
src/damage_nonlocal.f90
View File

@ -78,13 +78,10 @@ subroutine damage_nonlocal_init
Nmaterialpoints = count(material_homogenizationAt == h)
damageState(h)%sizeState = 1
allocate(damageState(h)%state0 (1,Nmaterialpoints), source=damage_initialPhi(h))
allocate(damageState(h)%subState0(1,Nmaterialpoints), source=damage_initialPhi(h))
allocate(damageState(h)%state (1,Nmaterialpoints), source=damage_initialPhi(h))
allocate(damageState(h)%state0 (1,Nmaterialpoints), source=1.0_pReal)
allocate(damageState(h)%subState0(1,Nmaterialpoints), source=1.0_pReal)
allocate(damageState(h)%state (1,Nmaterialpoints), source=1.0_pReal)
nullify(damageMapping(h)%p)
damageMapping(h)%p => material_homogenizationMemberAt
deallocate(damage(h)%p)
damage(h)%p => damageState(h)%state(1,:)
end associate
@ -181,7 +178,7 @@ subroutine damage_nonlocal_putNonLocalDamage(phi,ip,el)
offset
homog = material_homogenizationAt(el)
offset = damageMapping(homog)%p(ip,el)
offset = material_homogenizationMemberAt(ip,el)
damage(homog)%p(offset) = phi
end subroutine damage_nonlocal_putNonLocalDamage

30
src/grid/grid_mech_FEM.f90
View File

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

4
src/grid/grid_mech_spectral_basic.f90
View File

@ -199,7 +199,7 @@ subroutine grid_mech_spectral_basic_init
F = reshape(F_lastInc,[9,grid(1),grid(2),grid3])
endif restartRead
homogenization_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,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
call utilities_updateCoords(reshape(F,shape(F_lastInc)))
call utilities_constitutiveResponse(P,P_av,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2
reshape(F,shape(F_lastInc)), & ! target F
@ -319,7 +319,7 @@ subroutine grid_mech_spectral_basic_forward(cutBack,guess,Delta_t,Delta_t_old,t_
rotation_BC%rotate(F_aimDot,active=.true.))
F_lastInc = reshape(F,[3,3,grid(1),grid(2),grid3])
homogenization_F0 = reshape(F,[3,3,1,product(grid(1:2))*grid3])
homogenization_F0 = reshape(F,[3,3,product(grid(1:2))*grid3])
endif
!--------------------------------------------------------------------------------------------------

6
src/grid/grid_mech_spectral_polarisation.f90
View File

@ -225,7 +225,7 @@ subroutine grid_mech_spectral_polarisation_init
F_tau_lastInc = 2.0_pReal*F_lastInc
endif restartRead
homogenization_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,product(grid(1:2))*grid3]) ! set starting condition for materialpoint_stressAndItsTangent
call utilities_updateCoords(reshape(F,shape(F_lastInc)))
call utilities_constitutiveResponse(P,P_av,C_volAvg,C_minMaxAvg, & ! stress field, stress avg, global average of stiffness and (min+max)/2
reshape(F,shape(F_lastInc)), & ! target F
@ -359,7 +359,7 @@ subroutine grid_mech_spectral_polarisation_forward(cutBack,guess,Delta_t,Delta_t
F_lastInc = reshape(F, [3,3,grid(1),grid(2),grid3])
F_tau_lastInc = reshape(F_tau,[3,3,grid(1),grid(2),grid3])
homogenization_F0 = reshape(F,[3,3,1,product(grid(1:2))*grid3])
homogenization_F0 = reshape(F,[3,3,product(grid(1:2))*grid3])
endif
!--------------------------------------------------------------------------------------------------
@ -604,7 +604,7 @@ subroutine formResidual(in, FandF_tau, &
do k = 1, grid3; do j = 1, grid(2); do i = 1, grid(1)
e = e + 1
residual_F(1:3,1:3,i,j,k) = &
math_mul3333xx33(math_invSym3333(homogenization_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,e) + C_scale), &
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))) &
+ residual_F_tau(1:3,1:3,i,j,k)

3
src/grid/grid_thermal_spectral.f90
View File

@ -131,8 +131,7 @@ subroutine grid_thermal_spectral_init
cell = 0
do k = 1, grid3; do j = 1, grid(2); do i = 1,grid(1)
cell = cell + 1
T_current(i,j,k) = temperature(material_homogenizationAt(cell))% &
p(thermalMapping(material_homogenizationAt(cell))%p(1,cell))
T_current(i,j,k) = temperature(material_homogenizationAt(cell))%p(material_homogenizationMemberAt(1,cell))
T_lastInc(i,j,k) = T_current(i,j,k)
T_stagInc(i,j,k) = T_current(i,j,k)
enddo; enddo; enddo

16
src/grid/spectral_utilities.f90
View File

@ -810,7 +810,7 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
print'(/,a)', ' ... evaluating constitutive response ......................................'
flush(IO_STDOUT)
homogenization_F = reshape(F,[3,3,1,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field
homogenization_F = reshape(F,[3,3,product(grid(1:2))*grid3]) ! set materialpoint target F to estimated field
call materialpoint_stressAndItsTangent(timeinc) ! calculate P field
@ -829,13 +829,13 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
dPdF_min = huge(1.0_pReal)
dPdF_norm_min = huge(1.0_pReal)
do i = 1, product(grid(1:2))*grid3
if (dPdF_norm_max < sum(homogenization_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal)) then
dPdF_max = homogenization_dPdF(1:3,1:3,1:3,1:3,1,i)
dPdF_norm_max = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal)
if (dPdF_norm_max < sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2.0_pReal)) then
dPdF_max = homogenization_dPdF(1:3,1:3,1:3,1:3,i)
dPdF_norm_max = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2.0_pReal)
endif
if (dPdF_norm_min > sum(homogenization_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal)) then
dPdF_min = homogenization_dPdF(1:3,1:3,1:3,1:3,1,i)
dPdF_norm_min = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,1,i)**2.0_pReal)
if (dPdF_norm_min > sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2.0_pReal)) then
dPdF_min = homogenization_dPdF(1:3,1:3,1:3,1:3,i)
dPdF_norm_min = sum(homogenization_dPdF(1:3,1:3,1:3,1:3,i)**2.0_pReal)
endif
end do
@ -853,7 +853,7 @@ subroutine utilities_constitutiveResponse(P,P_av,C_volAvg,C_minmaxAvg,&
C_minmaxAvg = 0.5_pReal*(dPdF_max + dPdF_min)
C_volAvg = sum(sum(homogenization_dPdF,dim=6),dim=5)
C_volAvg = sum(homogenization_dPdF,dim=5)
call MPI_Allreduce(MPI_IN_PLACE,C_volAvg,81,MPI_DOUBLE,MPI_SUM,PETSC_COMM_WORLD,ierr)
if (ierr /= 0) error stop 'MPI error'
C_volAvg = C_volAvg * wgt

231
src/homogenization.f90
View File

@ -30,14 +30,16 @@ module homogenization
!--------------------------------------------------------------------------------------------------
! General variables for the homogenization at a material point
real(pReal), dimension(:,:,:,:), allocatable, public :: &
homogenization_F0, & !< def grad of IP at start of FE increment
homogenization_F !< def grad of IP to be reached at end of FE increment
real(pReal), dimension(:,:,:,:), allocatable, public, protected :: &
homogenization_P !< first P--K stress of IP
real(pReal), dimension(:,:,:,:,:,:), allocatable, public, protected :: &
homogenization_dPdF !< tangent of first P--K stress at IP
real(pReal), dimension(:,:,:), allocatable, public :: &
homogenization_F0, & !< def grad of IP at start of FE increment
homogenization_F !< def grad of IP to be reached at end of FE increment
real(pReal), dimension(:,:,:), allocatable, public :: & !, protected :: & Issue with ifort
homogenization_P !< first P--K stress of IP
real(pReal), dimension(:,:,:,:,:), allocatable, public :: & !, protected :: &
homogenization_dPdF !< tangent of first P--K stress at IP
!--------------------------------------------------------------------------------------------------
type :: tNumerics
integer :: &
nMPstate !< materialpoint state loop limit
@ -62,52 +64,37 @@ module homogenization
type(tDebugOptions) :: debugHomog
!--------------------------------------------------------------------------------------------------
interface
module subroutine mech_none_init
end subroutine mech_none_init
module subroutine mech_isostrain_init
end subroutine mech_isostrain_init
module subroutine mech_RGC_init(num_homogMech)
module subroutine mech_init(num_homog)
class(tNode), pointer, intent(in) :: &
num_homogMech !< pointer to mechanical homogenization numerics data
end subroutine mech_RGC_init
num_homog !< pointer to mechanical homogenization numerics data
end subroutine mech_init
module subroutine mech_partition(subF,ip,el)
real(pReal), intent(in), dimension(3,3) :: &
subF
integer, intent(in) :: &
ip, & !< integration point
el !< element number
end subroutine mech_partition
module subroutine mech_isostrain_partitionDeformation(F,avgF)
real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned deformation gradient
real(pReal), dimension (3,3), intent(in) :: avgF !< average deformation gradient at material point
end subroutine mech_isostrain_partitionDeformation
module subroutine mech_homogenize(ip,el)
integer, intent(in) :: &
ip, & !< integration point
el !< element number
end subroutine mech_homogenize
module subroutine mech_RGC_partitionDeformation(F,avgF,instance,of)
real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned deformation gradient
real(pReal), dimension (3,3), intent(in) :: avgF !< average deformation gradient at material point
integer, intent(in) :: &
instance, &
of
end subroutine mech_RGC_partitionDeformation
module subroutine mech_results(group_base,h)
character(len=*), intent(in) :: group_base
integer, intent(in) :: h
module subroutine mech_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance)
real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point
real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point
real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses
real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses
integer, intent(in) :: instance
end subroutine mech_isostrain_averageStressAndItsTangent
module subroutine mech_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance)
real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point
real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point
real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses
real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses
integer, intent(in) :: instance
end subroutine mech_RGC_averageStressAndItsTangent
end subroutine mech_results
! -------- ToDo ---------------------------------------------------------
module function mech_RGC_updateState(P,F,F0,avgF,dt,dPdF,ip,el)
logical, dimension(2) :: mech_RGC_updateState
real(pReal), dimension(:,:,:), intent(in) :: &
@ -122,13 +109,8 @@ module homogenization
el !< element number
end function mech_RGC_updateState
module subroutine mech_RGC_results(instance,group)
integer, intent(in) :: instance !< homogenization instance
character(len=*), intent(in) :: group !< group name in HDF5 file
end subroutine mech_RGC_results
end interface
! -----------------------------------------------------------------------
public :: &
homogenization_init, &
@ -145,10 +127,11 @@ subroutine homogenization_init
class (tNode) , pointer :: &
num_homog, &
num_homogMech, &
num_homogGeneric, &
debug_homogenization
print'(/,a)', ' <<<+- homogenization init -+>>>'; flush(IO_STDOUT)
debug_homogenization => config_debug%get('homogenization', defaultVal=emptyList)
debugHomog%basic = debug_homogenization%contains('basic')
debugHomog%extensive = debug_homogenization%contains('extensive')
@ -163,31 +146,8 @@ subroutine homogenization_init
num_homog => config_numerics%get('homogenization',defaultVal=emptyDict)
num_homogMech => num_homog%get('mech',defaultVal=emptyDict)
num_homogGeneric => num_homog%get('generic',defaultVal=emptyDict)
if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) call mech_none_init
if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) call mech_isostrain_init
if (any(homogenization_type == HOMOGENIZATION_RGC_ID)) call mech_RGC_init(num_homogMech)
if (any(thermal_type == THERMAL_isothermal_ID)) call thermal_isothermal_init
if (any(thermal_type == THERMAL_adiabatic_ID)) call thermal_adiabatic_init
if (any(thermal_type == THERMAL_conduction_ID)) call thermal_conduction_init
if (any(damage_type == DAMAGE_none_ID)) call damage_none_init
if (any(damage_type == DAMAGE_local_ID)) call damage_local_init
if (any(damage_type == DAMAGE_nonlocal_ID)) call damage_nonlocal_init
!--------------------------------------------------------------------------------------------------
! allocate and initialize global variables
allocate(homogenization_dPdF(3,3,3,3,discretization_nIPs,discretization_Nelems), source=0.0_pReal)
homogenization_F0 = spread(spread(math_I3,3,discretization_nIPs),4,discretization_Nelems) ! initialize to identity
homogenization_F = homogenization_F0 ! initialize to identity
allocate(homogenization_P(3,3,discretization_nIPs,discretization_Nelems), source=0.0_pReal)
print'(/,a)', ' <<<+- homogenization init -+>>>'; flush(IO_STDOUT)
num%nMPstate = num_homogGeneric%get_asInt ('nMPstate', defaultVal=10)
num%subStepMinHomog = num_homogGeneric%get_asFloat('subStepMin', defaultVal=1.0e-3_pReal)
num%subStepSizeHomog = num_homogGeneric%get_asFloat('subStepSize', defaultVal=0.25_pReal)
@ -198,6 +158,18 @@ subroutine homogenization_init
if (num%subStepSizeHomog <= 0.0_pReal) call IO_error(301,ext_msg='subStepSizeHomog')
if (num%stepIncreaseHomog <= 0.0_pReal) call IO_error(301,ext_msg='stepIncreaseHomog')
call mech_init(num_homog)
if (any(thermal_type == THERMAL_isothermal_ID)) call thermal_isothermal_init
if (any(thermal_type == THERMAL_adiabatic_ID)) call thermal_adiabatic_init
if (any(thermal_type == THERMAL_conduction_ID)) call thermal_conduction_init
if (any(damage_type == DAMAGE_none_ID)) call damage_none_init
if (any(damage_type == DAMAGE_local_ID)) call damage_local_init
if (any(damage_type == DAMAGE_nonlocal_ID)) call damage_nonlocal_init
end subroutine homogenization_init
@ -221,6 +193,7 @@ subroutine materialpoint_stressAndItsTangent(dt)
converged
logical, dimension(2,discretization_nIPs,discretization_Nelems) :: &
doneAndHappy
integer :: m
!--------------------------------------------------------------------------------------------------
@ -255,7 +228,7 @@ subroutine materialpoint_stressAndItsTangent(dt)
any(subStep(FEsolving_execIP(1):FEsolving_execIP(2),&
FEsolving_execElem(1):FEsolving_execElem(2)) > num%subStepMinHomog))
!$OMP PARALLEL DO
!$OMP PARALLEL DO PRIVATE(m)
elementLooping1: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Nconstituents(material_homogenizationAt(e))
IpLooping1: do i = FEsolving_execIP(1),FEsolving_execIP(2)
@ -325,13 +298,14 @@ subroutine materialpoint_stressAndItsTangent(dt)
!--------------------------------------------------------------------------------------------------
! deformation partitioning
!$OMP PARALLEL DO PRIVATE(myNgrains)
!$OMP PARALLEL DO PRIVATE(myNgrains,m)
elementLooping2: do e = FEsolving_execElem(1),FEsolving_execElem(2)
myNgrains = homogenization_Nconstituents(material_homogenizationAt(e))
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
call partitionDeformation(homogenization_F0(1:3,1:3,i,e) &
+ (homogenization_F(1:3,1:3,i,e)-homogenization_F0(1:3,1:3,i,e))&
m = (e-1)*discretization_nIPs + i
call mech_partition(homogenization_F0(1:3,1:3,m) &
+ (homogenization_F(1:3,1:3,m)-homogenization_F0(1:3,1:3,m))&
*(subStep(i,e)+subFrac(i,e)), &
i,e)
crystallite_dt(1:myNgrains,i,e) = dt*subStep(i,e) ! propagate materialpoint dt to grains
@ -349,16 +323,17 @@ subroutine materialpoint_stressAndItsTangent(dt)
!--------------------------------------------------------------------------------------------------
! state update
!$OMP PARALLEL DO
!$OMP PARALLEL DO PRIVATE(m)
elementLooping3: do e = FEsolving_execElem(1),FEsolving_execElem(2)
IpLooping3: do i = FEsolving_execIP(1),FEsolving_execIP(2)
if (requested(i,e) .and. .not. doneAndHappy(1,i,e)) then
if (.not. converged(i,e)) then
doneAndHappy(1:2,i,e) = [.true.,.false.]
else
m = (e-1)*discretization_nIPs + i
doneAndHappy(1:2,i,e) = updateState(dt*subStep(i,e), &
homogenization_F0(1:3,1:3,i,e) &
+ (homogenization_F(1:3,1:3,i,e)-homogenization_F0(1:3,1:3,i,e)) &
homogenization_F0(1:3,1:3,m) &
+ (homogenization_F(1:3,1:3,m)-homogenization_F0(1:3,1:3,m)) &
*(subStep(i,e)+subFrac(i,e)), &
i,e)
converged(i,e) = all(doneAndHappy(1:2,i,e)) ! converged if done and happy
@ -379,7 +354,7 @@ subroutine materialpoint_stressAndItsTangent(dt)
!$OMP PARALLEL DO
elementLooping4: do e = FEsolving_execElem(1),FEsolving_execElem(2)
IpLooping4: do i = FEsolving_execIP(1),FEsolving_execIP(2)
call averageStressAndItsTangent(i,e)
call mech_homogenize(i,e)
enddo IpLooping4
enddo elementLooping4
!$OMP END PARALLEL DO
@ -390,38 +365,6 @@ subroutine materialpoint_stressAndItsTangent(dt)
end subroutine materialpoint_stressAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief partition material point def grad onto constituents
!--------------------------------------------------------------------------------------------------
subroutine partitionDeformation(subF,ip,el)
real(pReal), intent(in), dimension(3,3) :: &
subF
integer, intent(in) :: &
ip, & !< integration point
el !< element number
chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
case (HOMOGENIZATION_NONE_ID) chosenHomogenization
crystallite_partitionedF(1:3,1:3,1,ip,el) = subF
case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
call mech_isostrain_partitionDeformation(&
crystallite_partitionedF(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
subF)
case (HOMOGENIZATION_RGC_ID) chosenHomogenization
call mech_RGC_partitionDeformation(&
crystallite_partitionedF(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
subF,&
ip, &
el)
end select chosenHomogenization
end subroutine partitionDeformation
!--------------------------------------------------------------------------------------------------
!> @brief update the internal state of the homogenization scheme and tell whether "done" and
!> "happy" with result
@ -478,49 +421,6 @@ function updateState(subdt,subF,ip,el)
end function updateState
!--------------------------------------------------------------------------------------------------
!> @brief derive average stress and stiffness from constituent quantities
!--------------------------------------------------------------------------------------------------
subroutine averageStressAndItsTangent(ip,el)
integer, intent(in) :: &
ip, & !< integration point
el !< element number
integer :: c
real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
case (HOMOGENIZATION_NONE_ID) chosenHomogenization
homogenization_P(1:3,1:3,ip,el) = crystallite_P(1:3,1:3,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
do c = 1, homogenization_Nconstituents(material_homogenizationAt(el))
dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el)
enddo
call mech_isostrain_averageStressAndItsTangent(&
homogenization_P(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_Nconstituents(material_homogenizationAt(el)),ip,el), &
dPdFs, &
homogenization_typeInstance(material_homogenizationAt(el)))
case (HOMOGENIZATION_RGC_ID) chosenHomogenization
do c = 1, homogenization_Nconstituents(material_homogenizationAt(el))
dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el)
enddo
call mech_RGC_averageStressAndItsTangent(&
homogenization_P(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_Nconstituents(material_homogenizationAt(el)),ip,el), &
dPdFs, &
homogenization_typeInstance(material_homogenizationAt(el)))
end select chosenHomogenization
end subroutine averageStressAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief writes homogenization results to HDF5 output file
!--------------------------------------------------------------------------------------------------
@ -531,27 +431,12 @@ subroutine homogenization_results
integer :: p
character(len=:), allocatable :: group_base,group
!real(pReal), dimension(:,:,:), allocatable :: temp
do p=1,size(material_name_homogenization)
group_base = 'current/homogenization/'//trim(material_name_homogenization(p))
call results_closeGroup(results_addGroup(group_base))
group = trim(group_base)//'/generic'
call results_closeGroup(results_addGroup(group))
!temp = reshape(homogenization_F,[3,3,discretization_nIPs*discretization_Nelems])
!call results_writeDataset(group,temp,'F',&
! 'deformation gradient','1')
!temp = reshape(homogenization_P,[3,3,discretization_nIPs*discretization_Nelems])
!call results_writeDataset(group,temp,'P',&
! '1st Piola-Kirchhoff stress','Pa')
group = trim(group_base)//'/mech'
call results_closeGroup(results_addGroup(group))
select case(material_homogenization_type(p))
case(HOMOGENIZATION_rgc_ID)
call mech_RGC_results(homogenization_typeInstance(p),group)
end select
call mech_results(group_base,p)
group = trim(group_base)//'/damage'
call results_closeGroup(results_addGroup(group))

200
src/homogenization_mech.f90 Normal file
View File

@ -0,0 +1,200 @@
!--------------------------------------------------------------------------------------------------
!> @author Martin Diehl, KU Leuven
!> @brief Partition F and homogenize P/dPdF
!--------------------------------------------------------------------------------------------------
submodule(homogenization) homogenization_mech
interface
module subroutine mech_none_init
end subroutine mech_none_init
module subroutine mech_isostrain_init
end subroutine mech_isostrain_init
module subroutine mech_RGC_init(num_homogMech)
class(tNode), pointer, intent(in) :: &
num_homogMech !< pointer to mechanical homogenization numerics data
end subroutine mech_RGC_init
module subroutine mech_isostrain_partitionDeformation(F,avgF)
real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned deformation gradient
real(pReal), dimension (3,3), intent(in) :: avgF !< average deformation gradient at material point
end subroutine mech_isostrain_partitionDeformation
module subroutine mech_RGC_partitionDeformation(F,avgF,instance,of)
real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned deformation gradient
real(pReal), dimension (3,3), intent(in) :: avgF !< average deformation gradient at material point
integer, intent(in) :: &
instance, &
of
end subroutine mech_RGC_partitionDeformation
module subroutine mech_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance)
real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point
real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point
real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses
real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses
integer, intent(in) :: instance
end subroutine mech_isostrain_averageStressAndItsTangent
module subroutine mech_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,instance)
real(pReal), dimension (3,3), intent(out) :: avgP !< average stress at material point
real(pReal), dimension (3,3,3,3), intent(out) :: dAvgPdAvgF !< average stiffness at material point
real(pReal), dimension (:,:,:), intent(in) :: P !< partitioned stresses
real(pReal), dimension (:,:,:,:,:), intent(in) :: dPdF !< partitioned stiffnesses
integer, intent(in) :: instance
end subroutine mech_RGC_averageStressAndItsTangent
module subroutine mech_RGC_results(instance,group)
integer, intent(in) :: instance !< homogenization instance
character(len=*), intent(in) :: group !< group name in HDF5 file
end subroutine mech_RGC_results
end interface
contains
!--------------------------------------------------------------------------------------------------
!> @brief Allocate variables and set parameters.
!--------------------------------------------------------------------------------------------------
module subroutine mech_init(num_homog)
class(tNode), pointer, intent(in) :: &
num_homog
class(tNode), pointer :: &
num_homogMech
print'(/,a)', ' <<<+- homogenization_mech init -+>>>'
allocate(homogenization_dPdF(3,3,3,3,discretization_nIPs*discretization_Nelems), source=0.0_pReal)
homogenization_F0 = spread(math_I3,3,discretization_nIPs*discretization_Nelems) ! initialize to identity
homogenization_F = homogenization_F0 ! initialize to identity
allocate(homogenization_P(3,3,discretization_nIPs*discretization_Nelems), source=0.0_pReal)
num_homogMech => num_homog%get('mech',defaultVal=emptyDict)
if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) call mech_none_init
if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) call mech_isostrain_init
if (any(homogenization_type == HOMOGENIZATION_RGC_ID)) call mech_RGC_init(num_homogMech)
end subroutine mech_init
!--------------------------------------------------------------------------------------------------
!> @brief Partition F onto the individual constituents.
!--------------------------------------------------------------------------------------------------
module subroutine mech_partition(subF,ip,el)
real(pReal), intent(in), dimension(3,3) :: &
subF
integer, intent(in) :: &
ip, & !< integration point
el !< element number
chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
case (HOMOGENIZATION_NONE_ID) chosenHomogenization
crystallite_partitionedF(1:3,1:3,1,ip,el) = subF
case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
call mech_isostrain_partitionDeformation(&
crystallite_partitionedF(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
subF)
case (HOMOGENIZATION_RGC_ID) chosenHomogenization
call mech_RGC_partitionDeformation(&
crystallite_partitionedF(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
subF,&
ip, &
el)
end select chosenHomogenization
end subroutine mech_partition
!--------------------------------------------------------------------------------------------------
!> @brief Average P and dPdF from the individual constituents.
!--------------------------------------------------------------------------------------------------
module subroutine mech_homogenize(ip,el)
integer, intent(in) :: &
ip, & !< integration point
el !< element number
integer :: c,m
real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt(el)))
m = (el-1)* discretization_nIPs + ip
chosenHomogenization: select case(homogenization_type(material_homogenizationAt(el)))
case (HOMOGENIZATION_NONE_ID) chosenHomogenization
homogenization_P(1:3,1:3,m) = crystallite_P(1:3,1:3,1,ip,el)
homogenization_dPdF(1:3,1:3,1:3,1:3,m) = crystallite_stressTangent(1,ip,el)
case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization
do c = 1, homogenization_Nconstituents(material_homogenizationAt(el))
dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el)
enddo
call mech_isostrain_averageStressAndItsTangent(&
homogenization_P(1:3,1:3,m), &
homogenization_dPdF(1:3,1:3,1:3,1:3,m),&
crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
dPdFs, &
homogenization_typeInstance(material_homogenizationAt(el)))
case (HOMOGENIZATION_RGC_ID) chosenHomogenization
do c = 1, homogenization_Nconstituents(material_homogenizationAt(el))
dPdFs(:,:,:,:,c) = crystallite_stressTangent(c,ip,el)
enddo
call mech_RGC_averageStressAndItsTangent(&
homogenization_P(1:3,1:3,m), &
homogenization_dPdF(1:3,1:3,1:3,1:3,m),&
crystallite_P(1:3,1:3,1:homogenization_Nconstituents(material_homogenizationAt(el)),ip,el), &
dPdFs, &
homogenization_typeInstance(material_homogenizationAt(el)))
end select chosenHomogenization
end subroutine mech_homogenize
!--------------------------------------------------------------------------------------------------
!> @brief Write results to file.
!--------------------------------------------------------------------------------------------------
module subroutine mech_results(group_base,h)
use material, only: &
material_homogenization_type => homogenization_type
character(len=*), intent(in) :: group_base
integer, intent(in) :: h
character(len=:), allocatable :: group
group = trim(group_base)//'/mech'
call results_closeGroup(results_addGroup(group))
select case(material_homogenization_type(h))
case(HOMOGENIZATION_rgc_ID)
call mech_RGC_results(homogenization_typeInstance(h),group)
end select
!temp = reshape(homogenization_F,[3,3,discretization_nIPs*discretization_Nelems])
!call results_writeDataset(group,temp,'F',&
! 'deformation gradient','1')
!temp = reshape(homogenization_P,[3,3,discretization_nIPs*discretization_Nelems])
!call results_writeDataset(group,temp,'P',&
! '1st Piola-Kirchhoff stress','Pa')
end subroutine mech_results
end submodule homogenization_mech

2
src/homogenization_mech_RGC.f90
View File

@ -6,7 +6,7 @@
!> @brief Relaxed grain cluster (RGC) homogenization scheme
!> N_constituents is defined as p x q x r (cluster)
!--------------------------------------------------------------------------------------------------
submodule(homogenization) homogenization_mech_RGC
submodule(homogenization:homogenization_mech) homogenization_mech_RGC
use rotations
type :: tParameters

2
src/homogenization_mech_isostrain.f90
View File

@ -4,7 +4,7 @@
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Isostrain (full constraint Taylor assuption) homogenization scheme
!--------------------------------------------------------------------------------------------------
submodule(homogenization) homogenization_mech_isostrain
submodule(homogenization:homogenization_mech) homogenization_mech_isostrain
enum, bind(c); enumerator :: &
parallel_ID, &

2
src/homogenization_mech_none.f90
View File

@ -4,7 +4,7 @@
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
!> @brief dummy homogenization homogenization scheme for 1 constituent per material point
!--------------------------------------------------------------------------------------------------
submodule(homogenization) homogenization_mech_none
submodule(homogenization:homogenization_mech) homogenization_mech_none
contains

2
src/kinematics_cleavage_opening.f90
View File

@ -120,7 +120,7 @@ module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S,
udotd, dudotd_dt, udott, dudott_dt, udotn, dudotn_dt
homog = material_homogenizationAt(el)
damageOffset = damageMapping(homog)%p(ip,el)
damageOffset = material_homogenizationMemberAt(ip,el)
Ld = 0.0_pReal
dLd_dTstar = 0.0_pReal

2
src/kinematics_slipplane_opening.f90
View File

@ -141,7 +141,7 @@ module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S
phase = material_phaseAt(ipc,el)
instance = kinematics_slipplane_opening_instance(phase)
homog = material_homogenizationAt(el)
damageOffset = damageMapping(homog)%p(ip,el)
damageOffset = material_homogenizationMemberAt(ip,el)
associate(prm => param(instance))
Ld = 0.0_pReal

4
src/kinematics_thermal_expansion.f90
View File

@ -126,8 +126,8 @@ module subroutine kinematics_thermal_expansion_LiAndItsTangent(Li, dLi_dTstar, i
phase = material_phaseAt(ipc,el)
homog = material_homogenizationAt(el)
T = temperature(homog)%p(thermalMapping(homog)%p(ip,el))
TDot = temperatureRate(homog)%p(thermalMapping(homog)%p(ip,el))
T = temperature(homog)%p(material_homogenizationMemberAt(ip,el))
TDot = temperatureRate(homog)%p(material_homogenizationMemberAt(ip,el))
associate(prm => param(kinematics_thermal_expansion_instance(phase)))
Li = TDot * ( &

38
src/material.f90
View File

@ -64,21 +64,20 @@ module material
homogenization_type !< type of each homogenization
integer, public, protected :: &
homogenization_maxNconstituents !< max number of grains in any USED homogenization
homogenization_maxNconstituents !< max number of grains in any USED homogenization
integer, dimension(:), allocatable, public, protected :: &
homogenization_Nconstituents, & !< number of grains in each homogenization
homogenization_Nconstituents, & !< number of grains in each homogenization
homogenization_typeInstance, & !< instance of particular type of each homogenization
thermal_typeInstance, & !< instance of particular type of each thermal transport
damage_typeInstance !< instance of particular type of each nonlocal damage
real(pReal), dimension(:), allocatable, public, protected :: &
thermal_initialT, & !< initial temperature per each homogenization
damage_initialPhi !< initial damage per each homogenization
thermal_initialT !< initial temperature per each homogenization
integer, dimension(:), allocatable, public, protected :: & ! (elem)
material_homogenizationAt !< homogenization ID of each element
integer, dimension(:,:), allocatable, public, target :: & ! (ip,elem) ToDo: ugly target for mapping hack
integer, dimension(:,:), allocatable, public, protected :: & ! (ip,elem)
material_homogenizationMemberAt !< position of the element within its homogenization instance
integer, dimension(:,:), allocatable, public, protected :: & ! (constituent,elem)
material_phaseAt !< phase ID of each element
@ -93,14 +92,6 @@ module material
type(Rotation), dimension(:,:,:), allocatable, public, protected :: &
material_orientation0 !< initial orientation of each grain,IP,element
! BEGIN DEPRECATED
integer, dimension(:,:), allocatable, private, target :: mappingHomogenizationConst !< mapping from material points to offset in constant state/field
! END DEPRECATED
type(tHomogMapping), allocatable, dimension(:), public :: &
thermalMapping, & !< mapping for thermal state/fields
damageMapping !< mapping for damage state/fields
type(group_float), allocatable, dimension(:), public :: &
temperature, & !< temperature field
damage, & !< damage field
@ -149,7 +140,6 @@ contains
subroutine material_init(restart)
logical, intent(in) :: restart
integer :: myHomog
print'(/,a)', ' <<<+- material init -+>>>'; flush(IO_STDOUT)
@ -165,12 +155,8 @@ subroutine material_init(restart)
allocate(thermalState (size(material_name_homogenization)))
allocate(damageState (size(material_name_homogenization)))
allocate(thermalMapping (size(material_name_homogenization)))
allocate(damageMapping (size(material_name_homogenization)))
allocate(temperature (size(material_name_homogenization)))
allocate(damage (size(material_name_homogenization)))
allocate(temperatureRate (size(material_name_homogenization)))
@ -181,20 +167,6 @@ subroutine material_init(restart)
call results_closeJobFile
endif
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! BEGIN DEPRECATED
allocate(mappingHomogenizationConst( discretization_nIPs,discretization_Nelems),source=1)
! hack needed to initialize field values used during constitutive initialization
do myHomog = 1, size(material_name_homogenization)
thermalMapping (myHomog)%p => mappingHomogenizationConst
damageMapping (myHomog)%p => mappingHomogenizationConst
allocate(temperature (myHomog)%p(1), source=thermal_initialT(myHomog))
allocate(damage (myHomog)%p(1), source=damage_initialPhi(myHomog))
allocate(temperatureRate (myHomog)%p(1), source=0.0_pReal)
enddo
! END DEPRECATED
end subroutine material_init
@ -222,7 +194,6 @@ subroutine material_parseHomogenization
allocate(thermal_typeInstance(size(material_name_homogenization)), source=0)
allocate(damage_typeInstance(size(material_name_homogenization)), source=0)
allocate(thermal_initialT(size(material_name_homogenization)), source=300.0_pReal)
allocate(damage_initialPhi(size(material_name_homogenization)), source=1.0_pReal)
do h=1, size(material_name_homogenization)
homog => material_homogenization%get(h)
@ -258,7 +229,6 @@ subroutine material_parseHomogenization
if(homog%contains('damage')) then
homogDamage => homog%get('damage')
damage_initialPhi(h) = homogDamage%get_asFloat('phi_0',defaultVal=1.0_pReal)
select case (homogDamage%get_asString('type'))
case('none')
damage_type(h) = DAMAGE_none_ID

2
src/mesh/FEM_utilities.f90
View File

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

43
src/mesh/mesh_mech_FEM.f90
View File

@ -316,16 +316,16 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
Vec :: x_local, f_local, xx_local
PetscSection :: section
PetscScalar, dimension(:), pointer :: x_scal, pf_scal
PetscScalar, target :: f_scal(cellDof)
PetscReal :: detJ, IcellJMat(dimPlex,dimPlex)
PetscReal, pointer,dimension(:) :: pV0, pCellJ, pInvcellJ, basisField, basisFieldDer
PetscScalar, dimension(cellDof), target :: f_scal
PetscReal :: IcellJMat(dimPlex,dimPlex)
PetscReal, dimension(:),pointer :: pV0, pCellJ, pInvcellJ, basisField, basisFieldDer
PetscInt :: cellStart, cellEnd, cell, field, face, &
qPt, basis, comp, cidx, &
numFields
PetscReal :: detFAvg
PetscReal :: BMat(dimPlex*dimPlex,cellDof)
numFields, &
bcSize,m
PetscReal :: detFAvg, detJ
PetscReal, dimension(dimPlex*dimPlex,cellDof) :: BMat
PetscInt :: bcSize
IS :: bcPoints
@ -366,6 +366,7 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
CHKERRQ(ierr)
IcellJMat = reshape(pInvcellJ,shape=[dimPlex,dimPlex])
do qPt = 0, nQuadrature-1
m = cell*nQuadrature + qPt+1
BMat = 0.0
do basis = 0, nBasis-1
do comp = 0, dimPlex-1
@ -375,15 +376,14 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
(((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp+1)*dimPlex))
enddo
enddo
homogenization_F(1:dimPlex,1:dimPlex,qPt+1,cell+1) = &
reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex], order=[2,1])
homogenization_F(1:dimPlex,1:dimPlex,m) = reshape(matmul(BMat,x_scal),shape=[dimPlex,dimPlex], order=[2,1])
enddo
if (num%BBarStabilisation) then
detFAvg = math_det33(sum(homogenization_F(1:3,1:3,1:nQuadrature,cell+1),dim=3)/real(nQuadrature))
do qPt = 1, nQuadrature
homogenization_F(1:dimPlex,1:dimPlex,qPt,cell+1) = &
homogenization_F(1:dimPlex,1:dimPlex,qPt,cell+1)* &
(detFAvg/math_det33(homogenization_F(1:3,1:3,qPt,cell+1)))**(1.0/real(dimPlex))
detFAvg = math_det33(sum(homogenization_F(1:3,1:3,cell*nQuadrature+1:(cell+1)*nQuadrature),dim=3)/real(nQuadrature))
do qPt = 0, nQuadrature-1
m = cell*nQuadrature + qPt+1
homogenization_F(1:dimPlex,1:dimPlex,m) = homogenization_F(1:dimPlex,1:dimPlex,m) &
* (detFAvg/math_det33(homogenization_F(1:3,1:3,m)))**(1.0/real(dimPlex))
enddo
endif
@ -407,6 +407,7 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
IcellJMat = reshape(pInvcellJ,shape=[dimPlex,dimPlex])
f_scal = 0.0
do qPt = 0, nQuadrature-1
m = cell*nQuadrature + qPt+1
BMat = 0.0
do basis = 0, nBasis-1
do comp = 0, dimPlex-1
@ -418,7 +419,7 @@ subroutine FEM_mech_formResidual(dm_local,xx_local,f_local,dummy,ierr)
enddo
f_scal = f_scal + &
matmul(transpose(BMat), &
reshape(transpose(homogenization_P(1:dimPlex,1:dimPlex,qPt+1,cell+1)), &
reshape(transpose(homogenization_P(1:dimPlex,1:dimPlex,m)), &
shape=[dimPlex*dimPlex]))*qWeights(qPt+1)
enddo
f_scal = f_scal*abs(detJ)
@ -463,7 +464,7 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
K_eB
PetscInt :: cellStart, cellEnd, cell, field, face, &
qPt, basis, comp, cidx,bcSize
qPt, basis, comp, cidx,bcSize, m
IS :: bcPoints
@ -506,6 +507,7 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
FAvg = 0.0
BMatAvg = 0.0
do qPt = 0, nQuadrature-1
m = cell*nQuadrature + qPt + 1
BMat = 0.0
do basis = 0, nBasis-1
do comp = 0, dimPlex-1
@ -516,7 +518,7 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
(((qPt*nBasis + basis)*dimPlex + comp)*dimPlex+comp+1)*dimPlex))
enddo
enddo
MatA = matmul(reshape(reshape(homogenization_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,m), &
shape=[dimPlex,dimPlex,dimPlex,dimPlex], order=[2,1,4,3]), &
shape=[dimPlex*dimPlex,dimPlex*dimPlex]),BMat)*qWeights(qPt+1)
if (num%BBarStabilisation) then
@ -524,12 +526,11 @@ subroutine FEM_mech_formJacobian(dm_local,xx_local,Jac_pre,Jac,dummy,ierr)
FInv = math_inv33(F)
K_eA = K_eA + matmul(transpose(BMat),MatA)*math_det33(FInv)**(1.0/real(dimPlex))
K_eB = K_eB - &
matmul(transpose(matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,qPt+1,cell+1), &
shape=[dimPlex*dimPlex,1]), &
matmul(transpose(matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,m),shape=[dimPlex*dimPlex,1]), &
matmul(reshape(FInv(1:dimPlex,1:dimPlex), &
shape=[1,dimPlex*dimPlex],order=[2,1]),BMat))),MatA)
MatB = MatB + &
matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,qPt+1,cell+1),shape=[1,dimPlex*dimPlex]),MatA)
MatB = MatB &
+ matmul(reshape(homogenization_F(1:dimPlex,1:dimPlex,m),shape=[1,dimPlex*dimPlex]),MatA)
FAvg = FAvg + F
BMatAvg = BMatAvg + BMat
else

6
src/prec.f90
View File

@ -54,7 +54,7 @@ module prec
type, extends(tState) :: tPlasticState
logical :: &
nonlocal = .false.
real(pReal), pointer, dimension(:,:) :: &
real(pReal), pointer, dimension(:,:) :: &
slipRate !< slip rate
end type
@ -62,10 +62,6 @@ module prec
type(tState), dimension(:), allocatable :: p !< tState for each active source mechanism in a phase
end type
type :: tHomogMapping
integer, pointer, dimension(:,:) :: p
end type
real(pReal), private, parameter :: PREAL_EPSILON = epsilon(0.0_pReal) !< minimum positive number such that 1.0 + EPSILON /= 1.0.
real(pReal), private, parameter :: PREAL_MIN = tiny(0.0_pReal) !< smallest normalized floating point number

2
src/source_damage_anisoBrittle.f90
View File

@ -143,7 +143,7 @@ module subroutine source_damage_anisoBrittle_dotState(S, ipc, ip, el)
constituent = material_phasememberAt(ipc,ip,el)
sourceOffset = source_damage_anisoBrittle_offset(phase)
homog = material_homogenizationAt(el)
damageOffset = damageMapping(homog)%p(ip,el)
damageOffset = material_homogenizationMemberAt(ip,el)
associate(prm => param(source_damage_anisoBrittle_instance(phase)))
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = 0.0_pReal

2
src/source_damage_anisoDuctile.f90
View File

@ -125,7 +125,7 @@ module subroutine source_damage_anisoDuctile_dotState(ipc, ip, el)
constituent = material_phasememberAt(ipc,ip,el)
sourceOffset = source_damage_anisoDuctile_offset(phase)
homog = material_homogenizationAt(el)
damageOffset = damageMapping(homog)%p(ip,el)
damageOffset = material_homogenizationMemberAt(ip,el)
associate(prm => param(source_damage_anisoDuctile_instance(phase)))
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) &

2
src/source_damage_isoDuctile.f90
View File

@ -116,7 +116,7 @@ module subroutine source_damage_isoDuctile_dotState(ipc, ip, el)
constituent = material_phasememberAt(ipc,ip,el)
sourceOffset = source_damage_isoDuctile_offset(phase)
homog = material_homogenizationAt(el)
damageOffset = damageMapping(homog)%p(ip,el)
damageOffset = material_homogenizationMemberAt(ip,el)
associate(prm => param(source_damage_isoDuctile_instance(phase)))
sourceState(phase)%p(sourceOffset)%dotState(1,constituent) = &

9
src/thermal_adiabatic.f90
View File

@ -73,11 +73,8 @@ subroutine thermal_adiabatic_init
allocate(thermalState(h)%subState0(1,Nmaterialpoints), source=thermal_initialT(h))
allocate(thermalState(h)%state (1,Nmaterialpoints), source=thermal_initialT(h))
thermalMapping(h)%p => material_homogenizationMemberAt
deallocate(temperature(h)%p)
temperature(h)%p => thermalState(h)%state(1,:)
deallocate(temperatureRate(h)%p)
allocate (temperatureRate(h)%p(Nmaterialpoints), source=0.0_pReal)
allocate(temperatureRate(h)%p(Nmaterialpoints),source = 0.0_pReal)
end associate
enddo
@ -117,8 +114,8 @@ function thermal_adiabatic_updateState(subdt, ip, el)
<= 1.0e-6_pReal*abs(thermalState(homog)%state(1,offset)), &
.true.]
temperature (homog)%p(thermalMapping(homog)%p(ip,el)) = T
temperatureRate(homog)%p(thermalMapping(homog)%p(ip,el)) = &
temperature (homog)%p(material_homogenizationMemberAt(ip,el)) = T
temperatureRate(homog)%p(material_homogenizationMemberAt(ip,el)) = &
(thermalState(homog)%state(1,offset) - thermalState(homog)%subState0(1,offset))/(subdt+tiny(0.0_pReal))
end function thermal_adiabatic_updateState

5
src/thermal_conduction.f90
View File

@ -71,10 +71,7 @@ subroutine thermal_conduction_init
allocate(thermalState(h)%subState0(0,Nmaterialpoints))
allocate(thermalState(h)%state (0,Nmaterialpoints))
thermalMapping(h)%p => material_homogenizationMemberAt
deallocate(temperature (h)%p)
allocate (temperature (h)%p(Nmaterialpoints), source=thermal_initialT(h))
deallocate(temperatureRate(h)%p)
allocate (temperatureRate(h)%p(Nmaterialpoints), source=0.0_pReal)
end associate
@ -205,7 +202,7 @@ subroutine thermal_conduction_putTemperatureAndItsRate(T,Tdot,ip,el)
offset
homog = material_homogenizationAt(el)
offset = thermalMapping(homog)%p(ip,el)
offset = material_homogenizationMemberAt(ip,el)
temperature (homog)%p(offset) = T
temperatureRate(homog)%p(offset) = Tdot

7
src/thermal_isothermal.f90
View File

@ -3,6 +3,7 @@
!> @brief material subroutine for isothermal temperature field
!--------------------------------------------------------------------------------------------------
module thermal_isothermal
use prec
use config
use material
@ -29,10 +30,8 @@ subroutine thermal_isothermal_init
allocate(thermalState(h)%subState0(0,Nmaterialpoints))
allocate(thermalState(h)%state (0,Nmaterialpoints))
deallocate(temperature (h)%p)
allocate (temperature (h)%p(1), source=thermal_initialT(h))
deallocate(temperatureRate(h)%p)
allocate (temperatureRate(h)%p(1))
allocate(temperature (h)%p(Nmaterialpoints),source=thermal_initialT(h))
allocate(temperatureRate(h)%p(Nmaterialpoints),source = 0.0_pReal)
enddo