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

Merge branch 'Fortran-cleaning' into 'development' 

shorter names

See merge request damask/DAMASK!364
This commit is contained in:
Sharan Roongta 2021-04-07 08:41:46 +00:00
parent 6bc1d5ad98 5a361c12f8
commit 3c967abcde
32 changed files with 347 additions and 332 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
src/commercialFEM_fileList.f90
View File

@ -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"

14
src/homogenization.f90
View File

@ -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 = &

37
src/homogenization_damage.f90
View File

@ -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

14
src/homogenization_damage_pass.f90 Normal file
View File

@ -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

84
src/homogenization_mechanical.f90
View File

@ -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

114
src/homogenization_mechanical_RGC.f90
View File

@ -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

14
src/homogenization_mechanical_isostrain.f90
View File

@ -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

10
src/homogenization_mechanical_pass.f90
View File

@ -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

44
src/homogenization_thermal.f90
View File

@ -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

14
src/homogenization_thermal_isotemperature.f90 Normal file
View File

@ -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

14
src/homogenization_thermal_pass.f90 Normal file
View File

@ -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

24
src/material.f90
View File

@ -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

100
src/phase.f90
View File

@ -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)

31
src/phase_damage.f90
View File

@ -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

6
src/phase_damage_anisobrittle.f90
View File

@ -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

20
src/phase_damage_anisoductile.f90
View File

@ -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

20
src/phase_damage_isobrittle.f90
View File

@ -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

20
src/phase_damage_isoductile.f90
View File

@ -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

36
src/phase_mechanical.f90
View File

@ -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

16
src/phase_mechanical_eigen.f90
View File

@ -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.

4
src/phase_mechanical_eigen_cleavageopening.f90
View File

@ -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

10
src/phase_mechanical_eigen_slipplaneopening.f90
View File

@ -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

2
src/phase_mechanical_plastic_dislotungsten.f90
View File

@ -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

2
src/phase_mechanical_plastic_dislotwin.f90
View File

@ -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

2
src/phase_mechanical_plastic_isotropic.f90
View File

@ -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

2
src/phase_mechanical_plastic_kinehardening.f90
View File

@ -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

2
src/phase_mechanical_plastic_none.f90
View File

@ -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

8
src/phase_mechanical_plastic_nonlocal.f90
View File

@ -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

2
src/phase_mechanical_plastic_phenopowerlaw.f90
View File

@ -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

6
src/phase_thermal.f90
View File

@ -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

2
src/phase_thermal_dissipation.f90
View File

@ -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

2
src/phase_thermal_externalheat.f90
View File

@ -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