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

1720 lines
92 KiB
Fortran
Raw Blame History

!--------------------------------------------------------------------------------------------------
! $Id$
!--------------------------------------------------------------------------------------------------
!> @author Denny Tjahjanto, Max-Planck-Institut für Eisenforschung GmbH
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief Relaxed grain cluster (RGC) homogenization scheme
!> Ngrains is defined as p x q x r (cluster)
!--------------------------------------------------------------------------------------------------
module homogenization_RGC
use prec, only: &
pReal, &
pInt
implicit none
private
integer(pInt), dimension(:), allocatable, public :: &
homogenization_RGC_sizeState, &
homogenization_RGC_sizePostResults
integer(pInt), dimension(:,:), allocatable,target, public :: &
homogenization_RGC_sizePostResult
character(len=64), dimension(:,:), allocatable,target, public :: &
homogenization_RGC_output ! name of each post result output
integer(pInt), dimension(:,:), allocatable, private :: &
homogenization_RGC_Ngrains
real(pReal), dimension(:,:), allocatable, private :: &
homogenization_RGC_dAlpha, &
homogenization_RGC_angles
real(pReal), dimension(:,:,:,:), allocatable, private :: &
homogenization_RGC_orientation
real(pReal), dimension(:), allocatable, private :: &
homogenization_RGC_xiAlpha, &
homogenization_RGC_ciAlpha
enum, bind(c)
enumerator :: undefined_ID, &
temperature_ID, &
constitutivework_ID, &
penaltyenergy_ID, &
volumediscrepancy_ID, &
averagerelaxrate_ID,&
maximumrelaxrate_ID,&
ipcoords_ID,&
magnitudemismatch_ID,&
avgdefgrad_ID,&
avgfirstpiola_ID
end enum
integer(kind(undefined_ID)), dimension(:,:), allocatable, private :: &
homogenization_RGC_outputID !< ID of each post result output
public :: &
homogenization_RGC_init, &
homogenization_RGC_partitionDeformation, &
homogenization_RGC_averageStressAndItsTangent, &
homogenization_RGC_updateState, &
homogenization_RGC_postResults
private :: &
homogenization_RGC_stressPenalty, &
homogenization_RGC_volumePenalty, &
homogenization_RGC_grainDeformation, &
homogenization_RGC_surfaceCorrection, &
homogenization_RGC_equivalentModuli, &
homogenization_RGC_relaxationVector, &
homogenization_RGC_interfaceNormal, &
homogenization_RGC_getInterface, &
homogenization_RGC_grain1to3, &
homogenization_RGC_grain3to1, &
homogenization_RGC_interface4to1, &
homogenization_RGC_interface1to4
contains
!--------------------------------------------------------------------------------------------------
!> @brief allocates all neccessary fields, reads information from material configuration file
!--------------------------------------------------------------------------------------------------
subroutine homogenization_RGC_init(fileUnit)
use, intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use prec, only: &
pReal, &
pInt
use debug, only: &
debug_level, &
debug_homogenization, &
debug_levelBasic, &
debug_levelExtensive
use math, only: &
math_Mandel3333to66,&
math_Voigt66to3333, &
math_I3, &
math_sampleRandomOri,&
math_EulerToR,&
INRAD
use mesh, only: &
mesh_maxNips, &
mesh_NcpElems,&
mesh_element, &
FE_Nips, &
FE_geomtype
use IO
use material
implicit none
integer(pInt), intent(in) :: fileUnit !< file pointer to material configuration
integer(pInt), parameter :: MAXNCHUNKS = 4_pInt
integer(pInt), dimension(1_pInt+2_pInt*MAXNCHUNKS) :: positions
#ifdef NEWSTATE
integer :: &
homog, &
NofMyHomog, &
instance, &
sizeHState
#endif
integer(pInt) :: section=0_pInt, maxNinstance, i,j,e, output=-1_pInt, mySize, myInstance
character(len=65536) :: &
tag = '', &
line = ''
write(6,'(/,a)') ' <<<+- homogenization_'//HOMOGENIZATION_RGC_label//' init -+>>>'
write(6,'(a)') ' $Id$'
write(6,'(a15,a)') ' Current time: ',IO_timeStamp()
#include "compilation_info.f90"
maxNinstance = int(count(homogenization_type == HOMOGENIZATION_RGC_ID),pInt)
if (maxNinstance == 0_pInt) return
allocate(homogenization_RGC_sizeState(maxNinstance), source=0_pInt)
allocate(homogenization_RGC_sizePostResults(maxNinstance), source=0_pInt)
allocate(homogenization_RGC_Ngrains(3,maxNinstance), source=0_pInt)
allocate(homogenization_RGC_ciAlpha(maxNinstance), source=0.0_pReal)
allocate(homogenization_RGC_xiAlpha(maxNinstance), source=0.0_pReal)
allocate(homogenization_RGC_dAlpha(3,maxNinstance), source=0.0_pReal)
allocate(homogenization_RGC_angles(3,maxNinstance), source=400.0_pReal)
allocate(homogenization_RGC_output(maxval(homogenization_Noutput),maxNinstance))
homogenization_RGC_output=''
allocate(homogenization_RGC_outputID(maxval(homogenization_Noutput),maxNinstance),source=undefined_ID)
allocate(homogenization_RGC_sizePostResult(maxval(homogenization_Noutput),maxNinstance),&
source=0_pInt)
allocate(homogenization_RGC_orientation(3,3,mesh_maxNips,mesh_NcpElems), source=0.0_pReal)
homogenization_RGC_orientation = spread(spread(math_I3,3,mesh_maxNips),4,mesh_NcpElems) ! initialize to identity
rewind(fileUnit)
do while (trim(line) /= IO_EOF .and. IO_lc(IO_getTag(line,'<','>'))/=material_partHomogenization) ! wind forward to <homogenization>
line = IO_read(fileUnit)
enddo
do while (trim(line) /= IO_EOF) ! read through sections of homogenization part
line = IO_read(fileUnit)
if (IO_isBlank(line)) cycle ! skip empty lines
if (IO_getTag(line,'<','>') /= '') then ! stop at next part
line = IO_read(fileUnit, .true.) ! reset IO_read
exit
endif
if (IO_getTag(line,'[',']') /= '') then ! next section
section = section + 1_pInt
output = 0_pInt ! reset output counter
cycle
endif
if (section > 0_pInt ) then ! do not short-circuit here (.and. with next if-statement). It's not safe in Fortran
if (homogenization_type(section) == HOMOGENIZATION_RGC_ID) then ! one of my sections
i = homogenization_typeInstance(section) ! which instance of my type is present homogenization
positions = IO_stringPos(line,MAXNCHUNKS)
tag = IO_lc(IO_stringValue(line,positions,1_pInt)) ! extract key
select case(tag)
case ('nconstituents','ngrains','type')
case ('(output)')
output = output + 1_pInt
homogenization_RGC_output(output,i) = IO_lc(IO_stringValue(line,positions,2_pInt))
select case(homogenization_RGC_output(output,i))
case('constitutivework')
homogenization_RGC_outputID(output,i) = constitutivework_ID
case('penaltyenergy')
homogenization_RGC_outputID(output,i) = penaltyenergy_ID
case('volumediscrepancy')
homogenization_RGC_outputID(output,i) = volumediscrepancy_ID
case('averagerelaxrate')
homogenization_RGC_outputID(output,i) = averagerelaxrate_ID
case('maximumrelaxrate')
homogenization_RGC_outputID(output,i) = maximumrelaxrate_ID
case('magnitudemismatch')
homogenization_RGC_outputID(output,i) = magnitudemismatch_ID
case('temperature')
homogenization_RGC_outputID(output,i) = temperature_ID
case('ipcoords')
homogenization_RGC_outputID(output,i) = ipcoords_ID
case('avgdefgrad','avgf')
homogenization_RGC_outputID(output,i) = avgdefgrad_ID
case('avgp','avgfirstpiola','avg1stpiola')
homogenization_RGC_outputID(output,i) = avgfirstpiola_ID
case default
call IO_error(105_pInt,ext_msg=IO_stringValue(line,positions,2_pInt)//&
' ('//HOMOGENIZATION_RGC_label//')')
end select
case ('clustersize')
homogenization_RGC_Ngrains(1,i) = IO_intValue(line,positions,2_pInt)
homogenization_RGC_Ngrains(2,i) = IO_intValue(line,positions,3_pInt)
homogenization_RGC_Ngrains(3,i) = IO_intValue(line,positions,4_pInt)
if (homogenization_Ngrains(section) /= product(homogenization_RGC_Ngrains(1:3,i))) &
call IO_error(211_pInt,ext_msg=trim(tag)//' ('//HOMOGENIZATION_RGC_label//')')
case ('scalingparameter')
homogenization_RGC_xiAlpha(i) = IO_floatValue(line,positions,2_pInt)
case ('overproportionality')
homogenization_RGC_ciAlpha(i) = IO_floatValue(line,positions,2_pInt)
case ('grainsize')
homogenization_RGC_dAlpha(1,i) = IO_floatValue(line,positions,2_pInt)
homogenization_RGC_dAlpha(2,i) = IO_floatValue(line,positions,3_pInt)
homogenization_RGC_dAlpha(3,i) = IO_floatValue(line,positions,4_pInt)
case ('clusterorientation')
homogenization_RGC_angles(1,i) = IO_floatValue(line,positions,2_pInt)
homogenization_RGC_angles(2,i) = IO_floatValue(line,positions,3_pInt)
homogenization_RGC_angles(3,i) = IO_floatValue(line,positions,4_pInt)
case default
call IO_error(210_pInt,ext_msg=trim(tag)//' ('//HOMOGENIZATION_RGC_label//')')
end select
endif
endif
enddo
!--------------------------------------------------------------------------------------------------
! assigning cluster orientations
elementLooping: do e = 1_pInt,mesh_NcpElems
if (homogenization_type(mesh_element(3,e)) == HOMOGENIZATION_RGC_ID) then
myInstance = homogenization_typeInstance(mesh_element(3,e))
if (all (homogenization_RGC_angles(1:3,myInstance) >= 399.9_pReal)) then
homogenization_RGC_orientation(1:3,1:3,1,e) = math_EulerToR(math_sampleRandomOri())
do i = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,e)))
if (microstructure_elemhomo(mesh_element(4,e))) then
homogenization_RGC_orientation(1:3,1:3,i,e) = homogenization_RGC_orientation(1:3,1:3,1,e)
else
homogenization_RGC_orientation(1:3,1:3,i,e) = math_EulerToR(math_sampleRandomOri())
endif
enddo
else
do i = 1_pInt,FE_Nips(FE_geomtype(mesh_element(2,e)))
homogenization_RGC_orientation(1:3,1:3,i,e) = &
math_EulerToR(homogenization_RGC_angles(1:3,myInstance)*inRad)
enddo
endif
endif
enddo elementLooping
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt) then
do i = 1_pInt,maxNinstance
write(6,'(a15,1x,i4,/)') 'instance: ', i
write(6,'(a25,3(1x,i8))') 'cluster size: ',(homogenization_RGC_Ngrains(j,i),j=1_pInt,3_pInt)
write(6,'(a25,1x,e10.3)') 'scaling parameter: ', homogenization_RGC_xiAlpha(i)
write(6,'(a25,1x,e10.3)') 'over-proportionality: ', homogenization_RGC_ciAlpha(i)
write(6,'(a25,3(1x,e10.3))') 'grain size: ',(homogenization_RGC_dAlpha(j,i),j=1_pInt,3_pInt)
write(6,'(a25,3(1x,e10.3))') 'cluster orientation: ',(homogenization_RGC_angles(j,i),j=1_pInt,3_pInt)
enddo
endif
do i = 1_pInt,maxNinstance
do j = 1_pInt,maxval(homogenization_Noutput)
select case(homogenization_RGC_outputID(j,i))
case(temperature_ID,constitutivework_ID,penaltyenergy_ID,volumediscrepancy_ID, &
averagerelaxrate_ID,maximumrelaxrate_ID)
mySize = 1_pInt
case(ipcoords_ID,magnitudemismatch_ID)
mySize = 3_pInt
case(avgdefgrad_ID,avgfirstpiola_ID)
mySize = 9_pInt
case default
mySize = 0_pInt
end select
outputFound: if (mySize > 0_pInt) then
homogenization_RGC_sizePostResult(j,i) = mySize
homogenization_RGC_sizePostResults(i) = &
homogenization_RGC_sizePostResults(i) + mySize
endif outputFound
enddo
homogenization_RGC_sizeState(i) &
= 3_pInt*(homogenization_RGC_Ngrains(1,i)-1_pInt)*homogenization_RGC_Ngrains(2,i)*homogenization_RGC_Ngrains(3,i) &
+ 3_pInt*homogenization_RGC_Ngrains(1,i)*(homogenization_RGC_Ngrains(2,i)-1_pInt)*homogenization_RGC_Ngrains(3,i) &
+ 3_pInt*homogenization_RGC_Ngrains(1,i)*homogenization_RGC_Ngrains(2,i)*(homogenization_RGC_Ngrains(3,i)-1_pInt) &
+ 8_pInt ! (1) Average constitutive work, (2-4) Overall mismatch, (5) Average penalty energy,
! (6) Volume discrepancy, (7) Avg relaxation rate component, (8) Max relaxation rate component
enddo
#ifdef NEWSTATE
initializeInstances: do homog = 1_pInt, material_Nhomogenization
myhomog: if (homogenization_type(homog) == HOMOGENIZATION_RGC_ID) then
NofMyHomog = count(material_homog == homog)
! instance = phase_plasticityInstance(phase)
! allocate homogenization state arrays
sizeHState = homogenization_RGC_sizeState(homog)
homogState(homog)%sizeState = sizeHState
homogState(homog)%sizePostResults = homogenization_RGC_sizePostResults(homog)
allocate(homogState(homog)%state0 ( sizeHState,NofMyHomog), source=0.0_pReal)
allocate(homogState(homog)%subState0 ( sizeHState,NofMyHomog), source=0.0_pReal)
allocate(homogState(homog)%state ( sizeHState,NofMyHomog), source=0.0_pReal)
endif myhomog
enddo initializeInstances
#endif
end subroutine homogenization_RGC_init
!--------------------------------------------------------------------------------------------------
!> @brief partitions the deformation gradient onto the constituents
!--------------------------------------------------------------------------------------------------
#ifdef NEWSTATE
subroutine homogenization_RGC_partitionDeformation(F,avgF,ip,el)
#else
subroutine homogenization_RGC_partitionDeformation(F,avgF,state,ip,el)
#endif
use prec, only: &
p_vec
use debug, only: &
debug_level, &
debug_homogenization, &
debug_levelExtensive
use mesh, only: &
mesh_element
use material, only: &
homogenization_maxNgrains, &
homogenization_Ngrains,&
#ifdef NEWSTATE
homogState, &
mappingHomogenization, &
#endif
homogenization_typeInstance
use FEsolving, only: &
theInc,&
cycleCounter
implicit none
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: F !< partioned F per grain
real(pReal), dimension (3,3), intent(in) :: avgF !< averaged F
#ifndef NEWSTATE
type(p_vec), intent(in) :: state
#endif
integer(pInt), intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), dimension (3) :: aVect,nVect
integer(pInt), dimension (4) :: intFace
integer(pInt), dimension (3) :: iGrain3
integer(pInt) :: homID, iGrain,iFace,i,j
integer(pInt), parameter :: nFace = 6_pInt
!--------------------------------------------------------------------------------------------------
! debugging the overall deformation gradient
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a,i3,a,i3,a)')'========== Increment: ',theInc,' Cycle: ',cycleCounter,' =========='
write(6,'(1x,a32)')'Overall deformation gradient: '
do i = 1_pInt,3_pInt
write(6,'(1x,3(e15.8,1x))')(avgF(i,j), j = 1_pInt,3_pInt)
enddo
write(6,*)' '
flush(6)
!$OMP END CRITICAL (write2out)
endif
!--------------------------------------------------------------------------------------------------
! compute the deformation gradient of individual grains due to relaxations
homID = homogenization_typeInstance(mesh_element(3,el))
F = 0.0_pReal
do iGrain = 1_pInt,homogenization_Ngrains(mesh_element(3,el))
iGrain3 = homogenization_RGC_grain1to3(iGrain,homID)
do iFace = 1_pInt,nFace
intFace = homogenization_RGC_getInterface(iFace,iGrain3) ! identifying 6 interfaces of each grain
#ifdef NEWSTATE
aVect = homogenization_RGC_relaxationVector(intFace,homID, ip, el) ! get the relaxation vectors for each interface from global relaxation vector array
#else
aVect = homogenization_RGC_relaxationVector(intFace,state,homID) ! get the relaxation vectors for each interface from global relaxation vector array
#endif
nVect = homogenization_RGC_interfaceNormal(intFace,ip,el) ! get the normal of each interface
forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt) &
F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! calculating deformation relaxations due to interface relaxation
enddo
F(1:3,1:3,iGrain) = F(1:3,1:3,iGrain) + avgF ! resulting relaxed deformation gradient
!--------------------------------------------------------------------------------------------------
! debugging the grain deformation gradients
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a32,1x,i3)')'Deformation gradient of grain: ',iGrain
do i = 1_pInt,3_pInt
write(6,'(1x,3(e15.8,1x))')(F(i,j,iGrain), j = 1_pInt,3_pInt)
enddo
write(6,*)' '
flush(6)
!$OMP END CRITICAL (write2out)
endif
enddo
end subroutine homogenization_RGC_partitionDeformation
!--------------------------------------------------------------------------------------------------
!> @brief update the internal state of the homogenization scheme and tell whether "done" and
! "happy" with result
!--------------------------------------------------------------------------------------------------
#ifdef NEWSTATE
function homogenization_RGC_updateState( P,F,F0,avgF,dt,dPdF,ip,el)
#else
function homogenization_RGC_updateState( state, state0,P,F,F0,avgF,dt,dPdF,ip,el)
#endif
use prec, only: &
p_vec
use debug, only: &
debug_level, &
debug_homogenization,&
debug_levelExtensive, &
debug_e, &
debug_i
use math, only: &
math_invert
use mesh, only: &
mesh_element
use material, only: &
homogenization_maxNgrains, &
homogenization_typeInstance, &
#ifdef NEWSTATE
homogState, &
mappingHomogenization, &
#endif
homogenization_Ngrains
use numerics, only: &
absTol_RGC, &
relTol_RGC, &
absMax_RGC, &
relMax_RGC, &
pPert_RGC, &
maxdRelax_RGC, &
viscPower_RGC, &
viscModus_RGC, &
refRelaxRate_RGC
implicit none
#ifndef NEWSTATE
type(p_vec), intent(inout) :: state !< current state
type(p_vec), intent(in) :: state0 !< initial state
#endif
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: &
P,& !< array of P
F,& !< array of F
F0 !< array of initial F
real(pReal), dimension (3,3,3,3,homogenization_maxNgrains), intent(in) :: dPdF !< array of current grain stiffness
real(pReal), dimension (3,3), intent(in) :: avgF !< average F
real(pReal), intent(in) :: dt !< time increment
integer(pInt), intent(in) :: &
ip, & !< integration point number
el !< element number
logical, dimension(2) :: homogenization_RGC_updateState
integer(pInt), dimension (4) :: intFaceN,intFaceP,faceID
integer(pInt), dimension (3) :: nGDim,iGr3N,iGr3P,stresLoc
integer(pInt), dimension (2) :: residLoc
integer(pInt) homID,iNum,i,j,nIntFaceTot,iGrN,iGrP,iMun,iFace,k,l,ipert,iGrain,nGrain
real(pReal), dimension (3,3,homogenization_maxNgrains) :: R,pF,pR,D,pD
real(pReal), dimension (3,homogenization_maxNgrains) :: NN,pNN
real(pReal), dimension (3) :: normP,normN,mornP,mornN
real(pReal) :: residMax,stresMax,constitutiveWork,penaltyEnergy,volDiscrep
logical error
integer(pInt), parameter :: nFace = 6_pInt
real(pReal), dimension(:,:), allocatable :: tract,jmatrix,jnverse,smatrix,pmatrix,rmatrix
real(pReal), dimension(:), allocatable :: resid,relax,p_relax,p_resid,drelax
if(dt < tiny(1.0_pReal)) then ! zero time step
homogenization_RGC_updateState = .true. ! pretend everything is fine and return
return
endif
!--------------------------------------------------------------------------------------------------
! get the dimension of the cluster (grains and interfaces)
homID = homogenization_typeInstance(mesh_element(3,el))
nGDim = homogenization_RGC_Ngrains(1:3,homID)
nGrain = homogenization_Ngrains(mesh_element(3,el))
nIntFaceTot = (nGDim(1)-1_pInt)*nGDim(2)*nGDim(3) + nGDim(1)*(nGDim(2)-1_pInt)*nGDim(3) &
+ nGDim(1)*nGDim(2)*(nGDim(3)-1_pInt)
!--------------------------------------------------------------------------------------------------
! allocate the size of the global relaxation arrays/jacobian matrices depending on the size of the cluster
allocate(resid(3_pInt*nIntFaceTot), source=0.0_pReal)
allocate(tract(nIntFaceTot,3), source=0.0_pReal)
#ifdef NEWSTATE
allocate(relax(3_pInt*nIntFaceTot)); relax= homogState(mappingHomogenization(2,ip,el))% &
state(1:3_pInt*nIntFaceTot,mappingHomogenization(1,ip,el))
allocate(drelax(3_pInt*nIntFaceTot)); drelax= homogState(mappingHomogenization(2,ip,el))% &
state(1:3_pInt*nIntFaceTot,mappingHomogenization(1,ip,el)) - &
homogState(mappingHomogenization(2,ip,el))% &
state0(1:3_pInt*nIntFaceTot,mappingHomogenization(1,ip,el))
#else
allocate(relax(3_pInt*nIntFaceTot)); relax=state%p(1:3_pInt*nIntFaceTot)
allocate(drelax(3_pInt*nIntFaceTot)); drelax=state%p(1:3_pInt*nIntFaceTot)-state0%p(1:3_pInt*nIntFaceTot)
#endif
!--------------------------------------------------------------------------------------------------
! debugging the obtained state
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a30)')'Obtained state: '
do i = 1_pInt,3_pInt*nIntFaceTot
#ifdef NEWSTATE
write(6,'(1x,2(e15.8,1x))')homogState(mappingHomogenization(2,ip,el))%state(i,mappingHomogenization(1,ip,el))
#else
write(6,'(1x,2(e15.8,1x))')state%p(i)
#endif
enddo
write(6,*)' '
!$OMP END CRITICAL (write2out)
endif
!--------------------------------------------------------------------------------------------------
! computing interface mismatch and stress penalty tensor for all interfaces of all grains
call homogenization_RGC_stressPenalty(R,NN,avgF,F,ip,el,homID)
!--------------------------------------------------------------------------------------------------
! calculating volume discrepancy and stress penalty related to overall volume discrepancy
call homogenization_RGC_volumePenalty(D,volDiscrep,F,avgF,ip,el)
!--------------------------------------------------------------------------------------------------
! debugging the mismatch, stress and penalties of grains
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt) then
!$OMP CRITICAL (write2out)
do iGrain = 1_pInt,nGrain
write(6,'(1x,a30,1x,i3,1x,a4,3(1x,e15.8))')'Mismatch magnitude of grain(',iGrain,') :',&
NN(1,iGrain),NN(2,iGrain),NN(3,iGrain)
write(6,'(/,1x,a30,1x,i3)')'Stress and penalties of grain: ',iGrain
do i = 1_pInt,3_pInt
write(6,'(1x,3(e15.8,1x),1x,3(e15.8,1x),1x,3(e15.8,1x))')(P(i,j,iGrain), j = 1_pInt,3_pInt), &
(R(i,j,iGrain), j = 1_pInt,3_pInt), &
(D(i,j,iGrain), j = 1_pInt,3_pInt)
enddo
write(6,*)' '
enddo
!$OMP END CRITICAL (write2out)
endif
!!!------------------------------------------------------------------------------------------------
! computing the residual stress from the balance of traction at all (interior) interfaces
do iNum = 1_pInt,nIntFaceTot
faceID = homogenization_RGC_interface1to4(iNum,homID) ! identifying the interface ID in local coordinate system (4-dimensional index)
!--------------------------------------------------------------------------------------------------
! identify the left/bottom/back grain (-|N)
iGr3N = faceID(2:4) ! identifying the grain ID in local coordinate system (3-dimensional index)
iGrN = homogenization_RGC_grain3to1(iGr3N,homID) ! translate the local grain ID into global coordinate system (1-dimensional index)
intFaceN = homogenization_RGC_getInterface(2_pInt*faceID(1),iGr3N)
normN = homogenization_RGC_interfaceNormal(intFaceN,ip,el) ! get the interface normal
!--------------------------------------------------------------------------------------------------
! identify the right/up/front grain (+|P)
iGr3P = iGr3N
iGr3P(faceID(1)) = iGr3N(faceID(1))+1_pInt ! identifying the grain ID in local coordinate system (3-dimensional index)
iGrP = homogenization_RGC_grain3to1(iGr3P,homID) ! translate the local grain ID into global coordinate system (1-dimensional index)
intFaceP = homogenization_RGC_getInterface(2_pInt*faceID(1)-1_pInt,iGr3P)
normP = homogenization_RGC_interfaceNormal(intFaceP,ip,el) ! get the interface normal
!--------------------------------------------------------------------------------------------------
! compute the residual of traction at the interface (in local system, 4-dimensional index)
do i = 1_pInt,3_pInt
tract(iNum,i) = sign(viscModus_RGC*(abs(drelax(i+3*(iNum-1_pInt)))/(refRelaxRate_RGC*dt))**viscPower_RGC, &
drelax(i+3*(iNum-1_pInt))) ! contribution from the relaxation viscosity
do j = 1_pInt,3_pInt
tract(iNum,i) = tract(iNum,i) + (P(i,j,iGrP) + R(i,j,iGrP) + D(i,j,iGrP))*normP(j) & ! contribution from material stress P, mismatch penalty R, and volume penalty D projected into the interface
+ (P(i,j,iGrN) + R(i,j,iGrN) + D(i,j,iGrN))*normN(j)
resid(i+3_pInt*(iNum-1_pInt)) = tract(iNum,i) ! translate the local residual into global 1-dimensional residual array
enddo
enddo
!--------------------------------------------------------------------------------------------------
! debugging the residual stress
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a30,1x,i3)')'Traction at interface: ',iNum
write(6,'(1x,3(e15.8,1x))')(tract(iNum,j), j = 1_pInt,3_pInt)
write(6,*)' '
!$OMP END CRITICAL (write2out)
endif
enddo
!--------------------------------------------------------------------------------------------------
! convergence check for stress residual
stresMax = maxval(abs(P)) ! get the maximum of first Piola-Kirchhoff (material) stress
stresLoc = int(maxloc(abs(P)),pInt) ! get the location of the maximum stress
residMax = maxval(abs(tract)) ! get the maximum of the residual
residLoc = int(maxloc(abs(tract)),pInt) ! get the position of the maximum residual
!--------------------------------------------------------------------------------------------------
! Debugging the convergent criteria
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt &
.and. debug_e == el .and. debug_i == ip) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a)')' '
write(6,'(1x,a,1x,i2,1x,i4)')'RGC residual check ...',ip,el
write(6,'(1x,a15,1x,e15.8,1x,a7,i3,1x,a12,i2,i2)')'Max stress: ',stresMax, &
'@ grain',stresLoc(3),'in component',stresLoc(1),stresLoc(2)
write(6,'(1x,a15,1x,e15.8,1x,a7,i3,1x,a12,i2)')'Max residual: ',residMax, &
'@ iface',residLoc(1),'in direction',residLoc(2)
flush(6)
!$OMP END CRITICAL (write2out)
endif
homogenization_RGC_updateState = .false.
!--------------------------------------------------------------------------------------------------
! If convergence reached => done and happy
if (residMax < relTol_RGC*stresMax .or. residMax < absTol_RGC) then
homogenization_RGC_updateState = .true.
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt &
.and. debug_e == el .and. debug_i == ip) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a55,/)')'... done and happy'
flush(6)
!$OMP END CRITICAL (write2out)
endif
!--------------------------------------------------------------------------------------------------
! compute/update the state for postResult, i.e., all energy densities computed by time-integration
#ifdef NEWSTATE
constitutiveWork = homogState(mappingHomogenization(2,ip,el))%state(3*nIntFaceTot+1,mappingHomogenization(1,ip,el))
penaltyEnergy = homogState(mappingHomogenization(2,ip,el))%state(3*nIntFaceTot+5,mappingHomogenization(1,ip,el))
#else
constitutiveWork = state%p(3*nIntFaceTot+1)
penaltyEnergy = state%p(3*nIntFaceTot+5)
#endif
do iGrain = 1_pInt,homogenization_Ngrains(mesh_element(3,el)) ! time-integration loop for the calculating the work and energy
do i = 1_pInt,3_pInt
do j = 1_pInt,3_pInt
constitutiveWork = constitutiveWork + P(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain))/real(nGrain,pReal)
penaltyEnergy = penaltyEnergy + R(i,j,iGrain)*(F(i,j,iGrain) - F0(i,j,iGrain))/real(nGrain,pReal)
enddo
enddo
enddo
#ifdef NEWSTATE
homogState(mappingHomogenization(2,ip,el))% &
state(3*nIntFaceTot+1,mappingHomogenization(1,ip,el)) = constitutiveWork ! the bulk mechanical/constitutive work
homogState(mappingHomogenization(2,ip,el))% &
state(3*nIntFaceTot+2,mappingHomogenization(1,ip,el)) = sum(NN(1,:))/real(nGrain,pReal) ! the overall mismatch of all interface normal to e1-direction
homogState(mappingHomogenization(2,ip,el))% &
state(3*nIntFaceTot+3,mappingHomogenization(1,ip,el)) = sum(NN(2,:))/real(nGrain,pReal) ! the overall mismatch of all interface normal to e2-direction
homogState(mappingHomogenization(2,ip,el))% &
state(3*nIntFaceTot+4,mappingHomogenization(1,ip,el)) = sum(NN(3,:))/real(nGrain,pReal) ! the overall mismatch of all interface normal to e3-direction
homogState(mappingHomogenization(2,ip,el))% &
state(3*nIntFaceTot+5,mappingHomogenization(1,ip,el)) = penaltyEnergy ! the overall penalty energy
homogState(mappingHomogenization(2,ip,el))% &
state(3*nIntFaceTot+6,mappingHomogenization(1,ip,el)) = volDiscrep ! the overall volume discrepancy
homogState(mappingHomogenization(2,ip,el))% &
state(3*nIntFaceTot+7,mappingHomogenization(1,ip,el)) = &
sum(abs(drelax))/dt/real(3_pInt*nIntFaceTot,pReal) ! the average rate of relaxation vectors
homogState(mappingHomogenization(2,ip,el))% &
state(3*nIntFaceTot+8,mappingHomogenization(1,ip,el)) = maxval(abs(drelax))/dt ! the maximum rate of relaxation vectors
#else
state%p(3*nIntFaceTot+1) = constitutiveWork ! the bulk mechanical/constitutive work
state%p(3*nIntFaceTot+2) = sum(NN(1,:))/real(nGrain,pReal) ! the overall mismatch of all interface normal to e1-direction
state%p(3*nIntFaceTot+3) = sum(NN(2,:))/real(nGrain,pReal) ! the overall mismatch of all interface normal to e2-direction
state%p(3*nIntFaceTot+4) = sum(NN(3,:))/real(nGrain,pReal) ! the overall mismatch of all interface normal to e3-direction
state%p(3*nIntFaceTot+5) = penaltyEnergy ! the overall penalty energy
state%p(3*nIntFaceTot+6) = volDiscrep ! the overall volume discrepancy
state%p(3*nIntFaceTot+7) = sum(abs(drelax))/dt/real(3_pInt*nIntFaceTot,pReal) ! the average rate of relaxation vectors
state%p(3*nIntFaceTot+8) = maxval(abs(drelax))/dt ! the maximum rate of relaxation vectors
#endif
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt &
.and. debug_e == el .and. debug_i == ip) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a30,1x,e15.8)') 'Constitutive work: ',constitutiveWork
write(6,'(1x,a30,3(1x,e15.8))')'Magnitude mismatch: ',sum(NN(1,:))/real(nGrain,pReal), &
sum(NN(2,:))/real(nGrain,pReal), &
sum(NN(3,:))/real(nGrain,pReal)
write(6,'(1x,a30,1x,e15.8)') 'Penalty energy: ',penaltyEnergy
write(6,'(1x,a30,1x,e15.8,/)') 'Volume discrepancy: ',volDiscrep
write(6,'(1x,a30,1x,e15.8)') 'Maximum relaxation rate: ',maxval(abs(drelax))/dt
write(6,'(1x,a30,1x,e15.8,/)') 'Average relaxation rate: ',sum(abs(drelax))/dt/real(3_pInt*nIntFaceTot,pReal)
flush(6)
!$OMP END CRITICAL (write2out)
endif
deallocate(tract,resid,relax,drelax)
return
!--------------------------------------------------------------------------------------------------
! if residual blows-up => done but unhappy
elseif (residMax > relMax_RGC*stresMax .or. residMax > absMax_RGC) then ! try to restart when residual blows up exceeding maximum bound
homogenization_RGC_updateState = [.true.,.false.] ! with direct cut-back
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt &
.and. debug_e == el .and. debug_i == ip) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a55,/)')'... broken'
flush(6)
!$OMP END CRITICAL (write2out)
endif
deallocate(tract,resid,relax,drelax)
return
else ! proceed with computing the Jacobian and state update
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt &
.and. debug_e == el .and. debug_i == ip) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a55,/)')'... not yet done'
flush(6)
!$OMP END CRITICAL (write2out)
endif
endif
!!!------------------------------------------------------------------------------------------------
! construct the global Jacobian matrix for updating the global relaxation vector array when
! convergence is not yet reached ...
!--------------------------------------------------------------------------------------------------
! ... of the constitutive stress tangent, assembled from dPdF or material constitutive model "smatrix"
allocate(smatrix(3*nIntFaceTot,3*nIntFaceTot), source=0.0_pReal)
do iNum = 1_pInt,nIntFaceTot
faceID = homogenization_RGC_interface1to4(iNum,homID) ! assembling of local dPdF into global Jacobian matrix
!--------------------------------------------------------------------------------------------------
! identify the left/bottom/back grain (-|N)
iGr3N = faceID(2:4) ! identifying the grain ID in local coordinate sytem
iGrN = homogenization_RGC_grain3to1(iGr3N,homID) ! translate into global grain ID
intFaceN = homogenization_RGC_getInterface(2_pInt*faceID(1),iGr3N) ! identifying the connecting interface in local coordinate system
normN = homogenization_RGC_interfaceNormal(intFaceN,ip,el) ! get the interface normal
do iFace = 1_pInt,nFace
intFaceN = homogenization_RGC_getInterface(iFace,iGr3N) ! identifying all interfaces that influence relaxation of the above interface
mornN = homogenization_RGC_interfaceNormal(intFaceN,ip,el) ! get normal of the interfaces
iMun = homogenization_RGC_interface4to1(intFaceN,homID) ! translate the interfaces ID into local 4-dimensional index
if (iMun > 0) then ! get the corresponding tangent
do i=1_pInt,3_pInt; do j=1_pInt,3_pInt; do k=1_pInt,3_pInt; do l=1_pInt,3_pInt
smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) + dPdF(i,k,j,l,iGrN)*normN(k)*mornN(l)
enddo;enddo;enddo;enddo
! projecting the material tangent dPdF into the interface
! to obtain the Jacobian matrix contribution of dPdF
endif
enddo
!--------------------------------------------------------------------------------------------------
! identify the right/up/front grain (+|P)
iGr3P = iGr3N
iGr3P(faceID(1)) = iGr3N(faceID(1))+1_pInt ! identifying the grain ID in local coordinate sytem
iGrP = homogenization_RGC_grain3to1(iGr3P,homID) ! translate into global grain ID
intFaceP = homogenization_RGC_getInterface(2_pInt*faceID(1)-1_pInt,iGr3P) ! identifying the connecting interface in local coordinate system
normP = homogenization_RGC_interfaceNormal(intFaceP,ip,el) ! get the interface normal
do iFace = 1_pInt,nFace
intFaceP = homogenization_RGC_getInterface(iFace,iGr3P) ! identifying all interfaces that influence relaxation of the above interface
mornP = homogenization_RGC_interfaceNormal(intFaceP,ip,el) ! get normal of the interfaces
iMun = homogenization_RGC_interface4to1(intFaceP,homID) ! translate the interfaces ID into local 4-dimensional index
if (iMun > 0_pInt) then ! get the corresponding tangent
do i=1_pInt,3_pInt; do j=1_pInt,3_pInt; do k=1_pInt,3_pInt; do l=1_pInt,3_pInt
smatrix(3*(iNum-1)+i,3*(iMun-1)+j) = smatrix(3*(iNum-1)+i,3*(iMun-1)+j) + dPdF(i,k,j,l,iGrP)*normP(k)*mornP(l)
enddo;enddo;enddo;enddo
endif
enddo
enddo
!--------------------------------------------------------------------------------------------------
! debugging the global Jacobian matrix of stress tangent
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a30)')'Jacobian matrix of stress'
do i = 1_pInt,3_pInt*nIntFaceTot
write(6,'(1x,100(e11.4,1x))')(smatrix(i,j), j = 1_pInt,3_pInt*nIntFaceTot)
enddo
write(6,*)' '
flush(6)
!$OMP END CRITICAL (write2out)
endif
!--------------------------------------------------------------------------------------------------
! ... of the stress penalty tangent (mismatch penalty and volume penalty, computed using numerical
! perturbation method) "pmatrix"
allocate(pmatrix(3*nIntFaceTot,3*nIntFaceTot), source=0.0_pReal)
allocate(p_relax(3*nIntFaceTot), source=0.0_pReal)
allocate(p_resid(3*nIntFaceTot), source=0.0_pReal)
do ipert = 1_pInt,3_pInt*nIntFaceTot
p_relax = relax
p_relax(ipert) = relax(ipert) + pPert_RGC ! perturb the relaxation vector
#ifdef NEWSTATE
homogState(mappingHomogenization(2,ip,el))%state(1:3*nIntFaceTot,mappingHomogenization(1,ip,el)) = p_relax
call homogenization_RGC_grainDeformation(pF,avgF,ip,el) ! compute the grains deformation from perturbed state
#else
state%p(1:3*nIntFaceTot) = p_relax
call homogenization_RGC_grainDeformation(pF,avgF,state,ip,el) ! compute the grains deformation from perturbed state
#endif
call homogenization_RGC_stressPenalty(pR,pNN,avgF,pF,ip,el,homID) ! compute stress penalty due to interface mismatch from perturbed state
call homogenization_RGC_volumePenalty(pD,volDiscrep,pF,avgF,ip,el) ! compute stress penalty due to volume discrepancy from perturbed state
!--------------------------------------------------------------------------------------------------
! computing the global stress residual array from the perturbed state
p_resid = 0.0_pReal
do iNum = 1_pInt,nIntFaceTot
faceID = homogenization_RGC_interface1to4(iNum,homID) ! identifying the interface ID in local coordinate system (4-dimensional index)
!--------------------------------------------------------------------------------------------------
! identify the left/bottom/back grain (-|N)
iGr3N = faceID(2:4) ! identifying the grain ID in local coordinate system (3-dimensional index)
iGrN = homogenization_RGC_grain3to1(iGr3N,homID) ! translate the local grain ID into global coordinate system (1-dimensional index)
intFaceN = homogenization_RGC_getInterface(2_pInt*faceID(1),iGr3N) ! identifying the interface ID of the grain
normN = homogenization_RGC_interfaceNormal(intFaceN,ip,el) ! get the corresponding interface normal
!--------------------------------------------------------------------------------------------------
! identify the right/up/front grain (+|P)
iGr3P = iGr3N
iGr3P(faceID(1)) = iGr3N(faceID(1))+1_pInt ! identifying the grain ID in local coordinate system (3-dimensional index)
iGrP = homogenization_RGC_grain3to1(iGr3P,homID) ! translate the local grain ID into global coordinate system (1-dimensional index)
intFaceP = homogenization_RGC_getInterface(2_pInt*faceID(1)-1_pInt,iGr3P) ! identifying the interface ID of the grain
normP = homogenization_RGC_interfaceNormal(intFaceP,ip,el) ! get the corresponding normal
!--------------------------------------------------------------------------------------------------
! compute the residual stress (contribution of mismatch and volume penalties) from perturbed state
! at all interfaces
do i = 1_pInt,3_pInt; do j = 1_pInt,3_pInt
p_resid(i+3*(iNum-1)) = p_resid(i+3*(iNum-1)) + (pR(i,j,iGrP) - R(i,j,iGrP))*normP(j) &
+ (pR(i,j,iGrN) - R(i,j,iGrN))*normN(j) &
+ (pD(i,j,iGrP) - D(i,j,iGrP))*normP(j) &
+ (pD(i,j,iGrN) - D(i,j,iGrN))*normN(j)
enddo; enddo
enddo
pmatrix(:,ipert) = p_resid/pPert_RGC
enddo
!--------------------------------------------------------------------------------------------------
! debugging the global Jacobian matrix of penalty tangent
if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a30)')'Jacobian matrix of penalty'
do i = 1_pInt,3_pInt*nIntFaceTot
write(6,'(1x,100(e11.4,1x))')(pmatrix(i,j), j = 1_pInt,3_pInt*nIntFaceTot)
enddo
write(6,*)' '
flush(6)
!$OMP END CRITICAL (write2out)
endif
!--------------------------------------------------------------------------------------------------
! ... of the numerical viscosity traction "rmatrix"
allocate(rmatrix(3*nIntFaceTot,3*nIntFaceTot),source=0.0_pReal)
forall (i=1_pInt:3_pInt*nIntFaceTot) &
rmatrix(i,i) = viscModus_RGC*viscPower_RGC/(refRelaxRate_RGC*dt)* & ! tangent due to numerical viscosity traction appears
(abs(drelax(i))/(refRelaxRate_RGC*dt))**(viscPower_RGC - 1.0_pReal) ! only in the main diagonal term
!--------------------------------------------------------------------------------------------------
! debugging the global Jacobian matrix of numerical viscosity tangent
if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a30)')'Jacobian matrix of penalty'
do i = 1_pInt,3_pInt*nIntFaceTot
write(6,'(1x,100(e11.4,1x))')(rmatrix(i,j), j = 1_pInt,3_pInt*nIntFaceTot)
enddo
write(6,*)' '
flush(6)
!$OMP END CRITICAL (write2out)
endif
!--------------------------------------------------------------------------------------------------
! The overall Jacobian matrix summarizing contributions of smatrix, pmatrix, rmatrix
allocate(jmatrix(3*nIntFaceTot,3*nIntFaceTot)); jmatrix = smatrix + pmatrix + rmatrix
if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a30)')'Jacobian matrix (total)'
do i = 1_pInt,3_pInt*nIntFaceTot
write(6,'(1x,100(e11.4,1x))')(jmatrix(i,j), j = 1_pInt,3_pInt*nIntFaceTot)
enddo
write(6,*)' '
flush(6)
!$OMP END CRITICAL (write2out)
endif
!--------------------------------------------------------------------------------------------------
! computing the update of the state variable (relaxation vectors) using the Jacobian matrix
allocate(jnverse(3_pInt*nIntFaceTot,3_pInt*nIntFaceTot),source=0.0_pReal)
call math_invert(size(jmatrix,1),jmatrix,jnverse,error) ! Compute the inverse of the overall Jacobian matrix
!--------------------------------------------------------------------------------------------------
! debugging the inverse Jacobian matrix
if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a30)')'Jacobian inverse'
do i = 1_pInt,3_pInt*nIntFaceTot
write(6,'(1x,100(e11.4,1x))')(jnverse(i,j), j = 1_pInt,3_pInt*nIntFaceTot)
enddo
write(6,*)' '
flush(6)
!$OMP END CRITICAL (write2out)
endif
!--------------------------------------------------------------------------------------------------
! calculate the state update (global relaxation vectors) for the next Newton-Raphson iteration
drelax = 0.0_pReal
do i = 1_pInt,3_pInt*nIntFaceTot
do j = 1_pInt,3_pInt*nIntFaceTot
drelax(i) = drelax(i) - jnverse(i,j)*resid(j) ! Calculate the correction for the state variable
enddo
enddo
relax = relax + drelax ! Updateing the state variable for the next iteration
#ifdef NEWSTATE
homogState(mappingHomogenization(2,ip,el))%state(1:3*nIntFaceTot,mappingHomogenization(1,ip,el)) = relax
#else
state%p(1:3*nIntFaceTot) = relax
#endif
if (any(abs(drelax) > maxdRelax_RGC)) then ! Forcing cutback when the incremental change of relaxation vector becomes too large
homogenization_RGC_updateState = [.true.,.false.]
!$OMP CRITICAL (write2out)
write(6,'(1x,a,1x,i3,1x,a,1x,i3,1x,a)')'RGC_updateState: ip',ip,'| el',el,'enforces cutback'
write(6,'(1x,a,1x,e15.8)')'due to large relaxation change =',maxval(abs(drelax))
flush(6)
!$OMP END CRITICAL (write2out)
endif
!--------------------------------------------------------------------------------------------------
! debugging the return state
if (iand(debug_homogenization, debug_levelExtensive) > 0_pInt) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a30)')'Returned state: '
do i = 1_pInt,3_pInt*nIntFaceTot
#ifdef NEWSTATE
write(6,'(1x,2(e15.8,1x))')homogState(mappingHomogenization(2,ip,el))%state(i,mappingHomogenization(1,ip,el))
#else
write(6,'(1x,2(e15.8,1x))')state%p(i)
#endif
enddo
write(6,*)' '
flush(6)
!$OMP END CRITICAL (write2out)
endif
deallocate(tract,resid,jmatrix,jnverse,relax,drelax,pmatrix,smatrix,p_relax,p_resid)
end function homogenization_RGC_updateState
!--------------------------------------------------------------------------------------------------
!> @brief derive average stress and stiffness from constituent quantities
!--------------------------------------------------------------------------------------------------
subroutine homogenization_RGC_averageStressAndItsTangent(avgP,dAvgPdAvgF,P,dPdF,el)
use debug, only: &
debug_level, &
debug_homogenization,&
debug_levelExtensive
use mesh, only: mesh_element
use material, only: &
homogenization_maxNgrains, &
#ifdef NEWSTATE
homogState, &
mappingHomogenization, &
#endif
homogenization_Ngrains, &
homogenization_typeInstance
use math, only: math_Plain3333to99
implicit none
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 (3,3,homogenization_maxNgrains), intent(in) :: P !< array of current grain stresses
real(pReal), dimension (3,3,3,3,homogenization_maxNgrains), intent(in) :: dPdF !< array of current grain stiffnesses
integer(pInt), intent(in) :: el !< element number
real(pReal), dimension (9,9) :: dPdF99
integer(pInt) :: homID, i, j, Ngrains, iGrain
homID = homogenization_typeInstance(mesh_element(3,el))
Ngrains = homogenization_Ngrains(mesh_element(3,el))
!--------------------------------------------------------------------------------------------------
! debugging the grain tangent
if (iand(debug_level(debug_homogenization), debug_levelExtensive) /= 0_pInt) then
!$OMP CRITICAL (write2out)
do iGrain = 1_pInt,Ngrains
dPdF99 = math_Plain3333to99(dPdF(1:3,1:3,1:3,1:3,iGrain))
write(6,'(1x,a30,1x,i3)')'Stress tangent of grain: ',iGrain
do i = 1_pInt,9_pInt
write(6,'(1x,(e15.8,1x))') (dPdF99(i,j), j = 1_pInt,9_pInt)
enddo
write(6,*)' '
enddo
flush(6)
!$OMP END CRITICAL (write2out)
endif
!--------------------------------------------------------------------------------------------------
! computing the average first Piola-Kirchhoff stress P and the average tangent dPdF
avgP = sum(P,3)/real(Ngrains,pReal)
dAvgPdAvgF = sum(dPdF,5)/real(Ngrains,pReal)
end subroutine homogenization_RGC_averageStressAndItsTangent
!--------------------------------------------------------------------------------------------------
!> @brief return array of homogenization results for post file inclusion
!--------------------------------------------------------------------------------------------------
#ifdef NEWSTATE
pure function homogenization_RGC_postResults(ip,el,avgP,avgF)
#else
pure function homogenization_RGC_postResults(state,ip,el,avgP,avgF)
#endif
use prec, only: &
p_vec
use mesh, only: &
mesh_element, &
mesh_ipCoordinates
use material, only: &
homogenization_typeInstance,&
#ifdef NEWSTATE
homogState, &
mappingHomogenization, &
#endif
homogenization_Noutput
use crystallite, only: &
crystallite_temperature
implicit none
integer(pInt), intent(in) :: &
ip, & !< integration point number
el !< element number
real(pReal), dimension(3,3), intent(in) :: &
avgP, & !< average stress at material point
avgF !< average deformation gradient at material point
#ifndef NEWSTATE
type(p_vec), intent(in) :: &
state ! my State
#endif
integer(pInt) homID,o,c,nIntFaceTot
real(pReal), dimension(homogenization_RGC_sizePostResults(homogenization_typeInstance(mesh_element(3,el)))) :: &
homogenization_RGC_postResults
homID = homogenization_typeInstance(mesh_element(3,el))
nIntFaceTot=(homogenization_RGC_Ngrains(1,homID)-1_pInt)*homogenization_RGC_Ngrains(2,homID)*homogenization_RGC_Ngrains(3,homID)&
+ homogenization_RGC_Ngrains(1,homID)*(homogenization_RGC_Ngrains(2,homID)-1_pInt)*homogenization_RGC_Ngrains(3,homID)&
+ homogenization_RGC_Ngrains(1,homID)*homogenization_RGC_Ngrains(2,homID)*(homogenization_RGC_Ngrains(3,homID)-1_pInt)
c = 0_pInt
homogenization_RGC_postResults = 0.0_pReal
do o = 1_pInt,homogenization_Noutput(mesh_element(3,el))
select case(homogenization_RGC_outputID(o,homID))
case (temperature_ID)
homogenization_RGC_postResults(c+1_pInt) = crystallite_temperature(ip,el)
c = c + 1_pInt
case (avgdefgrad_ID)
homogenization_RGC_postResults(c+1_pInt:c+9_pInt) = reshape(avgF,[9])
c = c + 9_pInt
case (avgfirstpiola_ID)
homogenization_RGC_postResults(c+1_pInt:c+9_pInt) = reshape(avgP,[9])
c = c + 9_pInt
case (ipcoords_ID)
homogenization_RGC_postResults(c+1_pInt:c+3_pInt) = mesh_ipCoordinates(1:3,ip,el) ! current ip coordinates
c = c + 3_pInt
case (constitutivework_ID)
#ifdef NEWSTATE
homogenization_RGC_postResults(c+1) = homogState(mappingHomogenization(2,ip,el))% &
state(3*nIntFaceTot+1,mappingHomogenization(1,ip,el))
#else
homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+1)
#endif
c = c + 1_pInt
case (magnitudemismatch_ID)
#ifdef NEWSTATE
homogenization_RGC_postResults(c+1) = homogState(mappingHomogenization(2,ip,el))% &
state(3*nIntFaceTot+2,mappingHomogenization(1,ip,el))
homogenization_RGC_postResults(c+2) = homogState(mappingHomogenization(2,ip,el))% &
state(3*nIntFaceTot+3,mappingHomogenization(1,ip,el))
homogenization_RGC_postResults(c+3) = homogState(mappingHomogenization(2,ip,el))% &
state(3*nIntFaceTot+4,mappingHomogenization(1,ip,el))
#else
homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+2)
homogenization_RGC_postResults(c+2) = state%p(3*nIntFaceTot+3)
homogenization_RGC_postResults(c+3) = state%p(3*nIntFaceTot+4)
#endif
c = c + 3_pInt
case (penaltyenergy_ID)
#ifdef NEWSTATE
homogenization_RGC_postResults(c+1) = homogState(mappingHomogenization(2,ip,el))% &
state(3*nIntFaceTot+5,mappingHomogenization(1,ip,el))
c = c + 1_pInt
#else
homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+5)
c = c + 1_pInt
#endif
case (volumediscrepancy_ID)
#ifdef NEWSTATE
homogenization_RGC_postResults(c+1) = homogState(mappingHomogenization(2,ip,el))% &
state(3*nIntFaceTot+6,mappingHomogenization(1,ip,el))
c = c + 1_pInt
#else
homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+6)
c = c + 1_pInt
#endif
case (averagerelaxrate_ID)
#ifdef NEWSTATE
homogenization_RGC_postResults(c+1) = homogState(mappingHomogenization(2,ip,el))% &
state(3*nIntFaceTot+7,mappingHomogenization(1,ip,el))
c = c + 1_pInt
#else
homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+7)
c = c + 1_pInt
#endif
case (maximumrelaxrate_ID)
#ifdef NEWSTATE
homogenization_RGC_postResults(c+1) = homogState(mappingHomogenization(2,ip,el))% &
state(3*nIntFaceTot+8,mappingHomogenization(1,ip,el))
c = c + 1_pInt
#else
homogenization_RGC_postResults(c+1) = state%p(3*nIntFaceTot+8)
c = c + 1_pInt
#endif
end select
enddo
end function homogenization_RGC_postResults
!--------------------------------------------------------------------------------------------------
!> @brief calculate stress-like penalty due to deformation mismatch
!--------------------------------------------------------------------------------------------------
subroutine homogenization_RGC_stressPenalty(rPen,nMis,avgF,fDef,ip,el,homID)
use debug, only: &
debug_level, &
debug_homogenization,&
debug_levelExtensive, &
debug_e, &
debug_i
use mesh, only: &
mesh_element
use constitutive, only: &
constitutive_homogenizedC
use math, only: &
math_civita,&
math_invert33
use material, only: &
homogenization_maxNgrains,&
#ifdef NEWSTATE
homogState, &
mappingHomogenization, &
#endif
homogenization_Ngrains
use numerics, only: &
xSmoo_RGC
implicit none
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: rPen !< stress-like penalty
real(pReal), dimension (3,homogenization_maxNgrains), intent(out) :: nMis !< total amount of mismatch
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: fDef !< deformation gradients
real(pReal), dimension (3,3), intent(in) :: avgF !< initial effective stretch tensor
integer(pInt), intent(in) :: ip,el
integer(pInt), dimension (4) :: intFace
integer(pInt), dimension (3) :: iGrain3,iGNghb3,nGDim
real(pReal), dimension (3,3) :: gDef,nDef
real(pReal), dimension (3) :: nVect,surfCorr
real(pReal), dimension (2) :: Gmoduli
integer(pInt) :: homID,iGrain,iGNghb,iFace,i,j,k,l
real(pReal) :: muGrain,muGNghb,nDefNorm,bgGrain,bgGNghb
integer(pInt), parameter :: nFace = 6_pInt
real(pReal), parameter :: nDefToler = 1.0e-10_pReal
nGDim = homogenization_RGC_Ngrains(1:3,homID)
rPen = 0.0_pReal
nMis = 0.0_pReal
!--------------------------------------------------------------------------------------------------
! get the correction factor the modulus of penalty stress representing the evolution of area of
! the interfaces due to deformations
surfCorr = homogenization_RGC_surfaceCorrection(avgF,ip,el)
!--------------------------------------------------------------------------------------------------
! debugging the surface correction factor
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt &
.and. debug_e == el .and. debug_i == ip) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a20,2(1x,i3))')'Correction factor: ',ip,el
write(6,'(1x,3(e11.4,1x))')(surfCorr(i), i = 1,3)
!$OMP END CRITICAL (write2out)
endif
!--------------------------------------------------------------------------------------------------
! computing the mismatch and penalty stress tensor of all grains
do iGrain = 1_pInt,homogenization_Ngrains(mesh_element(3,el))
Gmoduli = homogenization_RGC_equivalentModuli(iGrain,ip,el)
muGrain = Gmoduli(1) ! collecting the equivalent shear modulus of grain
bgGrain = Gmoduli(2) ! and the lengthh of Burgers vector
iGrain3 = homogenization_RGC_grain1to3(iGrain,homID) ! get the grain ID in local 3-dimensional index (x,y,z)-position
!* Looping over all six interfaces of each grain
do iFace = 1_pInt,nFace
intFace = homogenization_RGC_getInterface(iFace,iGrain3) ! get the 4-dimensional index of the interface in local numbering system of the grain
nVect = homogenization_RGC_interfaceNormal(intFace,ip,el) ! get the interface normal
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),pInt)
if (iGNghb3(1) < 1) iGNghb3(1) = nGDim(1) ! with periodicity along e1 direction
if (iGNghb3(1) > nGDim(1)) iGNghb3(1) = 1_pInt
if (iGNghb3(2) < 1) iGNghb3(2) = nGDim(2) ! with periodicity along e2 direction
if (iGNghb3(2) > nGDim(2)) iGNghb3(2) = 1_pInt
if (iGNghb3(3) < 1) iGNghb3(3) = nGDim(3) ! with periodicity along e3 direction
if (iGNghb3(3) > nGDim(3)) iGNghb3(3) = 1_pInt
iGNghb = homogenization_RGC_grain3to1(iGNghb3,homID) ! get the ID of the neighboring grain
Gmoduli = homogenization_RGC_equivalentModuli(iGNghb,ip,el) ! collecting the shear modulus and Burgers vector of the neighbor
muGNghb = Gmoduli(1)
bgGNghb = Gmoduli(2)
gDef = 0.5_pReal*(fDef(1:3,1:3,iGNghb) - fDef(1:3,1:3,iGrain)) ! compute the difference/jump in deformation gradeint across the neighbor
!--------------------------------------------------------------------------------------------------
! compute the mismatch tensor of all interfaces
nDefNorm = 0.0_pReal
nDef = 0.0_pReal
do i = 1_pInt,3_pInt; do j = 1_pInt,3_pInt
do k = 1_pInt,3_pInt; do l = 1_pInt,3_pInt
nDef(i,j) = nDef(i,j) - nVect(k)*gDef(i,l)*math_civita(j,k,l) ! compute the interface mismatch tensor from the jump of deformation gradient
enddo; enddo
nDefNorm = nDefNorm + nDef(i,j)*nDef(i,j) ! compute the norm of the mismatch tensor
enddo; enddo
nDefNorm = max(nDefToler,sqrt(nDefNorm)) ! approximation to zero mismatch if mismatch is zero (singularity)
nMis(abs(intFace(1)),iGrain) = nMis(abs(intFace(1)),iGrain) + nDefNorm ! total amount of mismatch experienced by the grain (at all six interfaces)
!--------------------------------------------------------------------------------------------------
! debuggin the mismatch tensor
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt &
.and. debug_e == el .and. debug_i == ip) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a20,i2,1x,a20,1x,i3)')'Mismatch to face: ',intFace(1),'neighbor grain: ',iGNghb
do i = 1,3
write(6,'(1x,3(e11.4,1x))')(nDef(i,j), j = 1,3)
enddo
write(6,'(1x,a20,e11.4)')'with magnitude: ',nDefNorm
!$OMP END CRITICAL (write2out)
endif
!--------------------------------------------------------------------------------------------------
! compute the stress penalty of all interfaces
do i = 1_pInt,3_pInt; do j = 1_pInt,3_pInt
do k = 1_pInt,3_pInt; do l = 1_pInt,3_pInt
rPen(i,j,iGrain) = rPen(i,j,iGrain) + 0.5_pReal*(muGrain*bgGrain + muGNghb*bgGNghb)*homogenization_RGC_xiAlpha(homID) &
*surfCorr(abs(intFace(1)))/homogenization_RGC_dAlpha(abs(intFace(1)),homID) &
*cosh(homogenization_RGC_ciAlpha(homID)*nDefNorm) &
*0.5_pReal*nVect(l)*nDef(i,k)/nDefNorm*math_civita(k,l,j) &
*tanh(nDefNorm/xSmoo_RGC)
enddo; enddo
enddo; enddo
enddo
!--------------------------------------------------------------------------------------------------
! debugging the stress-like penalty
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt &
.and. debug_e == el .and. debug_i == ip) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a20,i2)')'Penalty of grain: ',iGrain
do i = 1,3
write(6,'(1x,3(e11.4,1x))')(rPen(i,j,iGrain), j = 1,3)
enddo
!$OMP END CRITICAL (write2out)
endif
enddo
end subroutine homogenization_RGC_stressPenalty
!--------------------------------------------------------------------------------------------------
!> @brief calculate stress-like penalty due to volume discrepancy
!--------------------------------------------------------------------------------------------------
subroutine homogenization_RGC_volumePenalty(vPen,vDiscrep,fDef,fAvg,ip,el)
use debug, only: &
debug_level, &
debug_homogenization,&
debug_levelExtensive, &
debug_e, &
debug_i
use mesh, only: &
mesh_element
use math, only: &
math_det33, &
math_inv33
use material, only: &
homogenization_maxNgrains,&
#ifdef NEWSTATE
homogState, &
mappingHomogenization, &
#endif
homogenization_Ngrains
use numerics, only: &
maxVolDiscr_RGC,&
volDiscrMod_RGC,&
volDiscrPow_RGC
implicit none
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: vPen ! stress-like penalty due to volume
real(pReal), intent(out) :: vDiscrep ! total volume discrepancy
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(in) :: fDef ! deformation gradients
real(pReal), dimension (3,3), intent(in) :: fAvg ! overall deformation gradient
integer(pInt), intent(in) :: ip,& ! integration point
el
real(pReal), dimension (homogenization_maxNgrains) :: gVol
integer(pInt) :: iGrain,nGrain,i,j
nGrain = homogenization_Ngrains(mesh_element(3,el))
!--------------------------------------------------------------------------------------------------
! compute the volumes of grains and of cluster
vDiscrep = math_det33(fAvg) ! compute the volume of the cluster
do iGrain = 1_pInt,nGrain
gVol(iGrain) = math_det33(fDef(1:3,1:3,iGrain)) ! compute the volume of individual grains
vDiscrep = vDiscrep - gVol(iGrain)/real(nGrain,pReal) ! calculate the difference/dicrepancy between
! the volume of the cluster and the the total volume of grains
enddo
!--------------------------------------------------------------------------------------------------
! calculate the stress and penalty due to volume discrepancy
vPen = 0.0_pReal
do iGrain = 1_pInt,nGrain
vPen(:,:,iGrain) = -1.0_pReal/real(nGrain,pReal)*volDiscrMod_RGC*volDiscrPow_RGC/maxVolDiscr_RGC* &
sign((abs(vDiscrep)/maxVolDiscr_RGC)**(volDiscrPow_RGC - 1.0),vDiscrep)* &
gVol(iGrain)*transpose(math_inv33(fDef(:,:,iGrain)))
!--------------------------------------------------------------------------------------------------
! debugging the stress-like penalty
if (iand(debug_level(debug_homogenization),debug_levelExtensive) /= 0_pInt &
.and. debug_e == el .and. debug_i == ip) then
!$OMP CRITICAL (write2out)
write(6,'(1x,a30,i2)')'Volume penalty of grain: ',iGrain
do i = 1,3
write(6,'(1x,3(e11.4,1x))')(vPen(i,j,iGrain), j = 1,3)
enddo
!$OMP END CRITICAL (write2out)
endif
enddo
end subroutine homogenization_RGC_volumePenalty
!--------------------------------------------------------------------------------------------------
!> @brief compute the correction factor accouted for surface evolution (area change) due to
! deformation
!--------------------------------------------------------------------------------------------------
function homogenization_RGC_surfaceCorrection(avgF,ip,el)
use math, only: &
math_invert33, &
math_mul33x33
#ifdef NEWSTATE
use material, only: &
homogState, &
mappingHomogenization
#endif
implicit none
real(pReal), dimension(3) :: homogenization_RGC_surfaceCorrection
real(pReal), dimension(3,3), intent(in) :: avgF !< average F
integer(pInt), intent(in) :: ip,& !< integration point number
el !< element number
real(pReal), dimension(3,3) :: invC,avgC
real(pReal), dimension(3) :: nVect
real(pReal) :: detF
integer(pInt), dimension(4) :: intFace
integer(pInt) :: i,j,iBase
logical :: error
avgC = math_mul33x33(transpose(avgF),avgF)
call math_invert33(avgC,invC,detF,error)
homogenization_RGC_surfaceCorrection = 0.0_pReal
do iBase = 1_pInt,3_pInt
intFace = [iBase,1_pInt,1_pInt,1_pInt]
nVect = homogenization_RGC_interfaceNormal(intFace,ip,el) ! get the normal of the interface
do i = 1_pInt,3_pInt; do j = 1_pInt,3_pInt
homogenization_RGC_surfaceCorrection(iBase) = & ! compute the component of (the inverse of) the stretch in the direction of the normal
homogenization_RGC_surfaceCorrection(iBase) + invC(i,j)*nVect(i)*nVect(j)
enddo; enddo
homogenization_RGC_surfaceCorrection(iBase) = & ! get the surface correction factor (area contraction/enlargement)
sqrt(homogenization_RGC_surfaceCorrection(iBase))*detF
enddo
end function homogenization_RGC_surfaceCorrection
!--------------------------------------------------------------------------------------------------
!> @brief compute the equivalent shear and bulk moduli from the elasticity tensor
!--------------------------------------------------------------------------------------------------
function homogenization_RGC_equivalentModuli(grainID,ip,el)
use constitutive, only: &
constitutive_homogenizedC
implicit none
integer(pInt), intent(in) :: &
grainID,&
ip, & !< integration point number
el !< element number
real(pReal), dimension (6,6) :: elasTens
real(pReal), dimension(2) :: homogenization_RGC_equivalentModuli
real(pReal) :: &
cEquiv_11, &
cEquiv_12, &
cEquiv_44
elasTens = constitutive_homogenizedC(grainID,ip,el)
!--------------------------------------------------------------------------------------------------
! compute the equivalent shear modulus after Turterltaub and Suiker, JMPS (2005)
cEquiv_11 = (elasTens(1,1) + elasTens(2,2) + elasTens(3,3))/3.0_pReal
cEquiv_12 = (elasTens(1,2) + elasTens(2,3) + elasTens(3,1) + &
elasTens(1,3) + elasTens(2,1) + elasTens(3,2))/6.0_pReal
cEquiv_44 = (elasTens(4,4) + elasTens(5,5) + elasTens(6,6))/3.0_pReal
homogenization_RGC_equivalentModuli(1) = 0.2_pReal*(cEquiv_11 - cEquiv_12) + 0.6_pReal*cEquiv_44
!--------------------------------------------------------------------------------------------------
! obtain the length of Burgers vector (could be model dependend)
homogenization_RGC_equivalentModuli(2) = 2.5e-10_pReal
end function homogenization_RGC_equivalentModuli
!--------------------------------------------------------------------------------------------------
!> @brief collect relaxation vectors of an interface
!--------------------------------------------------------------------------------------------------
#ifdef NEWSTATE
function homogenization_RGC_relaxationVector(intFace,homID, ip, el)
#else
function homogenization_RGC_relaxationVector(intFace,state,homID)
#endif
use prec, only: &
p_vec
#ifdef NEWSTATE
use material, only: &
homogState, &
mappingHomogenization
#endif
implicit none
#ifdef NEWSTATE
integer(pInt), intent(in) :: ip, el
#endif
real(pReal), dimension (3) :: homogenization_RGC_relaxationVector
integer(pInt), dimension (4), intent(in) :: intFace !< set of interface ID in 4D array (normal and position)
#ifndef NEWSTATE
type(p_vec), intent(in) :: state !< set of global relaxation vectors
#endif
integer(pInt), dimension (3) :: nGDim
integer(pInt) :: &
iNum, &
homID !< homogenization ID
!--------------------------------------------------------------------------------------------------
! collect the interface relaxation vector from the global state array
homogenization_RGC_relaxationVector = 0.0_pReal
nGDim = homogenization_RGC_Ngrains(1:3,homID)
iNum = homogenization_RGC_interface4to1(intFace,homID) ! identify the position of the interface in global state array
#ifdef NEWSTATE
if (iNum > 0_pInt) homogenization_RGC_relaxationVector = homogState(mappingHomogenization(2,ip,el))% &
state((3*iNum-2):(3*iNum),mappingHomogenization(1,ip,el)) ! get the corresponding entries
#else
if (iNum > 0_pInt) homogenization_RGC_relaxationVector = state%p((3*iNum-2):(3*iNum)) ! get the corresponding entries
#endif
end function homogenization_RGC_relaxationVector
!--------------------------------------------------------------------------------------------------
!> @brief identify the normal of an interface
!--------------------------------------------------------------------------------------------------
function homogenization_RGC_interfaceNormal(intFace,ip,el)
use debug, only: &
debug_level, &
debug_homogenization,&
debug_levelExtensive, &
debug_e, &
debug_i
use math, only: &
math_mul33x3
implicit none
real(pReal), dimension (3) :: homogenization_RGC_interfaceNormal
integer(pInt), dimension (4), intent(in) :: intFace !< interface ID in 4D array (normal and position)
integer(pInt), intent(in) :: &
ip, & !< integration point number
el !< element number
integer(pInt) :: nPos
!--------------------------------------------------------------------------------------------------
! get the normal of the interface, identified from the value of intFace(1)
homogenization_RGC_interfaceNormal = 0.0_pReal
nPos = abs(intFace(1)) ! identify the position of the interface in global state array
homogenization_RGC_interfaceNormal(nPos) = real(intFace(1)/abs(intFace(1)),pReal) ! get the normal vector w.r.t. cluster axis
homogenization_RGC_interfaceNormal = &
math_mul33x3(homogenization_RGC_orientation(1:3,1:3,ip,el),homogenization_RGC_interfaceNormal)
! map the normal vector into sample coordinate system (basis)
end function homogenization_RGC_interfaceNormal
!--------------------------------------------------------------------------------------------------
!> @brief collect six faces of a grain in 4D (normal and position)
!--------------------------------------------------------------------------------------------------
function homogenization_RGC_getInterface(iFace,iGrain3)
implicit none
integer(pInt), dimension (4) :: homogenization_RGC_getInterface
integer(pInt), dimension (3), intent(in) :: iGrain3 !< grain ID in 3D array
integer(pInt), intent(in) :: iFace !< face index (1..6) mapped like (-e1,-e2,-e3,+e1,+e2,+e3) or iDir = (-1,-2,-3,1,2,3)
integer(pInt) :: iDir
!* Direction of interface normal
iDir = (int(real(iFace-1_pInt,pReal)/2.0_pReal,pInt)+1_pInt)*(-1_pInt)**iFace
homogenization_RGC_getInterface(1) = iDir
!--------------------------------------------------------------------------------------------------
! identify the interface position by the direction of its normal
homogenization_RGC_getInterface(2:4) = iGrain3
if (iDir < 0_pInt) & ! to have a correlation with coordinate/position in real space
homogenization_RGC_getInterface(1_pInt-iDir) = homogenization_RGC_getInterface(1_pInt-iDir)-1_pInt
end function homogenization_RGC_getInterface
!--------------------------------------------------------------------------------------------------
!> @brief map grain ID from in 1D (global array) to in 3D (local position)
!--------------------------------------------------------------------------------------------------
function homogenization_RGC_grain1to3(grain1,homID)
implicit none
integer(pInt), dimension (3) :: homogenization_RGC_grain1to3
integer(pInt), intent(in) :: &
grain1,& !< grain ID in 1D array
homID !< homogenization ID
integer(pInt), dimension (3) :: nGDim
!--------------------------------------------------------------------------------------------------
! get the grain position
nGDim = homogenization_RGC_Ngrains(1:3,homID)
homogenization_RGC_grain1to3(3) = 1_pInt+(grain1-1_pInt)/(nGDim(1)*nGDim(2))
homogenization_RGC_grain1to3(2) = 1_pInt+mod((grain1-1_pInt)/nGDim(1),nGDim(2))
homogenization_RGC_grain1to3(1) = 1_pInt+mod((grain1-1_pInt),nGDim(1))
end function homogenization_RGC_grain1to3
!--------------------------------------------------------------------------------------------------
!> @brief map grain ID from in 3D (local position) to in 1D (global array)
!--------------------------------------------------------------------------------------------------
pure function homogenization_RGC_grain3to1(grain3,homID)
implicit none
integer(pInt), dimension (3), intent(in) :: grain3 !< grain ID in 3D array (pos.x,pos.y,pos.z)
integer(pInt) :: homogenization_RGC_grain3to1
integer(pInt), dimension (3) :: nGDim
integer(pInt), intent(in) :: homID ! homogenization ID
!--------------------------------------------------------------------------------------------------
! get the grain ID
nGDim = homogenization_RGC_Ngrains(1:3,homID)
homogenization_RGC_grain3to1 = grain3(1) + nGDim(1)*(grain3(2)-1_pInt) + nGDim(1)*nGDim(2)*(grain3(3)-1_pInt)
end function homogenization_RGC_grain3to1
!--------------------------------------------------------------------------------------------------
!> @brief maps interface ID from 4D (normal and local position) into 1D (global array)
!--------------------------------------------------------------------------------------------------
integer(pInt) pure function homogenization_RGC_interface4to1(iFace4D, homID)
implicit none
integer(pInt), dimension (4), intent(in) :: iFace4D !< interface ID in 4D array (n.dir,pos.x,pos.y,pos.z)
integer(pInt), dimension (3) :: nGDim,nIntFace
integer(pInt), intent(in) :: homID !< homogenization ID
nGDim = homogenization_RGC_Ngrains(1:3,homID)
!--------------------------------------------------------------------------------------------------
! compute the total number of interfaces, which ...
nIntFace(1) = (nGDim(1)-1_pInt)*nGDim(2)*nGDim(3) ! ... normal //e1
nIntFace(2) = nGDim(1)*(nGDim(2)-1_pInt)*nGDim(3) ! ... normal //e2
nIntFace(3) = nGDim(1)*nGDim(2)*(nGDim(3)-1_pInt) ! ... normal //e3
homogenization_RGC_interface4to1 = -1_pInt
!--------------------------------------------------------------------------------------------------
! get the corresponding interface ID in 1D global array
if (abs(iFace4D(1)) == 1_pInt) then ! interface with normal //e1
homogenization_RGC_interface4to1 = iFace4D(3) + nGDim(2)*(iFace4D(4)-1_pInt) &
+ nGDim(2)*nGDim(3)*(iFace4D(2)-1_pInt)
if ((iFace4D(2) == 0_pInt) .or. (iFace4D(2) == nGDim(1))) homogenization_RGC_interface4to1 = 0_pInt
elseif (abs(iFace4D(1)) == 2_pInt) then ! interface with normal //e2
homogenization_RGC_interface4to1 = iFace4D(4) + nGDim(3)*(iFace4D(2)-1_pInt) &
+ nGDim(3)*nGDim(1)*(iFace4D(3)-1_pInt) + nIntFace(1)
if ((iFace4D(3) == 0_pInt) .or. (iFace4D(3) == nGDim(2))) homogenization_RGC_interface4to1 = 0_pInt
elseif (abs(iFace4D(1)) == 3_pInt) then ! interface with normal //e3
homogenization_RGC_interface4to1 = iFace4D(2) + nGDim(1)*(iFace4D(3)-1_pInt) &
+ nGDim(1)*nGDim(2)*(iFace4D(4)-1_pInt) + nIntFace(1) + nIntFace(2)
if ((iFace4D(4) == 0_pInt) .or. (iFace4D(4) == nGDim(3))) homogenization_RGC_interface4to1 = 0_pInt
endif
end function homogenization_RGC_interface4to1
!--------------------------------------------------------------------------------------------------
!> @brief maps interface ID from 1D (global array) into 4D (normal and local position)
!--------------------------------------------------------------------------------------------------
function homogenization_RGC_interface1to4(iFace1D, homID)
implicit none
integer(pInt), dimension (4) :: homogenization_RGC_interface1to4
integer(pInt), intent(in) :: iFace1D !< interface ID in 1D array
integer(pInt), dimension (3) :: nGDim,nIntFace
integer(pInt), intent(in) :: homID !< homogenization ID
nGDim = homogenization_RGC_Ngrains(:,homID)
!--------------------------------------------------------------------------------------------------
! compute the total number of interfaces, which ...
nIntFace(1) = (nGDim(1)-1_pInt)*nGDim(2)*nGDim(3) ! ... normal //e1
nIntFace(2) = nGDim(1)*(nGDim(2)-1_pInt)*nGDim(3) ! ... normal //e2
nIntFace(3) = nGDim(1)*nGDim(2)*(nGDim(3)-1_pInt) ! ... normal //e3
!--------------------------------------------------------------------------------------------------
! get the corresponding interface ID in 4D (normal and local position)
if (iFace1D > 0 .and. iFace1D <= nIntFace(1)) then ! interface with normal //e1
homogenization_RGC_interface1to4(1) = 1_pInt
homogenization_RGC_interface1to4(3) = mod((iFace1D-1_pInt),nGDim(2))+1_pInt
homogenization_RGC_interface1to4(4) = mod(&
int(&
real(iFace1D-1_pInt,pReal)/&
real(nGDim(2),pReal)&
,pInt)&
,nGDim(3))+1_pInt
homogenization_RGC_interface1to4(2) = int(&
real(iFace1D-1_pInt,pReal)/&
real(nGDim(2),pReal)/&
real(nGDim(3),pReal)&
,pInt)+1_pInt
elseif (iFace1D > nIntFace(1) .and. iFace1D <= (nIntFace(2) + nIntFace(1))) then ! interface with normal //e2
homogenization_RGC_interface1to4(1) = 2_pInt
homogenization_RGC_interface1to4(4) = mod((iFace1D-nIntFace(1)-1_pInt),nGDim(3))+1_pInt
homogenization_RGC_interface1to4(2) = mod(&
int(&
real(iFace1D-nIntFace(1)-1_pInt,pReal)/&
real(nGDim(3),pReal)&
,pInt)&
,nGDim(1))+1_pInt
homogenization_RGC_interface1to4(3) = int(&
real(iFace1D-nIntFace(1)-1_pInt,pReal)/&
real(nGDim(3),pReal)/&
real(nGDim(1),pReal)&
,pInt)+1_pInt
elseif (iFace1D > nIntFace(2) + nIntFace(1) .and. iFace1D <= (nIntFace(3) + nIntFace(2) + nIntFace(1))) then ! interface with normal //e3
homogenization_RGC_interface1to4(1) = 3_pInt
homogenization_RGC_interface1to4(2) = mod((iFace1D-nIntFace(2)-nIntFace(1)-1_pInt),nGDim(1))+1_pInt
homogenization_RGC_interface1to4(3) = mod(&
int(&
real(iFace1D-nIntFace(2)-nIntFace(1)-1_pInt,pReal)/&
real(nGDim(1),pReal)&
,pInt)&
,nGDim(2))+1_pInt
homogenization_RGC_interface1to4(4) = int(&
real(iFace1D-nIntFace(2)-nIntFace(1)-1_pInt,pReal)/&
real(nGDim(1),pReal)/&
real(nGDim(2),pReal)&
,pInt)+1_pInt
endif
end function homogenization_RGC_interface1to4
!--------------------------------------------------------------------------------------------------
!> @brief calculating the grain deformation gradient (the same with
! homogenization_RGC_partionDeformation, but used only for perturbation scheme)
!--------------------------------------------------------------------------------------------------
#ifdef NEWSTATE
subroutine homogenization_RGC_grainDeformation(F, avgF, ip, el)
#else
subroutine homogenization_RGC_grainDeformation(F, avgF, state, ip, el)
#endif
use prec, only: &
p_vec
use mesh, only: &
mesh_element
use material, only: &
homogenization_maxNgrains,&
homogenization_Ngrains, &
#ifdef NEWSTATE
homogState, &
mappingHomogenization, &
#endif
homogenization_typeInstance
implicit none
real(pReal), dimension (3,3,homogenization_maxNgrains), intent(out) :: F !< partioned F per grain
real(pReal), dimension (3,3), intent(in) :: avgF !<
#ifndef NEWSTATE
type(p_vec), intent(in) :: state
#endif
integer(pInt), intent(in) :: &
el, & !< element number
ip !< integration point number
real(pReal), dimension (3) :: aVect,nVect
integer(pInt), dimension (4) :: intFace
integer(pInt), dimension (3) :: iGrain3
integer(pInt) :: homID, iGrain,iFace,i,j
integer(pInt), parameter :: nFace = 6_pInt
!--------------------------------------------------------------------------------------------------
! compute the deformation gradient of individual grains due to relaxations
homID = homogenization_typeInstance(mesh_element(3,el))
F = 0.0_pReal
do iGrain = 1_pInt,homogenization_Ngrains(mesh_element(3,el))
iGrain3 = homogenization_RGC_grain1to3(iGrain,homID)
do iFace = 1_pInt,nFace
intFace = homogenization_RGC_getInterface(iFace,iGrain3)
#ifdef NEWSTATE
aVect = homogenization_RGC_relaxationVector(intFace,homID, ip, el)
#else
aVect = homogenization_RGC_relaxationVector(intFace,state,homID)
#endif
nVect = homogenization_RGC_interfaceNormal(intFace,ip,el)
forall (i=1_pInt:3_pInt,j=1_pInt:3_pInt) &
F(i,j,iGrain) = F(i,j,iGrain) + aVect(i)*nVect(j) ! effective relaxations
enddo
F(1:3,1:3,iGrain) = F(1:3,1:3,iGrain) + avgF ! relaxed deformation gradient
enddo
end subroutine homogenization_RGC_grainDeformation
end module homogenization_RGC
Reference in New Issue View Git Blame Copy Permalink