diff --git a/src/CPFEM.f90 b/src/CPFEM.f90 index c49ecbcb6..a19a70432 100644 --- a/src/CPFEM.f90 +++ b/src/CPFEM.f90 @@ -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 diff --git a/src/CPFEM2.f90 b/src/CPFEM2.f90 index 994859758..54e381d34 100644 --- a/src/CPFEM2.f90 +++ b/src/CPFEM2.f90 @@ -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 diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index d161b36eb..beabfcae1 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -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" diff --git a/src/constitutive.f90 b/src/constitutive.f90 index 470ac4dc7..a60352f81 100644 --- a/src/constitutive.f90 +++ b/src/constitutive.f90 @@ -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) diff --git a/src/constitutive_mech.f90 b/src/constitutive_mech.f90 index dc3a935e3..6b3c6fce6 100644 --- a/src/constitutive_mech.f90 +++ b/src/constitutive_mech.f90 @@ -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) diff --git a/src/crystallite.f90 b/src/crystallite.f90 index 66c1df607..31a6bde2d 100644 --- a/src/crystallite.f90 +++ b/src/crystallite.f90 @@ -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 diff --git a/src/damage_local.f90 b/src/damage_local.f90 index e63db90b0..97eaf9a8c 100644 --- a/src/damage_local.f90 +++ b/src/damage_local.f90 @@ -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 diff --git a/src/damage_none.f90 b/src/damage_none.f90 index 2279bc06b..3f1144833 100644 --- a/src/damage_none.f90 +++ b/src/damage_none.f90 @@ -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 diff --git a/src/damage_nonlocal.f90 b/src/damage_nonlocal.f90 index 24a51cf54..c4426f185 100644 --- a/src/damage_nonlocal.f90 +++ b/src/damage_nonlocal.f90 @@ -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 diff --git a/src/grid/grid_mech_FEM.f90 b/src/grid/grid_mech_FEM.f90 index 741ce404a..cdf806b35 100644 --- a/src/grid/grid_mech_FEM.f90 +++ b/src/grid/grid_mech_FEM.f90 @@ -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) diff --git a/src/grid/grid_mech_spectral_basic.f90 b/src/grid/grid_mech_spectral_basic.f90 index d87a22fb7..ebaaf3b55 100644 --- a/src/grid/grid_mech_spectral_basic.f90 +++ b/src/grid/grid_mech_spectral_basic.f90 @@ -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 !-------------------------------------------------------------------------------------------------- diff --git a/src/grid/grid_mech_spectral_polarisation.f90 b/src/grid/grid_mech_spectral_polarisation.f90 index 9bbb40a53..9f2a17c97 100644 --- a/src/grid/grid_mech_spectral_polarisation.f90 +++ b/src/grid/grid_mech_spectral_polarisation.f90 @@ -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) diff --git a/src/grid/grid_thermal_spectral.f90 b/src/grid/grid_thermal_spectral.f90 index f5d1a33bc..259b45f33 100644 --- a/src/grid/grid_thermal_spectral.f90 +++ b/src/grid/grid_thermal_spectral.f90 @@ -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 diff --git a/src/grid/spectral_utilities.f90 b/src/grid/spectral_utilities.f90 index 989448dc3..c0c84233d 100644 --- a/src/grid/spectral_utilities.f90 +++ b/src/grid/spectral_utilities.f90 @@ -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 diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 347634212..57478e039 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -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)) diff --git a/src/homogenization_mech.f90 b/src/homogenization_mech.f90 new file mode 100644 index 000000000..b0641be07 --- /dev/null +++ b/src/homogenization_mech.f90 @@ -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 diff --git a/src/homogenization_mech_RGC.f90 b/src/homogenization_mech_RGC.f90 index 585752469..0a9d0ac92 100644 --- a/src/homogenization_mech_RGC.f90 +++ b/src/homogenization_mech_RGC.f90 @@ -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 diff --git a/src/homogenization_mech_isostrain.f90 b/src/homogenization_mech_isostrain.f90 index 751518e09..a56104647 100644 --- a/src/homogenization_mech_isostrain.f90 +++ b/src/homogenization_mech_isostrain.f90 @@ -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, & diff --git a/src/homogenization_mech_none.f90 b/src/homogenization_mech_none.f90 index 5b12247cd..d434d1ca0 100644 --- a/src/homogenization_mech_none.f90 +++ b/src/homogenization_mech_none.f90 @@ -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 @@ -28,7 +28,7 @@ module subroutine mech_none_init if(homogenization_Nconstituents(h) /= 1) & call IO_error(211,ext_msg='N_constituents (mech_none)') - + Nmaterialpoints = count(material_homogenizationAt == h) homogState(h)%sizeState = 0 allocate(homogState(h)%state0 (0,Nmaterialpoints)) diff --git a/src/kinematics_cleavage_opening.f90 b/src/kinematics_cleavage_opening.f90 index 66bce6a92..6fe8ed7f6 100644 --- a/src/kinematics_cleavage_opening.f90 +++ b/src/kinematics_cleavage_opening.f90 @@ -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 diff --git a/src/kinematics_slipplane_opening.f90 b/src/kinematics_slipplane_opening.f90 index aa0bdfbde..b7adb6807 100644 --- a/src/kinematics_slipplane_opening.f90 +++ b/src/kinematics_slipplane_opening.f90 @@ -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 diff --git a/src/kinematics_thermal_expansion.f90 b/src/kinematics_thermal_expansion.f90 index 2b8b04d85..36a882a48 100644 --- a/src/kinematics_thermal_expansion.f90 +++ b/src/kinematics_thermal_expansion.f90 @@ -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 * ( & diff --git a/src/material.f90 b/src/material.f90 index b05979298..bb5f484f6 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -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 diff --git a/src/mesh/FEM_utilities.f90 b/src/mesh/FEM_utilities.f90 index 118735e89..cb81f1f0c 100644 --- a/src/mesh/FEM_utilities.f90 +++ b/src/mesh/FEM_utilities.f90 @@ -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 diff --git a/src/mesh/mesh_mech_FEM.f90 b/src/mesh/mesh_mech_FEM.f90 index eb5f862c2..e19c35998 100644 --- a/src/mesh/mesh_mech_FEM.f90 +++ b/src/mesh/mesh_mech_FEM.f90 @@ -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 diff --git a/src/prec.f90 b/src/prec.f90 index 738775e3b..95b1116cd 100644 --- a/src/prec.f90 +++ b/src/prec.f90 @@ -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 diff --git a/src/source_damage_anisoBrittle.f90 b/src/source_damage_anisoBrittle.f90 index ca8d6ec2b..55d5546fc 100644 --- a/src/source_damage_anisoBrittle.f90 +++ b/src/source_damage_anisoBrittle.f90 @@ -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 diff --git a/src/source_damage_anisoDuctile.f90 b/src/source_damage_anisoDuctile.f90 index 601ec2531..912fe1387 100644 --- a/src/source_damage_anisoDuctile.f90 +++ b/src/source_damage_anisoDuctile.f90 @@ -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) & diff --git a/src/source_damage_isoDuctile.f90 b/src/source_damage_isoDuctile.f90 index 1bff20570..b66e220d9 100644 --- a/src/source_damage_isoDuctile.f90 +++ b/src/source_damage_isoDuctile.f90 @@ -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) = & diff --git a/src/thermal_adiabatic.f90 b/src/thermal_adiabatic.f90 index aa807924c..c67d004bf 100644 --- a/src/thermal_adiabatic.f90 +++ b/src/thermal_adiabatic.f90 @@ -72,12 +72,9 @@ subroutine thermal_adiabatic_init allocate(thermalState(h)%state0 (1,Nmaterialpoints), source=thermal_initialT(h)) 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 diff --git a/src/thermal_conduction.f90 b/src/thermal_conduction.f90 index daa7391a9..602bdab35 100644 --- a/src/thermal_conduction.f90 +++ b/src/thermal_conduction.f90 @@ -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 diff --git a/src/thermal_isothermal.f90 b/src/thermal_isothermal.f90 index 39c8efe91..adf2257de 100644 --- a/src/thermal_isothermal.f90 +++ b/src/thermal_isothermal.f90 @@ -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