diff --git a/src/commercialFEM_fileList.f90 b/src/commercialFEM_fileList.f90 index 585f242b5..1eac4f8fc 100644 --- a/src/commercialFEM_fileList.f90 +++ b/src/commercialFEM_fileList.f90 @@ -47,5 +47,8 @@ #include "homogenization_mechanical_isostrain.f90" #include "homogenization_mechanical_RGC.f90" #include "homogenization_thermal.f90" +#include "homogenization_thermal_pass.f90" +#include "homogenization_thermal_isotemperature.f90" #include "homogenization_damage.f90" +#include "homogenization_damage_pass.f90" #include "CPFEM.f90" diff --git a/src/homogenization.f90 b/src/homogenization.f90 index 7dfaf5b37..3e4e18f56 100644 --- a/src/homogenization.f90 +++ b/src/homogenization.f90 @@ -259,25 +259,25 @@ subroutine materialpoint_stressAndItsTangent(dt,FEsolving_execIP,FEsolving_execE NiterationMPstate, & ip, & !< integration point number el, & !< element number - myNgrains, co, ce, ho, me, ph + myNgrains, co, ce, ho, en, ph logical :: & converged logical, dimension(2) :: & doneAndHappy !$OMP PARALLEL - !$OMP DO PRIVATE(ce,me,ho,myNgrains,NiterationMPstate,converged,doneAndHappy) + !$OMP DO PRIVATE(ce,en,ho,myNgrains,NiterationMPstate,converged,doneAndHappy) do el = FEsolving_execElem(1),FEsolving_execElem(2) ho = material_homogenizationAt(el) myNgrains = homogenization_Nconstituents(ho) do ip = FEsolving_execIP(1),FEsolving_execIP(2) ce = (el-1)*discretization_nIPs + ip - me = material_homogenizationMemberAt2(ce) + en = material_homogenizationEntry(ce) call phase_restore(ce,.false.) ! wrong name (is more a forward function) - if(homogState(ho)%sizeState > 0) homogState(ho)%state(:,me) = homogState(ho)%state0(:,me) - if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%state(:,me) = damageState_h(ho)%state0(:,me) + if(homogState(ho)%sizeState > 0) homogState(ho)%state(:,en) = homogState(ho)%state0(:,en) + if(damageState_h(ho)%sizeState > 0) damageState_h(ho)%state(:,en) = damageState_h(ho)%state0(:,en) call damage_partition(ce) doneAndHappy = [.false.,.true.] @@ -505,12 +505,12 @@ function damage_nonlocal_getDiffusion(ce) ho, & co - ho = material_homogenizationAt2(ce) + ho = material_homogenizationID(ce) damage_nonlocal_getDiffusion = 0.0_pReal do co = 1, homogenization_Nconstituents(ho) damage_nonlocal_getDiffusion = damage_nonlocal_getDiffusion + & - crystallite_push33ToRef(co,ce,lattice_D(1:3,1:3,material_phaseAt2(co,ce))) + crystallite_push33ToRef(co,ce,lattice_D(1:3,1:3,material_phaseID(co,ce))) enddo damage_nonlocal_getDiffusion = & diff --git a/src/homogenization_damage.f90 b/src/homogenization_damage.f90 index 495cb0beb..94978ccc8 100644 --- a/src/homogenization_damage.f90 +++ b/src/homogenization_damage.f90 @@ -1,10 +1,17 @@ !-------------------------------------------------------------------------------------------------- !> @author Martin Diehl, KU Leuven !-------------------------------------------------------------------------------------------------- -submodule(homogenization) homogenization_damage +submodule(homogenization) damage use lattice + interface + + module subroutine pass_init + end subroutine pass_init + + end interface + type :: tDataContainer real(pReal), dimension(:), allocatable :: phi end type tDataContainer @@ -42,7 +49,7 @@ module subroutine damage_init() allocate(current(configHomogenizations%length)) do ho = 1, configHomogenizations%length - allocate(current(ho)%phi(count(material_homogenizationAt2==ho)), source=1.0_pReal) + allocate(current(ho)%phi(count(material_homogenizationID==ho)), source=1.0_pReal) configHomogenization => configHomogenizations%get(ho) associate(prm => param(ho)) if (configHomogenization%contains('damage')) then @@ -72,9 +79,9 @@ module subroutine damage_partition(ce) integer :: co - if(damageState_h(material_homogenizationAt2(ce))%sizeState < 1) return - phi = damagestate_h(material_homogenizationAt2(ce))%state(1,material_homogenizationMemberAt2(ce)) - do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) + if(damageState_h(material_homogenizationID(ce))%sizeState < 1) return + phi = damagestate_h(material_homogenizationID(ce))%state(1,material_homogenizationEntry(ce)) + do co = 1, homogenization_Nconstituents(material_homogenizationID(ce)) call phase_damage_set_phi(phi,co,ce) enddo @@ -94,11 +101,11 @@ module function damage_nonlocal_getMobility(ce) result(M) M = 0.0_pReal - do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) - M = M + lattice_M(material_phaseAt2(co,ce)) + do co = 1, homogenization_Nconstituents(material_homogenizationID(ce)) + M = M + lattice_M(material_phaseID(co,ce)) enddo - M = M/real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal) + M = M/real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal) end function damage_nonlocal_getMobility @@ -118,8 +125,8 @@ module subroutine damage_nonlocal_getSourceAndItsTangent(phiDot, dPhiDot_dPhi, p dPhiDot_dPhi = 0.0_pReal call phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, ce) - phiDot = phiDot/real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal) - dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal) + phiDot = phiDot/real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal) + dPhiDot_dPhi = dPhiDot_dPhi/real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal) end subroutine damage_nonlocal_getSourceAndItsTangent @@ -134,11 +141,11 @@ module subroutine damage_nonlocal_putNonLocalDamage(phi,ce) phi integer :: & ho, & - me + en - ho = material_homogenizationAt2(ce) - me = material_homogenizationMemberAt2(ce) - damagestate_h(ho)%state(1,me) = phi + ho = material_homogenizationID(ce) + en = material_homogenizationEntry(ce) + damagestate_h(ho)%state(1,en) = phi end subroutine damage_nonlocal_putNonLocalDamage @@ -165,4 +172,4 @@ module subroutine damage_nonlocal_results(ho,group) end subroutine damage_nonlocal_results -end submodule homogenization_damage +end submodule damage diff --git a/src/homogenization_damage_pass.f90 b/src/homogenization_damage_pass.f90 new file mode 100644 index 000000000..cbb7ec79f --- /dev/null +++ b/src/homogenization_damage_pass.f90 @@ -0,0 +1,14 @@ +!-------------------------------------------------------------------------------------------------- +!> @author Martin Diehl, KU Leuven +!> @brief Dummy homogenization scheme for 1 constituent per material point +!-------------------------------------------------------------------------------------------------- +submodule(homogenization:damage) damage_pass + +contains + +module subroutine pass_init + + +end subroutine pass_init + +end submodule damage_pass diff --git a/src/homogenization_mechanical.f90 b/src/homogenization_mechanical.f90 index f7074bc40..954aac403 100644 --- a/src/homogenization_mechanical.f90 +++ b/src/homogenization_mechanical.f90 @@ -7,51 +7,51 @@ submodule(homogenization) mechanical interface - module subroutine mechanical_pass_init - end subroutine mechanical_pass_init + module subroutine pass_init + end subroutine pass_init - module subroutine mechanical_isostrain_init - end subroutine mechanical_isostrain_init + module subroutine isostrain_init + end subroutine isostrain_init - module subroutine mechanical_RGC_init(num_homogMech) + module subroutine RGC_init(num_homogMech) class(tNode), pointer, intent(in) :: & num_homogMech !< pointer to mechanical homogenization numerics data - end subroutine mechanical_RGC_init + end subroutine RGC_init - module subroutine mechanical_isostrain_partitionDeformation(F,avgF) + module subroutine 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 mechanical_isostrain_partitionDeformation + end subroutine isostrain_partitionDeformation - module subroutine mechanical_RGC_partitionDeformation(F,avgF,ce) + module subroutine RGC_partitionDeformation(F,avgF,ce) 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) :: & ce - end subroutine mechanical_RGC_partitionDeformation + end subroutine RGC_partitionDeformation - module subroutine mechanical_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho) + module subroutine isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho) 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) :: ho - end subroutine mechanical_isostrain_averageStressAndItsTangent + end subroutine isostrain_averageStressAndItsTangent - module subroutine mechanical_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho) + module subroutine RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho) 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) :: ho - end subroutine mechanical_RGC_averageStressAndItsTangent + end subroutine RGC_averageStressAndItsTangent - module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy) + module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy) logical, dimension(2) :: doneAndHappy real(pReal), dimension(:,:,:), intent(in) :: & P,& !< partitioned stresses @@ -61,13 +61,13 @@ submodule(homogenization) mechanical real(pReal), intent(in) :: dt !< time increment integer, intent(in) :: & ce !< cell - end function mechanical_RGC_updateState + end function RGC_updateState - module subroutine mechanical_RGC_results(ho,group) + module subroutine RGC_results(ho,group) integer, intent(in) :: ho !< homogenization type character(len=*), intent(in) :: group !< group name in HDF5 file - end subroutine mechanical_RGC_results + end subroutine RGC_results end interface @@ -92,9 +92,9 @@ module subroutine mechanical_init(num_homog) 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 mechanical_pass_init - if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) call mechanical_isostrain_init - if (any(homogenization_type == HOMOGENIZATION_RGC_ID)) call mechanical_RGC_init(num_homogMech) + if (any(homogenization_type == HOMOGENIZATION_NONE_ID)) call pass_init + if (any(homogenization_type == HOMOGENIZATION_ISOSTRAIN_ID)) call isostrain_init + if (any(homogenization_type == HOMOGENIZATION_RGC_ID)) call RGC_init(num_homogMech) end subroutine mechanical_init @@ -110,23 +110,23 @@ module subroutine mechanical_partition(subF,ce) ce integer :: co - real(pReal), dimension (3,3,homogenization_Nconstituents(material_homogenizationAt2(ce))) :: Fs + real(pReal), dimension (3,3,homogenization_Nconstituents(material_homogenizationID(ce))) :: Fs - chosenHomogenization: select case(homogenization_type(material_homogenizationAt2(ce))) + chosenHomogenization: select case(homogenization_type(material_homogenizationID(ce))) case (HOMOGENIZATION_NONE_ID) chosenHomogenization Fs(1:3,1:3,1) = subF case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization - call mechanical_isostrain_partitionDeformation(Fs,subF) + call isostrain_partitionDeformation(Fs,subF) case (HOMOGENIZATION_RGC_ID) chosenHomogenization - call mechanical_RGC_partitionDeformation(Fs,subF,ce) + call RGC_partitionDeformation(Fs,subF,ce) end select chosenHomogenization - do co = 1,homogenization_Nconstituents(material_homogenizationAt2(ce)) + do co = 1,homogenization_Nconstituents(material_homogenizationID(ce)) call phase_mechanical_setF(Fs(1:3,1:3,co),co,ce) enddo @@ -143,37 +143,37 @@ module subroutine mechanical_homogenize(dt,ce) integer, intent(in) :: ce integer :: co - real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt2(ce))) - real(pReal) :: Ps(3,3,homogenization_Nconstituents(material_homogenizationAt2(ce))) + real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationID(ce))) + real(pReal) :: Ps(3,3,homogenization_Nconstituents(material_homogenizationID(ce))) - chosenHomogenization: select case(homogenization_type(material_homogenizationAt2(ce))) + chosenHomogenization: select case(homogenization_type(material_homogenizationID(ce))) case (HOMOGENIZATION_NONE_ID) chosenHomogenization homogenization_P(1:3,1:3,ce) = phase_mechanical_getP(1,ce) homogenization_dPdF(1:3,1:3,1:3,1:3,ce) = phase_mechanical_dPdF(dt,1,ce) case (HOMOGENIZATION_ISOSTRAIN_ID) chosenHomogenization - do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) + do co = 1, homogenization_Nconstituents(material_homogenizationID(ce)) dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(dt,co,ce) Ps(:,:,co) = phase_mechanical_getP(co,ce) enddo - call mechanical_isostrain_averageStressAndItsTangent(& + call isostrain_averageStressAndItsTangent(& homogenization_P(1:3,1:3,ce), & homogenization_dPdF(1:3,1:3,1:3,1:3,ce),& Ps,dPdFs, & - material_homogenizationAt2(ce)) + material_homogenizationID(ce)) case (HOMOGENIZATION_RGC_ID) chosenHomogenization - do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) + do co = 1, homogenization_Nconstituents(material_homogenizationID(ce)) dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(dt,co,ce) Ps(:,:,co) = phase_mechanical_getP(co,ce) enddo - call mechanical_RGC_averageStressAndItsTangent(& + call RGC_averageStressAndItsTangent(& homogenization_P(1:3,1:3,ce), & homogenization_dPdF(1:3,1:3,1:3,1:3,ce),& Ps,dPdFs, & - material_homogenizationAt2(ce)) + material_homogenizationID(ce)) end select chosenHomogenization @@ -195,18 +195,18 @@ module function mechanical_updateState(subdt,subF,ce) result(doneAndHappy) logical, dimension(2) :: doneAndHappy integer :: co - real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationAt2(ce))) - real(pReal) :: Fs(3,3,homogenization_Nconstituents(material_homogenizationAt2(ce))) - real(pReal) :: Ps(3,3,homogenization_Nconstituents(material_homogenizationAt2(ce))) + real(pReal) :: dPdFs(3,3,3,3,homogenization_Nconstituents(material_homogenizationID(ce))) + real(pReal) :: Fs(3,3,homogenization_Nconstituents(material_homogenizationID(ce))) + real(pReal) :: Ps(3,3,homogenization_Nconstituents(material_homogenizationID(ce))) - if (homogenization_type(material_homogenizationAt2(ce)) == HOMOGENIZATION_RGC_ID) then - do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) + if (homogenization_type(material_homogenizationID(ce)) == HOMOGENIZATION_RGC_ID) then + do co = 1, homogenization_Nconstituents(material_homogenizationID(ce)) dPdFs(:,:,:,:,co) = phase_mechanical_dPdF(subdt,co,ce) Fs(:,:,co) = phase_mechanical_getF(co,ce) Ps(:,:,co) = phase_mechanical_getP(co,ce) enddo - doneAndHappy = mechanical_RGC_updateState(Ps,Fs,subF,subdt,dPdFs,ce) + doneAndHappy = RGC_updateState(Ps,Fs,subF,subdt,dPdFs,ce) else doneAndHappy = .true. endif @@ -230,7 +230,7 @@ module subroutine mechanical_results(group_base,ho) select case(homogenization_type(ho)) case(HOMOGENIZATION_rgc_ID) - call mechanical_RGC_results(ho,group) + call RGC_results(ho,group) end select diff --git a/src/homogenization_mechanical_RGC.f90 b/src/homogenization_mechanical_RGC.f90 index 0de7a470b..772d1c4f5 100644 --- a/src/homogenization_mechanical_RGC.f90 +++ b/src/homogenization_mechanical_RGC.f90 @@ -71,7 +71,7 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief allocates all necessary fields, reads information from material configuration file !-------------------------------------------------------------------------------------------------- -module subroutine mechanical_RGC_init(num_homogMech) +module subroutine RGC_init(num_homogMech) class(tNode), pointer, intent(in) :: & num_homogMech !< pointer to mechanical homogenization numerics data @@ -152,7 +152,7 @@ module subroutine mechanical_RGC_init(num_homogMech) prm%N_constituents = homogMech%get_as1dInt('cluster_size',requiredSize=3) if (homogenization_Nconstituents(ho) /= product(prm%N_constituents)) & - call IO_error(211,ext_msg='N_constituents (mechanical_RGC)') + call IO_error(211,ext_msg='N_constituents (RGC)') prm%xi_alpha = homogMech%get_asFloat('xi_alpha') prm%c_alpha = homogMech%get_asFloat('c_alpha') @@ -187,13 +187,13 @@ module subroutine mechanical_RGC_init(num_homogMech) enddo -end subroutine mechanical_RGC_init +end subroutine RGC_init !-------------------------------------------------------------------------------------------------- !> @brief partitions the deformation gradient onto the constituents !-------------------------------------------------------------------------------------------------- -module subroutine mechanical_RGC_partitionDeformation(F,avgF,ce) +module subroutine RGC_partitionDeformation(F,avgF,ce) real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned F per grain @@ -204,12 +204,12 @@ module subroutine mechanical_RGC_partitionDeformation(F,avgF,ce) real(pReal), dimension(3) :: aVect,nVect integer, dimension(4) :: intFace integer, dimension(3) :: iGrain3 - integer :: iGrain,iFace,i,j,ho,me + integer :: iGrain,iFace,i,j,ho,en - associate(prm => param(material_homogenizationAt2(ce))) + associate(prm => param(material_homogenizationID(ce))) - ho = material_homogenizationAt2(ce) - me = material_homogenizationMemberAt2(ce) + ho = material_homogenizationID(ce) + en = material_homogenizationEntry(ce) !-------------------------------------------------------------------------------------------------- ! compute the deformation gradient of individual grains due to relaxations F = 0.0_pReal @@ -217,8 +217,8 @@ module subroutine mechanical_RGC_partitionDeformation(F,avgF,ce) iGrain3 = grain1to3(iGrain,prm%N_constituents) do iFace = 1,6 intFace = getInterface(iFace,iGrain3) ! identifying 6 interfaces of each grain - aVect = relaxationVector(intFace,ho,me) ! get the relaxation vectors for each interface from global relaxation vector array - nVect = interfaceNormal(intFace,ho,me) + aVect = relaxationVector(intFace,ho,en) ! get the relaxation vectors for each interface from global relaxation vector array + nVect = interfaceNormal(intFace,ho,en) forall (i=1:3,j=1:3) & F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! calculating deformation relaxations due to interface relaxation enddo @@ -227,14 +227,14 @@ module subroutine mechanical_RGC_partitionDeformation(F,avgF,ce) end associate -end subroutine mechanical_RGC_partitionDeformation +end subroutine RGC_partitionDeformation !-------------------------------------------------------------------------------------------------- !> @brief update the internal state of the homogenization scheme and tell whether "done" and ! "happy" with result !-------------------------------------------------------------------------------------------------- -module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy) +module function RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHappy) logical, dimension(2) :: doneAndHappy real(pReal), dimension(:,:,:), intent(in) :: & P,& !< partitioned stresses @@ -247,7 +247,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa integer, dimension(4) :: intFaceN,intFaceP,faceID integer, dimension(3) :: nGDim,iGr3N,iGr3P - integer :: ho,iNum,i,j,nIntFaceTot,iGrN,iGrP,iMun,iFace,k,l,ipert,iGrain,nGrain, me + integer :: ho,iNum,i,j,nIntFaceTot,iGrN,iGrP,iMun,iFace,k,l,ipert,nGrain, en real(pReal), dimension(3,3,size(P,3)) :: R,pF,pR,D,pD real(pReal), dimension(3,size(P,3)) :: NN,devNull real(pReal), dimension(3) :: normP,normN,mornP,mornN @@ -261,9 +261,9 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa return endif zeroTimeStep - ho = material_homogenizationAt2(ce) + ho = material_homogenizationID(ce) + en = material_homogenizationEntry(ce) - me = material_homogenizationMemberAt2(ce) associate(stt => state(ho), st0 => state0(ho), dst => dependentState(ho), prm => param(ho)) !-------------------------------------------------------------------------------------------------- @@ -278,16 +278,16 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa ! allocate the size of the global relaxation arrays/jacobian matrices depending on the size of the cluster allocate(resid(3*nIntFaceTot), source=0.0_pReal) allocate(tract(nIntFaceTot,3), source=0.0_pReal) - relax = stt%relaxationVector(:,me) - drelax = stt%relaxationVector(:,me) - st0%relaxationVector(:,me) + relax = stt%relaxationVector(:,en) + drelax = stt%relaxationVector(:,en) - st0%relaxationVector(:,en) !-------------------------------------------------------------------------------------------------- ! computing interface mismatch and stress penalty tensor for all interfaces of all grains - call stressPenalty(R,NN,avgF,F,ho,me) + call stressPenalty(R,NN,avgF,F,ho,en) !-------------------------------------------------------------------------------------------------- ! calculating volume discrepancy and stress penalty related to overall volume discrepancy - call volumePenalty(D,dst%volumeDiscrepancy(me),avgF,F,nGrain) + call volumePenalty(D,dst%volumeDiscrepancy(en),avgF,F,nGrain) !------------------------------------------------------------------------------------------------ ! computing the residual stress from the balance of traction at all (interior) interfaces @@ -299,7 +299,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa iGr3N = faceID(2:4) ! identifying the grain ID in local coordinate system (3-dimensional index) iGrN = grain3to1(iGr3N,param(ho)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index) intFaceN = getInterface(2*faceID(1),iGr3N) - normN = interfaceNormal(intFaceN,ho,me) + normN = interfaceNormal(intFaceN,ho,en) !-------------------------------------------------------------------------------------------------- ! identify the right/up/front grain (+|P) @@ -307,7 +307,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identifying the grain ID in local coordinate system (3-dimensional index) iGrP = grain3to1(iGr3P,param(ho)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index) intFaceP = getInterface(2*faceID(1)-1,iGr3P) - normP = interfaceNormal(intFaceP,ho,me) + normP = interfaceNormal(intFaceP,ho,en) !-------------------------------------------------------------------------------------------------- ! compute the residual of traction at the interface (in local system, 4-dimensional index) @@ -335,9 +335,9 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa if (residMax < num%rtol*stresMax .or. residMax < num%atol) then doneAndHappy = .true. - dst%mismatch(1:3,me) = sum(NN,2)/real(nGrain,pReal) - dst%relaxationRate_avg(me) = sum(abs(drelax))/dt/real(3*nIntFaceTot,pReal) - dst%relaxationRate_max(me) = maxval(abs(drelax))/dt + dst%mismatch(1:3,en) = sum(NN,2)/real(nGrain,pReal) + dst%relaxationRate_avg(en) = sum(abs(drelax))/dt/real(3*nIntFaceTot,pReal) + dst%relaxationRate_max(en) = maxval(abs(drelax))/dt return @@ -363,10 +363,10 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa iGr3N = faceID(2:4) ! identifying the grain ID in local coordinate sytem iGrN = grain3to1(iGr3N,param(ho)%N_constituents) ! translate into global grain ID intFaceN = getInterface(2*faceID(1),iGr3N) ! identifying the connecting interface in local coordinate system - normN = interfaceNormal(intFaceN,ho,me) + normN = interfaceNormal(intFaceN,ho,en) do iFace = 1,6 intFaceN = getInterface(iFace,iGr3N) ! identifying all interfaces that influence relaxation of the above interface - mornN = interfaceNormal(intFaceN,ho,me) + mornN = interfaceNormal(intFaceN,ho,en) iMun = interface4to1(intFaceN,param(ho)%N_constituents) ! translate the interfaces ID into local 4-dimensional index if (iMun > 0) then ! get the corresponding tangent do i=1,3; do j=1,3; do k=1,3; do l=1,3 @@ -384,10 +384,10 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identifying the grain ID in local coordinate sytem iGrP = grain3to1(iGr3P,param(ho)%N_constituents) ! translate into global grain ID intFaceP = getInterface(2*faceID(1)-1,iGr3P) ! identifying the connecting interface in local coordinate system - normP = interfaceNormal(intFaceP,ho,me) + normP = interfaceNormal(intFaceP,ho,en) do iFace = 1,6 intFaceP = getInterface(iFace,iGr3P) ! identifying all interfaces that influence relaxation of the above interface - mornP = interfaceNormal(intFaceP,ho,me) + mornP = interfaceNormal(intFaceP,ho,en) iMun = interface4to1(intFaceP,param(ho)%N_constituents) ! translate the interfaces ID into local 4-dimensional index if (iMun > 0) then ! get the corresponding tangent do i=1,3; do j=1,3; do k=1,3; do l=1,3 @@ -408,9 +408,9 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa do ipert = 1,3*nIntFaceTot p_relax = relax p_relax(ipert) = relax(ipert) + num%pPert ! perturb the relaxation vector - stt%relaxationVector(:,me) = p_relax - call grainDeformation(pF,avgF,ho,me) ! rain deformation from perturbed state - call stressPenalty(pR,DevNull, avgF,pF,ho,me) ! stress penalty due to interface mismatch from perturbed state + stt%relaxationVector(:,en) = p_relax + call grainDeformation(pF,avgF,ho,en) ! rain deformation from perturbed state + call stressPenalty(pR,DevNull, avgF,pF,ho,en) ! stress penalty due to interface mismatch from perturbed state call volumePenalty(pD,devNull(1,1), avgF,pF,nGrain) ! stress penalty due to volume discrepancy from perturbed state !-------------------------------------------------------------------------------------------------- @@ -424,7 +424,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa iGr3N = faceID(2:4) ! identify the grain ID in local coordinate system (3-dimensional index) iGrN = grain3to1(iGr3N,param(ho)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index) intFaceN = getInterface(2*faceID(1),iGr3N) ! identify the interface ID of the grain - normN = interfaceNormal(intFaceN,ho,me) + normN = interfaceNormal(intFaceN,ho,en) !-------------------------------------------------------------------------------------------------- ! identify the right/up/front grain (+|P) @@ -432,7 +432,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa iGr3P(faceID(1)) = iGr3N(faceID(1))+1 ! identify the grain ID in local coordinate system (3-dimensional index) iGrP = grain3to1(iGr3P,param(ho)%N_constituents) ! translate the local grain ID into global coordinate system (1-dimensional index) intFaceP = getInterface(2*faceID(1)-1,iGr3P) ! identify the interface ID of the grain - normP = interfaceNormal(intFaceP,ho,me) + normP = interfaceNormal(intFaceP,ho,en) !-------------------------------------------------------------------------------------------------- ! compute the residual stress (contribution of mismatch and volume penalties) from perturbed state @@ -472,7 +472,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa do i = 1,3*nIntFaceTot;do j = 1,3*nIntFaceTot drelax(i) = drelax(i) - jnverse(i,j)*resid(j) ! Calculate the correction for the state variable enddo; enddo - stt%relaxationVector(:,me) = relax + drelax ! Updateing the state variable for the next iteration + stt%relaxationVector(:,en) = relax + drelax ! Updateing the state variable for the next iteration if (any(abs(drelax) > num%maxdRelax)) then ! Forcing cutback when the incremental change of relaxation vector becomes too large doneAndHappy = [.true.,.false.] !$OMP CRITICAL (write2out) @@ -488,14 +488,14 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa !------------------------------------------------------------------------------------------------ !> @brief calculate stress-like penalty due to deformation mismatch !------------------------------------------------------------------------------------------------ - subroutine stressPenalty(rPen,nMis,avgF,fDef,ho,me) + subroutine stressPenalty(rPen,nMis,avgF,fDef,ho,en) real(pReal), dimension (:,:,:), intent(out) :: rPen !< stress-like penalty real(pReal), dimension (:,:), intent(out) :: nMis !< total amount of mismatch real(pReal), dimension (:,:,:), intent(in) :: fDef !< deformation gradients real(pReal), dimension (3,3), intent(in) :: avgF !< initial effective stretch tensor - integer, intent(in) :: ho, me + integer, intent(in) :: ho, en integer, dimension (4) :: intFace integer, dimension (3) :: iGrain3,iGNghb3,nGDim @@ -515,7 +515,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa ! get the correction factor the modulus of penalty stress representing the evolution of area of ! the interfaces due to deformations - surfCorr = surfaceCorrection(avgF,ho,me) + surfCorr = surfaceCorrection(avgF,ho,en) associate(prm => param(ho)) @@ -527,7 +527,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa interfaceLoop: do iFace = 1,6 intFace = getInterface(iFace,iGrain3) ! get the 4-dimensional index of the interface in local numbering system of the grain - nVect = interfaceNormal(intFace,ho,me) + nVect = interfaceNormal(intFace,ho,en) iGNghb3 = iGrain3 ! identify the neighboring grain across the interface iGNghb3(abs(intFace(1))) = iGNghb3(abs(intFace(1))) & + int(real(intFace(1),pReal)/real(abs(intFace(1)),pReal)) @@ -611,14 +611,14 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa !> @brief compute the correction factor accouted for surface evolution (area change) due to ! deformation !-------------------------------------------------------------------------------------------------- - function surfaceCorrection(avgF,ho,me) + function surfaceCorrection(avgF,ho,en) real(pReal), dimension(3) :: surfaceCorrection real(pReal), dimension(3,3), intent(in) :: avgF !< average F integer, intent(in) :: & ho, & - me + en real(pReal), dimension(3,3) :: invC real(pReal), dimension(3) :: nVect real(pReal) :: detF @@ -629,7 +629,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa surfaceCorrection = 0.0_pReal do iBase = 1,3 - nVect = interfaceNormal([iBase,1,1,1],ho,me) + nVect = interfaceNormal([iBase,1,1,1],ho,en) do i = 1,3; do j = 1,3 surfaceCorrection(iBase) = surfaceCorrection(iBase) + invC(i,j)*nVect(i)*nVect(j) ! compute the component of (the inverse of) the stretch in the direction of the normal enddo; enddo @@ -651,7 +651,7 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa real(pReal), dimension(6,6) :: C - C = phase_homogenizedC(material_phaseAt2(grainID,ce),material_phaseMemberAt2(grainID,ce)) + C = phase_homogenizedC(material_phaseID(grainID,ce),material_phaseEntry(grainID,ce)) equivalentMu = lattice_equivalent_mu(C,'voigt') end function equivalentMu @@ -661,14 +661,14 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa !> @brief calculating the grain deformation gradient (the same with ! homogenization_RGC_partitionDeformation, but used only for perturbation scheme) !------------------------------------------------------------------------------------------------- - subroutine grainDeformation(F, avgF, ho, me) + subroutine grainDeformation(F, avgF, ho, en) real(pReal), dimension(:,:,:), intent(out) :: F !< partitioned F per grain real(pReal), dimension(:,:), intent(in) :: avgF !< averaged F integer, intent(in) :: & ho, & - me + en real(pReal), dimension(3) :: aVect,nVect integer, dimension(4) :: intFace @@ -685,8 +685,8 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa iGrain3 = grain1to3(iGrain,prm%N_constituents) do iFace = 1,6 intFace = getInterface(iFace,iGrain3) - aVect = relaxationVector(intFace,ho,me) - nVect = interfaceNormal(intFace,ho,me) + aVect = relaxationVector(intFace,ho,en) + nVect = interfaceNormal(intFace,ho,en) forall (i=1:3,j=1:3) & F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! effective relaxations enddo @@ -697,13 +697,13 @@ module function mechanical_RGC_updateState(P,F,avgF,dt,dPdF,ce) result(doneAndHa end subroutine grainDeformation -end function mechanical_RGC_updateState +end function RGC_updateState !-------------------------------------------------------------------------------------------------- !> @brief derive average stress and stiffness from constituent quantities !-------------------------------------------------------------------------------------------------- -module subroutine mechanical_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho) +module subroutine RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho) 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 @@ -715,13 +715,13 @@ module subroutine mechanical_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dP avgP = sum(P,3) /real(product(param(ho)%N_constituents),pReal) dAvgPdAvgF = sum(dPdF,5)/real(product(param(ho)%N_constituents),pReal) -end subroutine mechanical_RGC_averageStressAndItsTangent +end subroutine RGC_averageStressAndItsTangent !-------------------------------------------------------------------------------------------------- !> @brief writes results to HDF5 output file !-------------------------------------------------------------------------------------------------- -module subroutine mechanical_RGC_results(ho,group) +module subroutine RGC_results(ho,group) integer, intent(in) :: ho character(len=*), intent(in) :: group @@ -747,17 +747,17 @@ module subroutine mechanical_RGC_results(ho,group) enddo outputsLoop end associate -end subroutine mechanical_RGC_results +end subroutine RGC_results !-------------------------------------------------------------------------------------------------- !> @brief collect relaxation vectors of an interface !-------------------------------------------------------------------------------------------------- -pure function relaxationVector(intFace,ho,me) +pure function relaxationVector(intFace,ho,en) real(pReal), dimension (3) :: relaxationVector - integer, intent(in) :: ho,me + integer, intent(in) :: ho,en integer, dimension(4), intent(in) :: intFace !< set of interface ID in 4D array (normal and position) integer :: iNum @@ -770,7 +770,7 @@ pure function relaxationVector(intFace,ho,me) iNum = interface4to1(intFace,prm%N_constituents) ! identify the position of the interface in global state array if (iNum > 0) then - relaxationVector = stt%relaxationVector((3*iNum-2):(3*iNum),me) + relaxationVector = stt%relaxationVector((3*iNum-2):(3*iNum),en) else relaxationVector = 0.0_pReal endif @@ -783,14 +783,14 @@ end function relaxationVector !-------------------------------------------------------------------------------------------------- !> @brief identify the normal of an interface !-------------------------------------------------------------------------------------------------- -pure function interfaceNormal(intFace,ho,me) +pure function interfaceNormal(intFace,ho,en) real(pReal), dimension(3) :: interfaceNormal integer, dimension(4), intent(in) :: intFace !< interface ID in 4D array (normal and position) integer, intent(in) :: & ho, & - me + en integer :: nPos associate (dst => dependentState(ho)) @@ -801,7 +801,7 @@ pure function interfaceNormal(intFace,ho,me) nPos = abs(intFace(1)) ! identify the position of the interface in global state array interfaceNormal(nPos) = real(intFace(1)/abs(intFace(1)),pReal) ! get the normal vector w.r.t. cluster axis - interfaceNormal = matmul(dst%orientation(1:3,1:3,me),interfaceNormal) ! map the normal vector into sample coordinate system (basis) + interfaceNormal = matmul(dst%orientation(1:3,1:3,en),interfaceNormal) ! map the normal vector into sample coordinate system (basis) end associate diff --git a/src/homogenization_mechanical_isostrain.f90 b/src/homogenization_mechanical_isostrain.f90 index e7d4cff4c..9a3704575 100644 --- a/src/homogenization_mechanical_isostrain.f90 +++ b/src/homogenization_mechanical_isostrain.f90 @@ -26,7 +26,7 @@ contains !-------------------------------------------------------------------------------------------------- !> @brief allocates all neccessary fields, reads information from material configuration file !-------------------------------------------------------------------------------------------------- -module subroutine mechanical_isostrain_init +module subroutine isostrain_init integer :: & h, & @@ -56,7 +56,7 @@ module subroutine mechanical_isostrain_init case ('avg') prm%mapping = average_ID case default - call IO_error(211,ext_msg='sum'//' (mechanical_isostrain)') + call IO_error(211,ext_msg='sum'//' (isostrain)') end select Nmaterialpoints = count(material_homogenizationAt == h) @@ -68,13 +68,13 @@ module subroutine mechanical_isostrain_init enddo -end subroutine mechanical_isostrain_init +end subroutine isostrain_init !-------------------------------------------------------------------------------------------------- !> @brief partitions the deformation gradient onto the constituents !-------------------------------------------------------------------------------------------------- -module subroutine mechanical_isostrain_partitionDeformation(F,avgF) +module subroutine isostrain_partitionDeformation(F,avgF) real(pReal), dimension (:,:,:), intent(out) :: F !< partitioned deformation gradient @@ -82,13 +82,13 @@ module subroutine mechanical_isostrain_partitionDeformation(F,avgF) F = spread(avgF,3,size(F,3)) -end subroutine mechanical_isostrain_partitionDeformation +end subroutine isostrain_partitionDeformation !-------------------------------------------------------------------------------------------------- !> @brief derive average stress and stiffness from constituent quantities !-------------------------------------------------------------------------------------------------- -module subroutine mechanical_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho) +module subroutine isostrain_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,ho) 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 @@ -110,6 +110,6 @@ module subroutine mechanical_isostrain_averageStressAndItsTangent(avgP,dAvgPdAvg end associate -end subroutine mechanical_isostrain_averageStressAndItsTangent +end subroutine isostrain_averageStressAndItsTangent end submodule isostrain diff --git a/src/homogenization_mechanical_pass.f90 b/src/homogenization_mechanical_pass.f90 index 6217e6836..23fe74f44 100644 --- a/src/homogenization_mechanical_pass.f90 +++ b/src/homogenization_mechanical_pass.f90 @@ -4,14 +4,14 @@ !> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH !> @brief dummy homogenization homogenization scheme for 1 constituent per material point !-------------------------------------------------------------------------------------------------- -submodule(homogenization:mechanical) none +submodule(homogenization:mechanical) mechanical_pass contains !-------------------------------------------------------------------------------------------------- !> @brief allocates all necessary fields, reads information from material configuration file !-------------------------------------------------------------------------------------------------- -module subroutine mechanical_pass_init +module subroutine pass_init integer :: & Ninstances, & @@ -27,7 +27,7 @@ module subroutine mechanical_pass_init if(homogenization_type(h) /= HOMOGENIZATION_NONE_ID) cycle if(homogenization_Nconstituents(h) /= 1) & - call IO_error(211,ext_msg='N_constituents (mechanical_pass)') + call IO_error(211,ext_msg='N_constituents (pass)') Nmaterialpoints = count(material_homogenizationAt == h) homogState(h)%sizeState = 0 @@ -36,6 +36,6 @@ module subroutine mechanical_pass_init enddo -end subroutine mechanical_pass_init +end subroutine pass_init -end submodule none +end submodule mechanical_pass diff --git a/src/homogenization_thermal.f90 b/src/homogenization_thermal.f90 index 7602e50d2..7bb81d24a 100644 --- a/src/homogenization_thermal.f90 +++ b/src/homogenization_thermal.f90 @@ -5,6 +5,16 @@ submodule(homogenization) thermal use lattice + interface + + module subroutine pass_init + end subroutine pass_init + + module subroutine isotemperature_init + end subroutine isotemperature_init + + end interface + type :: tDataContainer real(pReal), dimension(:), allocatable :: T, dot_T end type tDataContainer @@ -44,8 +54,8 @@ module subroutine thermal_init() allocate(current(configHomogenizations%length)) do ho = 1, configHomogenizations%length - allocate(current(ho)%T(count(material_homogenizationAt2==ho)), source=300.0_pReal) - allocate(current(ho)%dot_T(count(material_homogenizationAt2==ho)), source=0.0_pReal) + allocate(current(ho)%T(count(material_homogenizationID==ho)), source=300.0_pReal) + allocate(current(ho)%dot_T(count(material_homogenizationID==ho)), source=0.0_pReal) configHomogenization => configHomogenizations%get(ho) associate(prm => param(ho)) if (configHomogenization%contains('thermal')) then @@ -75,9 +85,9 @@ module subroutine thermal_partition(ce) integer :: co - T = current(material_homogenizationAt2(ce))%T(material_homogenizationMemberAt2(ce)) - dot_T = current(material_homogenizationAt2(ce))%dot_T(material_homogenizationMemberAt2(ce)) - do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) + T = current(material_homogenizationID(ce))%T(material_homogenizationEntry(ce)) + dot_T = current(material_homogenizationID(ce))%dot_T(material_homogenizationEntry(ce)) + do co = 1, homogenization_Nconstituents(material_homogenizationID(ce)) call phase_thermal_setField(T,dot_T,co,ce) enddo @@ -109,11 +119,11 @@ module function thermal_conduction_getConductivity(ce) result(K) K = 0.0_pReal - do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) - K = K + crystallite_push33ToRef(co,ce,lattice_K(:,:,material_phaseAt2(co,ce))) + do co = 1, homogenization_Nconstituents(material_homogenizationID(ce)) + K = K + crystallite_push33ToRef(co,ce,lattice_K(:,:,material_phaseID(co,ce))) enddo - K = K / real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal) + K = K / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal) end function thermal_conduction_getConductivity @@ -131,11 +141,11 @@ module function thermal_conduction_getSpecificHeat(ce) result(c_P) c_P = 0.0_pReal - do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) - c_P = c_P + lattice_c_p(material_phaseAt2(co,ce)) + do co = 1, homogenization_Nconstituents(material_homogenizationID(ce)) + c_P = c_P + lattice_c_p(material_phaseID(co,ce)) enddo - c_P = c_P / real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal) + c_P = c_P / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal) end function thermal_conduction_getSpecificHeat @@ -153,11 +163,11 @@ module function thermal_conduction_getMassDensity(ce) result(rho) rho = 0.0_pReal - do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) - rho = rho + lattice_rho(material_phaseAt2(co,ce)) + do co = 1, homogenization_Nconstituents(material_homogenizationID(ce)) + rho = rho + lattice_rho(material_phaseID(co,ce)) enddo - rho = rho / real(homogenization_Nconstituents(material_homogenizationAt2(ce)),pReal) + rho = rho / real(homogenization_Nconstituents(material_homogenizationID(ce)),pReal) end function thermal_conduction_getMassDensity @@ -172,8 +182,8 @@ module subroutine homogenization_thermal_setField(T,dot_T, ce) real(pReal), intent(in) :: T, dot_T - current(material_homogenizationAt2(ce))%T(material_homogenizationMemberAt2(ce)) = T - current(material_homogenizationAt2(ce))%dot_T(material_homogenizationMemberAt2(ce)) = dot_T + current(material_homogenizationID(ce))%T(material_homogenizationEntry(ce)) = T + current(material_homogenizationID(ce))%dot_T(material_homogenizationEntry(ce)) = dot_T end subroutine homogenization_thermal_setField @@ -207,7 +217,7 @@ module function homogenization_thermal_T(ce) result(T) integer, intent(in) :: ce real(pReal) :: T - T = current(material_homogenizationAt2(ce))%T(material_homogenizationMemberAt2(ce)) + T = current(material_homogenizationID(ce))%T(material_homogenizationEntry(ce)) end function homogenization_thermal_T diff --git a/src/homogenization_thermal_isotemperature.f90 b/src/homogenization_thermal_isotemperature.f90 new file mode 100644 index 000000000..7358ecf08 --- /dev/null +++ b/src/homogenization_thermal_isotemperature.f90 @@ -0,0 +1,14 @@ +!-------------------------------------------------------------------------------------------------- +!> @author Martin Diehl, KU Leuven +!> @brief Dummy homogenization scheme for 1 constituent per material point +!-------------------------------------------------------------------------------------------------- +submodule(homogenization:thermal) isotemperature + +contains + +module subroutine isotemperature_init + + +end subroutine isotemperature_init + +end submodule isotemperature diff --git a/src/homogenization_thermal_pass.f90 b/src/homogenization_thermal_pass.f90 new file mode 100644 index 000000000..940a15024 --- /dev/null +++ b/src/homogenization_thermal_pass.f90 @@ -0,0 +1,14 @@ +!-------------------------------------------------------------------------------------------------- +!> @author Martin Diehl, KU Leuven +!> @brief Dummy homogenization scheme for 1 constituent per material point +!-------------------------------------------------------------------------------------------------- +submodule(homogenization:thermal) thermal_pass + +contains + +module subroutine pass_init + + +end subroutine pass_init + +end submodule thermal_pass diff --git a/src/material.f90 b/src/material.f90 index aecfaa1dd..6575872ed 100644 --- a/src/material.f90 +++ b/src/material.f90 @@ -28,14 +28,14 @@ module material integer, dimension(:), allocatable, public, protected :: & ! (elem) material_homogenizationAt, & !< homogenization ID of each element - material_homogenizationAt2, & !< per cell - material_homogenizationMemberAt2 !< cell + material_homogenizationID, & !< per cell + material_homogenizationEntry !< cell integer, dimension(:,:), allocatable :: & ! (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 - material_phaseAt2, & !< per constituent,cell - material_phaseMemberAt2 !< per constituent, cell + material_phaseID, & !< per constituent,cell + material_phaseEntry !< per constituent, cell integer, dimension(:,:,:), allocatable, public, protected :: & ! (constituent,IP,elem) material_phaseMemberAt !< position of the element within its phase instance @@ -117,10 +117,10 @@ subroutine material_parseMaterial allocate(material_phaseMemberAt(homogenization_maxNconstituents,discretization_nIPs,discretization_Nelems),source=0) - allocate(material_homogenizationAt2(discretization_nIPs*discretization_Nelems),source=0) - allocate(material_homogenizationMemberAt2(discretization_nIPs*discretization_Nelems),source=0) - allocate(material_phaseAt2(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0) - allocate(material_phaseMemberAt2(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0) + allocate(material_homogenizationID(discretization_nIPs*discretization_Nelems),source=0) + allocate(material_homogenizationEntry(discretization_nIPs*discretization_Nelems),source=0) + allocate(material_phaseID(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0) + allocate(material_phaseEntry(homogenization_maxNconstituents,discretization_nIPs*discretization_Nelems),source=0) do el = 1, discretization_Nelems material => materials%get(discretization_materialAt(el)) @@ -131,8 +131,8 @@ subroutine material_parseMaterial ce = (el-1)*discretization_nIPs + ip counterHomogenization(material_homogenizationAt(el)) = counterHomogenization(material_homogenizationAt(el)) + 1 material_homogenizationMemberAt(ip,el) = counterHomogenization(material_homogenizationAt(el)) - material_homogenizationAt2(ce) = material_homogenizationAt(el) - material_homogenizationMemberAt2(ce) = material_homogenizationMemberAt(ip,el) + material_homogenizationID(ce) = material_homogenizationAt(el) + material_homogenizationEntry(ce) = material_homogenizationMemberAt(ip,el) enddo frac = 0.0_pReal @@ -146,8 +146,8 @@ subroutine material_parseMaterial counterPhase(material_phaseAt(co,el)) = counterPhase(material_phaseAt(co,el)) + 1 material_phaseMemberAt(co,ip,el) = counterPhase(material_phaseAt(co,el)) - material_phaseAt2(co,ce) = material_phaseAt(co,el) - material_phaseMemberAt2(co,ce) = material_phaseMemberAt(co,ip,el) + material_phaseID(co,ce) = material_phaseAt(co,el) + material_phaseEntry(co,ce) = material_phaseMemberAt(co,ip,el) enddo enddo diff --git a/src/phase.f90 b/src/phase.f90 index 9470ab8af..d10df84ba 100644 --- a/src/phase.f90 +++ b/src/phase.f90 @@ -100,14 +100,12 @@ module phase integer, intent(in) :: ph end subroutine damage_results - module subroutine mechanical_windForward(ph,me) - integer, intent(in) :: ph, me - end subroutine mechanical_windForward - - module subroutine mechanical_forward() end subroutine mechanical_forward + module subroutine damage_forward() + end subroutine damage_forward + module subroutine thermal_forward() end subroutine thermal_forward @@ -261,7 +259,7 @@ module phase end subroutine plastic_dependentState - module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) + module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) integer, intent(in) :: ph, me real(pReal), intent(in), dimension(3,3) :: & S @@ -269,9 +267,9 @@ module phase Ld !< damage velocity gradient real(pReal), intent(out), dimension(3,3,3,3) :: & dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) - end subroutine kinematics_cleavage_opening_LiAndItsTangent + end subroutine damage_anisobrittle_LiAndItsTangent - module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) + module subroutine damage_isoductile_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) integer, intent(in) :: ph, me real(pReal), intent(in), dimension(3,3) :: & S @@ -279,7 +277,7 @@ module phase Ld !< damage velocity gradient real(pReal), intent(out), dimension(3,3,3,3) :: & dLd_dTstar !< derivative of Ld with respect to Tstar (4th-order tensor) - end subroutine kinematics_slipplane_opening_LiAndItsTangent + end subroutine damage_isoductile_LiAndItsTangent end interface @@ -325,13 +323,12 @@ module phase phase_damage_get_phi, & phase_mechanical_getP, & phase_mechanical_setF, & - phase_mechanical_getF, & - phase_windForward + phase_mechanical_getF contains !-------------------------------------------------------------------------------------------------- -!> @brief Initialze constitutive models for individual physics +!> @brief Initialize constitutive models for individual physics !-------------------------------------------------------------------------------------------------- subroutine phase_init @@ -382,12 +379,12 @@ end subroutine phase_init !> @brief Allocate the components of the state structure for a given phase !-------------------------------------------------------------------------------------------------- subroutine phase_allocateState(state, & - Nconstituents,sizeState,sizeDotState,sizeDeltaState) + NEntries,sizeState,sizeDotState,sizeDeltaState) class(tState), intent(out) :: & state integer, intent(in) :: & - Nconstituents, & + NEntries, & sizeState, & sizeDotState, & sizeDeltaState @@ -398,13 +395,13 @@ subroutine phase_allocateState(state, & state%sizeDeltaState = sizeDeltaState state%offsetDeltaState = sizeState-sizeDeltaState ! deltaState occupies latter part of state by definition - allocate(state%atol (sizeState), source=0.0_pReal) - allocate(state%state0 (sizeState,Nconstituents), source=0.0_pReal) - allocate(state%state (sizeState,Nconstituents), source=0.0_pReal) + allocate(state%atol (sizeState), source=0.0_pReal) + allocate(state%state0 (sizeState,NEntries), source=0.0_pReal) + allocate(state%state (sizeState,NEntries), source=0.0_pReal) - allocate(state%dotState (sizeDotState,Nconstituents), source=0.0_pReal) + allocate(state%dotState (sizeDotState,NEntries), source=0.0_pReal) - allocate(state%deltaState (sizeDeltaState,Nconstituents), source=0.0_pReal) + allocate(state%deltaState (sizeDeltaState,NEntries), source=0.0_pReal) end subroutine phase_allocateState @@ -422,10 +419,10 @@ subroutine phase_restore(ce,includeL) co - do co = 1,homogenization_Nconstituents(material_homogenizationAt2(ce)) - if (damageState(material_phaseAt2(co,ce))%sizeState > 0) & - damageState(material_phaseAt2(co,ce))%state( :,material_phasememberAt2(co,ce)) = & - damageState(material_phaseAt2(co,ce))%state0(:,material_phasememberAt2(co,ce)) + do co = 1,homogenization_Nconstituents(material_homogenizationID(ce)) + if (damageState(material_phaseID(co,ce))%sizeState > 0) & + damageState(material_phaseID(co,ce))%state( :,material_phaseEntry(co,ce)) = & + damageState(material_phaseID(co,ce))%state0(:,material_phaseEntry(co,ce)) enddo call mechanical_restore(ce,includeL) @@ -435,21 +432,13 @@ end subroutine phase_restore !-------------------------------------------------------------------------------------------------- !> @brief Forward data after successful increment. -! ToDo: Any guessing for the current states possible? !-------------------------------------------------------------------------------------------------- subroutine phase_forward() - integer :: ph - - call mechanical_forward() + call damage_forward() call thermal_forward() - do ph = 1, size(damageState) - if (damageState(ph)%sizeState > 0) & - damageState(ph)%state0 = damageState(ph)%state - enddo - end subroutine phase_forward @@ -484,7 +473,6 @@ subroutine crystallite_init() integer :: & ph, & - me, & co, & !< counter in integration point component loop ip, & !< counter in integration point loop el, & !< counter in element loop @@ -550,12 +538,10 @@ subroutine crystallite_init() flush(IO_STDOUT) - !$OMP PARALLEL DO PRIVATE(ph,me) + !$OMP PARALLEL DO do el = 1, size(material_phaseMemberAt,3) do ip = 1, size(material_phaseMemberAt,2) do co = 1,homogenization_Nconstituents(material_homogenizationAt(el)) - ph = material_phaseAt(co,el) - me = material_phaseMemberAt(co,ip,el) call crystallite_orientations(co,ip,el) call plastic_dependentState(co,ip,el) ! update dependent state variables to be consistent with basic states enddo @@ -567,34 +553,6 @@ subroutine crystallite_init() end subroutine crystallite_init -!-------------------------------------------------------------------------------------------------- -!> @brief Wind homog inc forward. -!-------------------------------------------------------------------------------------------------- -subroutine phase_windForward(ip,el) - - integer, intent(in) :: & - ip, & !< integration point number - el !< element number - - integer :: & - co, & !< constituent number - so, ph, me - - - do co = 1,homogenization_Nconstituents(material_homogenizationAt(el)) - ph = material_phaseAt(co,el) - me = material_phaseMemberAt(co,ip,el) - - call mechanical_windForward(ph,me) - - if(damageState(ph)%sizeState > 0) damageState(ph)%state0(:,me) = damageState(ph)%state(:,me) - - - enddo - -end subroutine phase_windForward - - !-------------------------------------------------------------------------------------------------- !> @brief calculates orientations !-------------------------------------------------------------------------------------------------- @@ -629,11 +587,11 @@ function crystallite_push33ToRef(co,ce, tensor33) real(pReal), dimension(3,3) :: crystallite_push33ToRef real(pReal), dimension(3,3) :: T - integer :: ph, me + integer :: ph, en - ph = material_phaseAt2(co,ce) - me = material_phaseMemberAt2(co,ce) - T = matmul(material_orientation0(co,ph,me)%asMatrix(),transpose(math_inv33(phase_mechanical_getF(co,ce)))) ! ToDo: initial orientation correct? + ph = material_phaseID(co,ce) + en = material_phaseEntry(co,ce) + T = matmul(material_orientation0(co,ph,en)%asMatrix(),transpose(math_inv33(phase_mechanical_getF(co,ce)))) ! ToDo: initial orientation correct? crystallite_push33ToRef = matmul(transpose(T),matmul(tensor33,T)) @@ -658,8 +616,7 @@ end function converged !-------------------------------------------------------------------------------------------------- -!> @brief Write current restart information (Field and constitutive data) to file. -! ToDo: Merge data into one file for MPI +!> @brief Write restart data to file. !-------------------------------------------------------------------------------------------------- subroutine phase_restartWrite(fileHandle) @@ -687,8 +644,7 @@ end subroutine phase_restartWrite !-------------------------------------------------------------------------------------------------- -!> @brief Read data for restart -! ToDo: Merge data into one file for MPI +!> @brief Read restart data from file. !-------------------------------------------------------------------------------------------------- subroutine phase_restartRead(fileHandle) diff --git a/src/phase_damage.f90 b/src/phase_damage.f90 index c6b430b74..4397d5cad 100644 --- a/src/phase_damage.f90 +++ b/src/phase_damage.f90 @@ -1,7 +1,7 @@ !---------------------------------------------------------------------------------------------------- !> @brief internal microstructure state for all damage sources and kinematics constitutive models !---------------------------------------------------------------------------------------------------- -submodule(phase) damagee +submodule(phase) damage enum, bind(c); enumerator :: & DAMAGE_UNDEFINED_ID, & DAMAGE_ISOBRITTLE_ID, & @@ -151,7 +151,7 @@ module subroutine damage_init do ph = 1,phases%length - Nmembers = count(material_phaseAt2 == ph) + Nmembers = count(material_phaseID == ph) allocate(current(ph)%phi(Nmembers),source=1.0_pReal) allocate(current(ph)%d_phi_d_dot_phi(Nmembers),source=0.0_pReal) @@ -199,9 +199,9 @@ module subroutine phase_damage_getRateAndItsTangents(phiDot, dPhiDot_dPhi, phi, phiDot = 0.0_pReal dPhiDot_dPhi = 0.0_pReal - do co = 1, homogenization_Nconstituents(material_homogenizationAt2(ce)) - ph = material_phaseAt2(co,ce) - me = material_phasememberAt2(co,ce) + do co = 1, homogenization_Nconstituents(material_homogenizationID(ce)) + ph = material_phaseID(co,ce) + me = material_phaseEntry(co,ce) select case(phase_source(ph)) case (DAMAGE_ISOBRITTLE_ID) @@ -477,7 +477,7 @@ module subroutine phase_damage_set_phi(phi,co,ce) integer, intent(in) :: ce, co - current(material_phaseAt2(co,ce))%phi(material_phaseMemberAt2(co,ce)) = phi + current(material_phaseID(co,ce))%phi(material_phaseEntry(co,ce)) = phi end subroutine phase_damage_set_phi @@ -503,4 +503,21 @@ module function damage_phi(ph,me) result(phi) end function damage_phi -end submodule damagee + +!-------------------------------------------------------------------------------------------------- +!> @brief Forward data after successful increment. +!-------------------------------------------------------------------------------------------------- +module subroutine damage_forward() + + integer :: ph + + + do ph = 1, size(damageState) + if (damageState(ph)%sizeState > 0) & + damageState(ph)%state0 = damageState(ph)%state + enddo + +end subroutine damage_forward + + +end submodule damage diff --git a/src/phase_damage_anisobrittle.f90 b/src/phase_damage_anisobrittle.f90 index 3c573f33a..e9672dd3b 100644 --- a/src/phase_damage_anisobrittle.f90 +++ b/src/phase_damage_anisobrittle.f90 @@ -4,7 +4,7 @@ !> @brief material subroutine incorporating anisotropic brittle damage source mechanism !> @details to be done !-------------------------------------------------------------------------------------------------- -submodule (phase:damagee) anisobrittle +submodule (phase:damage) anisobrittle type :: tParameters !< container type for internal constitutive parameters real(pReal) :: & @@ -197,7 +197,7 @@ end subroutine anisobrittle_results !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- -module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) +module subroutine damage_anisobrittle_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) integer, intent(in) :: & ph,me @@ -253,6 +253,6 @@ module subroutine kinematics_cleavage_opening_LiAndItsTangent(Ld, dLd_dTstar, S, enddo end associate -end subroutine kinematics_cleavage_opening_LiAndItsTangent +end subroutine damage_anisobrittle_LiAndItsTangent end submodule anisobrittle diff --git a/src/phase_damage_anisoductile.f90 b/src/phase_damage_anisoductile.f90 index f047a5dac..54b63278f 100644 --- a/src/phase_damage_anisoductile.f90 +++ b/src/phase_damage_anisoductile.f90 @@ -4,7 +4,7 @@ !> @brief material subroutine incorporating anisotropic ductile damage source mechanism !> @details to be done !-------------------------------------------------------------------------------------------------- -submodule(phase:damagee) anisoductile +submodule(phase:damage) anisoductile type :: tParameters !< container type for internal constitutive parameters real(pReal) :: & @@ -35,7 +35,7 @@ module function anisoductile_init() result(mySources) pl, & sources, & src - integer :: Ninstances,Nmembers,p + integer :: Ninstances,Nmembers,ph integer, dimension(:), allocatable :: N_sl character(len=pStringLen) :: extmsg = '' @@ -50,15 +50,15 @@ module function anisoductile_init() result(mySources) phases => config_material%get('phase') allocate(param(phases%length)) - do p = 1, phases%length - if(mySources(p)) then - phase => phases%get(p) + do ph = 1, phases%length + if(mySources(ph)) then + phase => phases%get(ph) mech => phase%get('mechanical') pl => mech%get('plastic') sources => phase%get('damage') - associate(prm => param(p)) + associate(prm => param(ph)) src => sources%get(1) N_sl = pl%get_as1dInt('N_sl',defaultVal=emptyIntArray) @@ -78,10 +78,10 @@ module function anisoductile_init() result(mySources) if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q' if (any(prm%gamma_crit < 0.0_pReal)) extmsg = trim(extmsg)//' gamma_crit' - Nmembers=count(material_phaseAt2==p) - call phase_allocateState(damageState(p),Nmembers,1,1,0) - damageState(p)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal) - if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol' + Nmembers=count(material_phaseID==ph) + call phase_allocateState(damageState(ph),Nmembers,1,1,0) + damageState(ph)%atol = src%get_asFloat('anisoDuctile_atol',defaultVal=1.0e-3_pReal) + if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' anisoductile_atol' end associate diff --git a/src/phase_damage_isobrittle.f90 b/src/phase_damage_isobrittle.f90 index 3c7866e66..59cedb554 100644 --- a/src/phase_damage_isobrittle.f90 +++ b/src/phase_damage_isobrittle.f90 @@ -4,7 +4,7 @@ !> @brief material subroutine incoprorating isotropic brittle damage source mechanism !> @details to be done !-------------------------------------------------------------------------------------------------- -submodule(phase:damagee) isobrittle +submodule(phase:damage) isobrittle type :: tParameters !< container type for internal constitutive parameters real(pReal) :: & @@ -31,7 +31,7 @@ module function isobrittle_init() result(mySources) phase, & sources, & src - integer :: Nmembers,p + integer :: Nmembers,ph character(len=pStringLen) :: extmsg = '' @@ -45,12 +45,12 @@ module function isobrittle_init() result(mySources) phases => config_material%get('phase') allocate(param(phases%length)) - do p = 1, phases%length - if(mySources(p)) then - phase => phases%get(p) + do ph = 1, phases%length + if(mySources(ph)) then + phase => phases%get(ph) sources => phase%get('damage') - associate(prm => param(p)) + associate(prm => param(ph)) src => sources%get(1) prm%W_crit = src%get_asFloat('W_crit') @@ -64,10 +64,10 @@ module function isobrittle_init() result(mySources) ! sanity checks if (prm%W_crit <= 0.0_pReal) extmsg = trim(extmsg)//' W_crit' - Nmembers = count(material_phaseAt2==p) - call phase_allocateState(damageState(p),Nmembers,1,1,1) - damageState(p)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal) - if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol' + Nmembers = count(material_phaseID==ph) + call phase_allocateState(damageState(ph),Nmembers,1,1,1) + damageState(ph)%atol = src%get_asFloat('isoBrittle_atol',defaultVal=1.0e-3_pReal) + if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isobrittle_atol' end associate diff --git a/src/phase_damage_isoductile.f90 b/src/phase_damage_isoductile.f90 index b5ca9a1c2..1f1bba847 100644 --- a/src/phase_damage_isoductile.f90 +++ b/src/phase_damage_isoductile.f90 @@ -4,7 +4,7 @@ !> @brief material subroutine incorporating isotropic ductile damage source mechanism !> @details to be done !-------------------------------------------------------------------------------------------------- -submodule(phase:damagee) isoductile +submodule(phase:damage) isoductile type:: tParameters !< container type for internal constitutive parameters real(pReal) :: & @@ -33,7 +33,7 @@ module function isoductile_init() result(mySources) phase, & sources, & src - integer :: Ninstances,Nmembers,p + integer :: Ninstances,Nmembers,ph character(len=pStringLen) :: extmsg = '' @@ -47,12 +47,12 @@ module function isoductile_init() result(mySources) phases => config_material%get('phase') allocate(param(phases%length)) - do p = 1, phases%length - if(mySources(p)) then - phase => phases%get(p) + do ph = 1, phases%length + if(mySources(ph)) then + phase => phases%get(ph) sources => phase%get('damage') - associate(prm => param(p)) + associate(prm => param(ph)) src => sources%get(1) prm%q = src%get_asFloat('q') @@ -68,10 +68,10 @@ module function isoductile_init() result(mySources) if (prm%q <= 0.0_pReal) extmsg = trim(extmsg)//' q' if (prm%gamma_crit <= 0.0_pReal) extmsg = trim(extmsg)//' gamma_crit' - Nmembers=count(material_phaseAt2==p) - call phase_allocateState(damageState(p),Nmembers,1,1,0) - damageState(p)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal) - if(any(damageState(p)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol' + Nmembers=count(material_phaseID==ph) + call phase_allocateState(damageState(ph),Nmembers,1,1,0) + damageState(ph)%atol = src%get_asFloat('isoDuctile_atol',defaultVal=1.0e-3_pReal) + if(any(damageState(ph)%atol < 0.0_pReal)) extmsg = trim(extmsg)//' isoductile_atol' end associate diff --git a/src/phase_mechanical.f90 b/src/phase_mechanical.f90 index 99bde64a5..ddd6e1fcb 100644 --- a/src/phase_mechanical.f90 +++ b/src/phase_mechanical.f90 @@ -1060,26 +1060,6 @@ subroutine crystallite_results(group,ph) end subroutine crystallite_results -!-------------------------------------------------------------------------------------------------- -!> @brief Wind homog inc forward. -!-------------------------------------------------------------------------------------------------- -module subroutine mechanical_windForward(ph,me) - - integer, intent(in) :: ph, me - - - phase_mechanical_Fp0(ph)%data(1:3,1:3,me) = phase_mechanical_Fp(ph)%data(1:3,1:3,me) - phase_mechanical_Fi0(ph)%data(1:3,1:3,me) = phase_mechanical_Fi(ph)%data(1:3,1:3,me) - phase_mechanical_F0(ph)%data(1:3,1:3,me) = phase_mechanical_F(ph)%data(1:3,1:3,me) - phase_mechanical_Li0(ph)%data(1:3,1:3,me) = phase_mechanical_Li(ph)%data(1:3,1:3,me) - phase_mechanical_Lp0(ph)%data(1:3,1:3,me) = phase_mechanical_Lp(ph)%data(1:3,1:3,me) - phase_mechanical_S0(ph)%data(1:3,1:3,me) = phase_mechanical_S(ph)%data(1:3,1:3,me) - - plasticState(ph)%State0(:,me) = plasticState(ph)%state(:,me) - -end subroutine mechanical_windForward - - !-------------------------------------------------------------------------------------------------- !> @brief Forward data after successful increment. ! ToDo: Any guessing for the current states possible? @@ -1235,9 +1215,9 @@ module subroutine mechanical_restore(ce,includeL) co, ph, me - do co = 1,homogenization_Nconstituents(material_homogenizationAt2(ce)) - ph = material_phaseAt2(co,ce) - me = material_phaseMemberAt2(co,ce) + do co = 1,homogenization_Nconstituents(material_homogenizationID(ce)) + ph = material_phaseID(co,ce) + me = material_phaseEntry(co,ce) if (includeL) then phase_mechanical_Lp(ph)%data(1:3,1:3,me) = phase_mechanical_Lp0(ph)%data(1:3,1:3,me) phase_mechanical_Li(ph)%data(1:3,1:3,me) = phase_mechanical_Li0(ph)%data(1:3,1:3,me) @@ -1285,8 +1265,8 @@ module function phase_mechanical_dPdF(dt,co,ce) result(dPdF) logical :: error - ph = material_phaseAt2(co,ce) - me = material_phaseMemberAt2(co,ce) + ph = material_phaseID(co,ce) + me = material_phaseEntry(co,ce) call phase_hooke_SandItsTangents(devNull,dSdFe,dSdFi, & phase_mechanical_Fe(ph)%data(1:3,1:3,me), & @@ -1450,7 +1430,7 @@ module function phase_mechanical_getF(co,ce) result(F) real(pReal), dimension(3,3) :: F - F = phase_mechanical_F(material_phaseAt2(co,ce))%data(1:3,1:3,material_phaseMemberAt2(co,ce)) + F = phase_mechanical_F(material_phaseID(co,ce))%data(1:3,1:3,material_phaseEntry(co,ce)) end function phase_mechanical_getF @@ -1479,7 +1459,7 @@ module function phase_mechanical_getP(co,ce) result(P) real(pReal), dimension(3,3) :: P - P = phase_mechanical_P(material_phaseAt2(co,ce))%data(1:3,1:3,material_phaseMemberAt2(co,ce)) + P = phase_mechanical_P(material_phaseID(co,ce))%data(1:3,1:3,material_phaseEntry(co,ce)) end function phase_mechanical_getP @@ -1491,7 +1471,7 @@ module subroutine phase_mechanical_setF(F,co,ce) integer, intent(in) :: co, ce - phase_mechanical_F(material_phaseAt2(co,ce))%data(1:3,1:3,material_phaseMemberAt2(co,ce)) = F + phase_mechanical_F(material_phaseID(co,ce))%data(1:3,1:3,material_phaseEntry(co,ce)) = F end subroutine phase_mechanical_setF diff --git a/src/phase_mechanical_eigen.f90 b/src/phase_mechanical_eigen.f90 index 955f6ca5d..3f2aec58c 100644 --- a/src/phase_mechanical_eigen.f90 +++ b/src/phase_mechanical_eigen.f90 @@ -9,13 +9,13 @@ submodule(phase:mechanical) eigen model_damage interface - module function kinematics_cleavage_opening_init() result(myKinematics) + module function damage_anisobrittle_init() result(myKinematics) logical, dimension(:), allocatable :: myKinematics - end function kinematics_cleavage_opening_init + end function damage_anisobrittle_init - module function kinematics_slipplane_opening_init() result(myKinematics) + module function damage_isoductile_init() result(myKinematics) logical, dimension(:), allocatable :: myKinematics - end function kinematics_slipplane_opening_init + end function damage_isoductile_init module function thermalexpansion_init(kinematics_length) result(myKinematics) integer, intent(in) :: kinematics_length @@ -70,8 +70,8 @@ module subroutine eigendeformation_init(phases) allocate(model_damage(phases%length), source = KINEMATICS_UNDEFINED_ID) - where(kinematics_cleavage_opening_init()) model_damage = KINEMATICS_cleavage_opening_ID - where(kinematics_slipplane_opening_init()) model_damage = KINEMATICS_slipplane_opening_ID + where(damage_anisobrittle_init()) model_damage = KINEMATICS_cleavage_opening_ID + where(damage_isoductile_init()) model_damage = KINEMATICS_slipplane_opening_ID end subroutine eigendeformation_init @@ -198,12 +198,12 @@ module subroutine phase_LiAndItsTangents(Li, dLi_dS, dLi_dFi, & select case (model_damage(ph)) case (KINEMATICS_cleavage_opening_ID) - call kinematics_cleavage_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, me) + call damage_anisobrittle_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, me) Li = Li + my_Li dLi_dS = dLi_dS + my_dLi_dS active = .true. case (KINEMATICS_slipplane_opening_ID) - call kinematics_slipplane_opening_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, me) + call damage_isoductile_LiAndItsTangent(my_Li, my_dLi_dS, S, ph, me) Li = Li + my_Li dLi_dS = dLi_dS + my_dLi_dS active = .true. diff --git a/src/phase_mechanical_eigen_cleavageopening.f90 b/src/phase_mechanical_eigen_cleavageopening.f90 index bccf1cb5f..4243d09a4 100644 --- a/src/phase_mechanical_eigen_cleavageopening.f90 +++ b/src/phase_mechanical_eigen_cleavageopening.f90 @@ -13,7 +13,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -module function kinematics_cleavage_opening_init() result(myKinematics) +module function damage_anisobrittle_init() result(myKinematics) logical, dimension(:), allocatable :: myKinematics @@ -24,7 +24,7 @@ module function kinematics_cleavage_opening_init() result(myKinematics) print'(/,a)', ' <<<+- phase:mechanical:eigen:cleavageopening init -+>>>' print'(a,i2)', ' # phases: ',count(myKinematics); flush(IO_STDOUT) -end function kinematics_cleavage_opening_init +end function damage_anisobrittle_init end submodule cleavageopening diff --git a/src/phase_mechanical_eigen_slipplaneopening.f90 b/src/phase_mechanical_eigen_slipplaneopening.f90 index 18d99f750..2fd14700d 100644 --- a/src/phase_mechanical_eigen_slipplaneopening.f90 +++ b/src/phase_mechanical_eigen_slipplaneopening.f90 @@ -6,7 +6,7 @@ !-------------------------------------------------------------------------------------------------- submodule(phase:eigen) slipplaneopening - integer, dimension(:), allocatable :: kinematics_slipplane_opening_instance + integer, dimension(:), allocatable :: damage_isoductile_instance type :: tParameters !< container type for internal constitutive parameters integer :: & @@ -32,7 +32,7 @@ contains !> @brief module initialization !> @details reads in material parameters, allocates arrays, and does sanity checks !-------------------------------------------------------------------------------------------------- -module function kinematics_slipplane_opening_init() result(myKinematics) +module function damage_isoductile_init() result(myKinematics) logical, dimension(:), allocatable :: myKinematics @@ -107,13 +107,13 @@ module function kinematics_slipplane_opening_init() result(myKinematics) enddo -end function kinematics_slipplane_opening_init +end function damage_isoductile_init !-------------------------------------------------------------------------------------------------- !> @brief contains the constitutive equation for calculating the velocity gradient !-------------------------------------------------------------------------------------------------- -module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) +module subroutine damage_isoductile_LiAndItsTangent(Ld, dLd_dTstar, S, ph,me) integer, intent(in) :: & ph, me @@ -179,6 +179,6 @@ module subroutine kinematics_slipplane_opening_LiAndItsTangent(Ld, dLd_dTstar, S end associate -end subroutine kinematics_slipplane_opening_LiAndItsTangent +end subroutine damage_isoductile_LiAndItsTangent end submodule slipplaneopening diff --git a/src/phase_mechanical_plastic_dislotungsten.f90 b/src/phase_mechanical_plastic_dislotungsten.f90 index b7ade48a8..2e7130756 100644 --- a/src/phase_mechanical_plastic_dislotungsten.f90 +++ b/src/phase_mechanical_plastic_dislotungsten.f90 @@ -220,7 +220,7 @@ module function plastic_dislotungsten_init() result(myPlasticity) !-------------------------------------------------------------------------------------------------- ! allocate state arrays - Nmembers = count(material_phaseAt2 == ph) + Nmembers = count(material_phaseID == ph) sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl sizeState = sizeDotState diff --git a/src/phase_mechanical_plastic_dislotwin.f90 b/src/phase_mechanical_plastic_dislotwin.f90 index 1e67db59e..43be614ac 100644 --- a/src/phase_mechanical_plastic_dislotwin.f90 +++ b/src/phase_mechanical_plastic_dislotwin.f90 @@ -406,7 +406,7 @@ module function plastic_dislotwin_init() result(myPlasticity) !-------------------------------------------------------------------------------------------------- ! allocate state arrays - Nmembers = count(material_phaseAt2 == ph) + Nmembers = count(material_phaseID == ph) sizeDotState = size(['rho_mob ','rho_dip ','gamma_sl']) * prm%sum_N_sl & + size(['f_tw']) * prm%sum_N_tw & + size(['f_tr']) * prm%sum_N_tr diff --git a/src/phase_mechanical_plastic_isotropic.f90 b/src/phase_mechanical_plastic_isotropic.f90 index 5a37369f0..548e2fd75 100644 --- a/src/phase_mechanical_plastic_isotropic.f90 +++ b/src/phase_mechanical_plastic_isotropic.f90 @@ -119,7 +119,7 @@ module function plastic_isotropic_init() result(myPlasticity) !-------------------------------------------------------------------------------------------------- ! allocate state arrays - Nmembers = count(material_phaseAt2 == ph) + Nmembers = count(material_phaseID == ph) sizeDotState = size(['xi ','gamma']) sizeState = sizeDotState diff --git a/src/phase_mechanical_plastic_kinehardening.f90 b/src/phase_mechanical_plastic_kinehardening.f90 index 0c02e11db..936a7f884 100644 --- a/src/phase_mechanical_plastic_kinehardening.f90 +++ b/src/phase_mechanical_plastic_kinehardening.f90 @@ -165,7 +165,7 @@ module function plastic_kinehardening_init() result(myPlasticity) !-------------------------------------------------------------------------------------------------- ! allocate state arrays - Nmembers = count(material_phaseAt2 == ph) + Nmembers = count(material_phaseID == ph) sizeDotState = size(['crss ','crss_back', 'accshear ']) * prm%sum_N_sl !ToDo: adjust names like in material.yaml sizeDeltaState = size(['sense ', 'chi0 ', 'gamma0' ]) * prm%sum_N_sl !ToDo: adjust names like in material.yaml sizeState = sizeDotState + sizeDeltaState diff --git a/src/phase_mechanical_plastic_none.f90 b/src/phase_mechanical_plastic_none.f90 index 28e9fbc7c..544de31b1 100644 --- a/src/phase_mechanical_plastic_none.f90 +++ b/src/phase_mechanical_plastic_none.f90 @@ -30,7 +30,7 @@ module function plastic_none_init() result(myPlasticity) phases => config_material%get('phase') do ph = 1, phases%length if(.not. myPlasticity(ph)) cycle - call phase_allocateState(plasticState(ph),count(material_phaseAt2 == ph),0,0,0) + call phase_allocateState(plasticState(ph),count(material_phaseID == ph),0,0,0) enddo end function plastic_none_init diff --git a/src/phase_mechanical_plastic_nonlocal.f90 b/src/phase_mechanical_plastic_nonlocal.f90 index d86f7dc7e..ba7f76881 100644 --- a/src/phase_mechanical_plastic_nonlocal.f90 +++ b/src/phase_mechanical_plastic_nonlocal.f90 @@ -398,7 +398,7 @@ module function plastic_nonlocal_init() result(myPlasticity) !-------------------------------------------------------------------------------------------------- ! allocate state arrays - Nmembers = count(material_phaseAt2 == ph) + Nmembers = count(material_phaseID == ph) sizeDotState = size([ 'rhoSglEdgePosMobile ','rhoSglEdgeNegMobile ', & 'rhoSglScrewPosMobile ','rhoSglScrewNegMobile ', & 'rhoSglEdgePosImmobile ','rhoSglEdgeNegImmobile ', & @@ -527,7 +527,7 @@ module function plastic_nonlocal_init() result(myPlasticity) if(.not. myPlasticity(ph)) cycle phase => phases%get(ph) - Nmembers = count(material_phaseAt2 == ph) + Nmembers = count(material_phaseID == ph) l = 0 do t = 1,4 do s = 1,param(ph)%sum_N_sl @@ -1824,9 +1824,9 @@ subroutine storeGeometry(ph) V = reshape(IPvolume,[product(shape(IPvolume))]) - do ce = 1, size(material_homogenizationMemberAt2,1) + do ce = 1, size(material_homogenizationEntry,1) do co = 1, homogenization_maxNconstituents - if (material_phaseAt2(co,ce) == ph) geom(ph)%V_0(material_phaseMemberAt2(co,ce)) = V(ce) + if (material_phaseID(co,ce) == ph) geom(ph)%V_0(material_phaseEntry(co,ce)) = V(ce) enddo enddo diff --git a/src/phase_mechanical_plastic_phenopowerlaw.f90 b/src/phase_mechanical_plastic_phenopowerlaw.f90 index 6be46f6a1..77c47907f 100644 --- a/src/phase_mechanical_plastic_phenopowerlaw.f90 +++ b/src/phase_mechanical_plastic_phenopowerlaw.f90 @@ -223,7 +223,7 @@ module function plastic_phenopowerlaw_init() result(myPlasticity) !-------------------------------------------------------------------------------------------------- ! allocate state arrays - Nmembers = count(material_phaseAt2 == ph) + Nmembers = count(material_phaseID == ph) sizeDotState = size(['xi_sl ','gamma_sl']) * prm%sum_N_sl & + size(['xi_tw ','gamma_tw']) * prm%sum_N_tw sizeState = sizeDotState diff --git a/src/phase_thermal.f90 b/src/phase_thermal.f90 index 21ec93c9c..1e157f338 100644 --- a/src/phase_thermal.f90 +++ b/src/phase_thermal.f90 @@ -89,7 +89,7 @@ module subroutine thermal_init(phases) allocate(thermal_Nsources(phases%length),source = 0) do ph = 1, phases%length - Nmembers = count(material_phaseAt2 == ph) + Nmembers = count(material_phaseID == ph) allocate(current(ph)%T(Nmembers),source=300.0_pReal) allocate(current(ph)%dot_T(Nmembers),source=0.0_pReal) phase => phases%get(ph) @@ -268,8 +268,8 @@ module subroutine phase_thermal_setField(T,dot_T, co,ce) integer, intent(in) :: ce, co - current(material_phaseAt2(co,ce))%T(material_phaseMemberAt2(co,ce)) = T - current(material_phaseAt2(co,ce))%dot_T(material_phaseMemberAt2(co,ce)) = dot_T + current(material_phaseID(co,ce))%T(material_phaseEntry(co,ce)) = T + current(material_phaseID(co,ce))%dot_T(material_phaseEntry(co,ce)) = dot_T end subroutine phase_thermal_setField diff --git a/src/phase_thermal_dissipation.f90 b/src/phase_thermal_dissipation.f90 index d3e7094a1..b16ed3cf7 100644 --- a/src/phase_thermal_dissipation.f90 +++ b/src/phase_thermal_dissipation.f90 @@ -54,7 +54,7 @@ module function dissipation_init(source_length) result(mySources) src => sources%get(so) prm%kappa = src%get_asFloat('kappa') - Nmembers = count(material_phaseAt2 == ph) + Nmembers = count(material_phaseID == ph) call phase_allocateState(thermalState(ph)%p(so),Nmembers,0,0,0) end associate diff --git a/src/phase_thermal_externalheat.f90 b/src/phase_thermal_externalheat.f90 index 760326044..d5bbc7c38 100644 --- a/src/phase_thermal_externalheat.f90 +++ b/src/phase_thermal_externalheat.f90 @@ -67,7 +67,7 @@ module function externalheat_init(source_length) result(mySources) prm%f_T = src%get_as1dFloat('f_T',requiredSize = size(prm%t_n)) - Nmembers = count(material_phaseAt2 == ph) + Nmembers = count(material_phaseID == ph) call phase_allocateState(thermalState(ph)%p(so),Nmembers,1,1,0) end associate endif