2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
2018-08-28 16:45:38 +05:30
!> @author Martin Diehl, Max-Planck-Institut für Eisenforschung GmbH
2013-02-21 03:36:15 +05:30
!> @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
2018-08-28 16:45:38 +05:30
!> Nconstituents is defined as p x q x r (cluster)
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
2019-05-18 10:53:46 +05:30
submodule ( homogenization ) homogenization_mech_RGC
2020-01-29 17:46:49 +05:30
use rotations
2018-08-28 16:27:22 +05:30
2019-06-15 21:57:38 +05:30
type :: tParameters
integer , dimension ( : ) , allocatable :: &
Nconstituents
real ( pReal ) :: &
xiAlpha , &
ciAlpha
real ( pReal ) , dimension ( : ) , allocatable :: &
dAlpha , &
angles
integer :: &
of_debug = 0
2020-02-15 03:51:58 +05:30
character ( len = pStringLen ) , allocatable , dimension ( : ) :: &
output
2019-06-15 21:57:38 +05:30
end type tParameters
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
type :: tRGCstate
real ( pReal ) , pointer , dimension ( : ) :: &
work , &
penaltyEnergy
real ( pReal ) , pointer , dimension ( : , : ) :: &
relaxationVector
end type tRGCstate
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
type :: tRGCdependentState
real ( pReal ) , allocatable , dimension ( : ) :: &
volumeDiscrepancy , &
relaxationRate_avg , &
relaxationRate_max
real ( pReal ) , allocatable , dimension ( : , : ) :: &
mismatch
real ( pReal ) , allocatable , dimension ( : , : , : ) :: &
orientation
end type tRGCdependentState
2020-03-17 02:09:53 +05:30
2020-03-17 05:09:32 +05:30
type :: tNumerics_RGC
2020-03-17 02:09:53 +05:30
real ( pReal ) :: &
atol , & !< absolute tolerance of RGC residuum
rtol , & !< relative tolerance of RGC residuum
absMax , & !< absolute maximum of RGC residuum
relMax , & !< relative maximum of RGC residuum
pPert , & !< perturbation for computing RGC penalty tangent
xSmoo , & !< RGC penalty smoothing parameter (hyperbolic tangent)
viscPower , & !< power (sensitivity rate) of numerical viscosity in RGC scheme, Default 1.0e0: Newton viscosity (linear model)
viscModus , & !< stress modulus of RGC numerical viscosity, Default 0.0e0: No viscosity is applied
refRelaxRate , & !< reference relaxation rate in RGC viscosity
maxdRelax , & !< threshold of maximum relaxation vector increment (if exceed this then cutback)
maxVolDiscr , & !< threshold of maximum volume discrepancy allowed
volDiscrMod , & !< stiffness of RGC volume discrepancy (zero = without volume discrepancy constraint)
volDiscrPow !< powerlaw penalty for volume discrepancy
2020-03-17 05:09:32 +05:30
end type tNumerics_RGC
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
type ( tparameters ) , dimension ( : ) , allocatable :: &
param
type ( tRGCstate ) , dimension ( : ) , allocatable :: &
state , &
state0
type ( tRGCdependentState ) , dimension ( : ) , allocatable :: &
dependentState
2020-03-17 05:09:32 +05:30
type ( tNumerics_RGC ) :: &
2020-03-17 02:09:53 +05:30
num ! numerics parameters. Better name?
2018-11-04 02:06:49 +05:30
2013-02-21 03:36:15 +05:30
contains
!--------------------------------------------------------------------------------------------------
2018-04-17 18:12:04 +05:30
!> @brief allocates all necessary fields, reads information from material configuration file
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
2019-05-18 10:53:46 +05:30
module subroutine mech_RGC_init
2012-03-09 01:55:28 +05:30
2019-06-15 21:57:38 +05:30
integer :: &
Ninstance , &
2020-02-15 03:51:58 +05:30
h , &
2019-06-15 21:57:38 +05:30
NofMyHomog , &
sizeState , nIntFaceTot
2020-03-17 02:09:53 +05:30
2020-02-15 03:51:58 +05:30
write ( 6 , '(/,a)' ) ' <<<+- homogenization_' / / HOMOGENIZATION_RGC_label / / ' init -+>>>' ; flush ( 6 )
2020-03-17 02:09:53 +05:30
2020-02-15 03:51:58 +05:30
write ( 6 , '(/,a)' ) ' Tjahjanto et al., International Journal of Material Forming 2(1):939– 942, 2009'
write ( 6 , '(a)' ) ' https://doi.org/10.1007/s12289-009-0619-1'
2020-03-17 02:09:53 +05:30
2020-02-15 03:51:58 +05:30
write ( 6 , '(/,a)' ) ' Tjahjanto et al., Modelling and Simulation in Materials Science and Engineering 18:015006, 2010'
write ( 6 , '(a)' ) ' https://doi.org/10.1088/0965-0393/18/1/015006'
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
Ninstance = count ( homogenization_type == HOMOGENIZATION_RGC_ID )
if ( iand ( debug_level ( debug_HOMOGENIZATION ) , debug_levelBasic ) / = 0 ) &
write ( 6 , '(a16,1x,i5,/)' ) '# instances:' , Ninstance
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
allocate ( param ( Ninstance ) )
allocate ( state ( Ninstance ) )
allocate ( state0 ( Ninstance ) )
allocate ( dependentState ( Ninstance ) )
2020-03-17 02:09:53 +05:30
num % atol = config_numerics % getFloat ( 'atol_rgc' , defaultVal = 1.0e+4_pReal )
num % rtol = config_numerics % getFloat ( 'rtol_rgc' , defaultVal = 1.0e-3_pReal )
num % absMax = config_numerics % getFloat ( 'amax_rgc' , defaultVal = 1.0e+10_pReal )
num % relMax = config_numerics % getFloat ( 'rmax_rgc' , defaultVal = 1.0e+2_pReal )
num % pPert = config_numerics % getFloat ( 'perturbpenalty_rgc' , defaultVal = 1.0e-7_pReal )
num % xSmoo = config_numerics % getFloat ( 'relvantmismatch_rgc' , defaultVal = 1.0e-5_pReal )
num % viscPower = config_numerics % getFloat ( 'viscositypower_rgc' , defaultVal = 1.0e+0_pReal )
num % viscModus = config_numerics % getFloat ( 'viscositymodulus_rgc' , defaultVal = 0.0e+0_pReal )
num % refRelaxRate = config_numerics % getFloat ( 'refrelaxationrate_rgc' , defaultVal = 1.0e-3_pReal )
num % maxdRelax = config_numerics % getFloat ( 'maxrelaxationrate_rgc' , defaultVal = 1.0e+0_pReal )
num % maxVolDiscr = config_numerics % getFloat ( 'maxvoldiscrepancy_rgc' , defaultVal = 1.0e-5_pReal )
num % volDiscrMod = config_numerics % getFloat ( 'voldiscrepancymod_rgc' , defaultVal = 1.0e+12_pReal )
num % volDiscrPow = config_numerics % getFloat ( 'dicrepancypower_rgc' , defaultVal = 5.0_pReal )
if ( num % atol < = 0.0_pReal ) call IO_error ( 301 , ext_msg = 'absTol_RGC' )
if ( num % rtol < = 0.0_pReal ) call IO_error ( 301 , ext_msg = 'relTol_RGC' )
if ( num % absMax < = 0.0_pReal ) call IO_error ( 301 , ext_msg = 'absMax_RGC' )
if ( num % relMax < = 0.0_pReal ) call IO_error ( 301 , ext_msg = 'relMax_RGC' )
if ( num % pPert < = 0.0_pReal ) call IO_error ( 301 , ext_msg = 'pPert_RGC' )
if ( num % xSmoo < = 0.0_pReal ) call IO_error ( 301 , ext_msg = 'xSmoo_RGC' )
if ( num % viscPower < 0.0_pReal ) call IO_error ( 301 , ext_msg = 'viscPower_RGC' )
if ( num % viscModus < 0.0_pReal ) call IO_error ( 301 , ext_msg = 'viscModus_RGC' )
if ( num % refRelaxRate < = 0.0_pReal ) call IO_error ( 301 , ext_msg = 'refRelaxRate_RGC' )
if ( num % maxdRelax < = 0.0_pReal ) call IO_error ( 301 , ext_msg = 'maxdRelax_RGC' )
if ( num % maxVolDiscr < = 0.0_pReal ) call IO_error ( 301 , ext_msg = 'maxVolDiscr_RGC' )
if ( num % volDiscrMod < 0.0_pReal ) call IO_error ( 301 , ext_msg = 'volDiscrMod_RGC' )
if ( num % volDiscrPow < = 0.0_pReal ) call IO_error ( 301 , ext_msg = 'volDiscrPw_RGC' )
2019-06-15 21:57:38 +05:30
do h = 1 , size ( homogenization_type )
if ( homogenization_type ( h ) / = HOMOGENIZATION_RGC_ID ) cycle
associate ( prm = > param ( homogenization_typeInstance ( h ) ) , &
stt = > state ( homogenization_typeInstance ( h ) ) , &
st0 = > state0 ( homogenization_typeInstance ( h ) ) , &
dst = > dependentState ( homogenization_typeInstance ( h ) ) , &
config = > config_homogenization ( h ) )
2020-03-17 02:09:53 +05:30
2019-01-13 02:34:03 +05:30
#ifdef DEBUG
2019-06-15 21:57:38 +05:30
if ( h == material_homogenizationAt ( debug_e ) ) then
2020-01-23 17:18:18 +05:30
prm % of_debug = material_homogenizationMemberAt ( debug_i , debug_e )
2019-06-15 21:57:38 +05:30
endif
2019-01-13 02:34:03 +05:30
#endif
2018-08-31 04:54:28 +05:30
2020-03-17 02:09:53 +05:30
prm % output = config % getStrings ( '(output)' , defaultVal = emptyStringArray )
2019-06-15 21:57:38 +05:30
prm % Nconstituents = config % getInts ( 'clustersize' , requiredSize = 3 )
if ( homogenization_Ngrains ( h ) / = product ( prm % Nconstituents ) ) &
call IO_error ( 211 , ext_msg = 'clustersize (' / / HOMOGENIZATION_RGC_label / / ')' )
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
prm % xiAlpha = config % getFloat ( 'scalingparameter' )
prm % ciAlpha = config % getFloat ( 'overproportionality' )
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
prm % dAlpha = config % getFloats ( 'grainsize' , requiredSize = 3 )
prm % angles = config % getFloats ( 'clusterorientation' , requiredSize = 3 )
2020-02-15 03:51:58 +05:30
2019-06-15 21:57:38 +05:30
NofMyHomog = count ( material_homogenizationAt == h )
nIntFaceTot = 3 * ( ( prm % Nconstituents ( 1 ) - 1 ) * prm % Nconstituents ( 2 ) * prm % Nconstituents ( 3 ) &
+ prm % Nconstituents ( 1 ) * ( prm % Nconstituents ( 2 ) - 1 ) * prm % Nconstituents ( 3 ) &
+ prm % Nconstituents ( 1 ) * prm % Nconstituents ( 2 ) * ( prm % Nconstituents ( 3 ) - 1 ) )
sizeState = nIntFaceTot &
+ size ( [ 'avg constitutive work ' , 'average penalty energy' ] )
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
homogState ( h ) % sizeState = sizeState
allocate ( homogState ( h ) % state0 ( sizeState , NofMyHomog ) , source = 0.0_pReal )
allocate ( homogState ( h ) % subState0 ( sizeState , NofMyHomog ) , source = 0.0_pReal )
allocate ( homogState ( h ) % state ( sizeState , NofMyHomog ) , source = 0.0_pReal )
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
stt % relaxationVector = > homogState ( h ) % state ( 1 : nIntFaceTot , : )
st0 % relaxationVector = > homogState ( h ) % state0 ( 1 : nIntFaceTot , : )
stt % work = > homogState ( h ) % state ( nIntFaceTot + 1 , : )
stt % penaltyEnergy = > homogState ( h ) % state ( nIntFaceTot + 2 , : )
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
allocate ( dst % volumeDiscrepancy ( NofMyHomog ) )
allocate ( dst % relaxationRate_avg ( NofMyHomog ) )
allocate ( dst % relaxationRate_max ( NofMyHomog ) )
allocate ( dst % mismatch ( 3 , NofMyHomog ) )
2018-11-04 02:06:49 +05:30
!--------------------------------------------------------------------------------------------------
2019-01-13 02:34:03 +05:30
! assigning cluster orientations
2020-01-29 13:55:39 +05:30
dependentState ( homogenization_typeInstance ( h ) ) % orientation = spread ( eu2om ( prm % angles * inRad ) , 3 , NofMyHomog )
!dst%orientation = spread(eu2om(prm%angles*inRad),3,NofMyHomog) ifort version 18.0.1 crashes (for whatever reason)
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
end associate
2020-03-17 02:09:53 +05:30
enddo
2019-05-18 10:53:46 +05:30
end subroutine mech_RGC_init
2009-10-12 21:31:49 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief partitions the deformation gradient onto the constituents
!--------------------------------------------------------------------------------------------------
2019-05-18 10:53:46 +05:30
module subroutine mech_RGC_partitionDeformation ( F , avgF , instance , of )
2009-10-12 21:31:49 +05:30
2019-06-15 21:57:38 +05:30
real ( pReal ) , dimension ( : , : , : ) , intent ( out ) :: F !< partioned F per grain
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: avgF !< averaged F
integer , intent ( in ) :: &
instance , &
of
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
real ( pReal ) , dimension ( 3 ) :: aVect , nVect
integer , dimension ( 4 ) :: intFace
integer , dimension ( 3 ) :: iGrain3
integer :: iGrain , iFace , i , j
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
associate ( prm = > param ( instance ) )
2020-03-17 02:09:53 +05:30
2019-01-13 13:44:23 +05:30
!--------------------------------------------------------------------------------------------------
! compute the deformation gradient of individual grains due to relaxations
2019-06-15 21:57:38 +05:30
F = 0.0_pReal
do iGrain = 1 , product ( prm % Nconstituents )
iGrain3 = grain1to3 ( iGrain , prm % Nconstituents )
do iFace = 1 , 6
intFace = getInterface ( iFace , iGrain3 ) ! identifying 6 interfaces of each grain
aVect = relaxationVector ( intFace , instance , of ) ! get the relaxation vectors for each interface from global relaxation vector array
nVect = interfaceNormal ( intFace , instance , of )
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
F ( 1 : 3 , 1 : 3 , iGrain ) = F ( 1 : 3 , 1 : 3 , iGrain ) + avgF ! resulting relaxed deformation gradient
2010-11-26 17:20:20 +05:30
2019-01-13 02:34:03 +05:30
#ifdef DEBUG
2019-06-15 21:57:38 +05:30
if ( iand ( debug_level ( debug_homogenization ) , debug_levelExtensive ) / = 0 ) then
write ( 6 , '(1x,a32,1x,i3)' ) 'Deformation gradient of grain: ' , iGrain
do i = 1 , 3
write ( 6 , '(1x,3(e15.8,1x))' ) ( F ( i , j , iGrain ) , j = 1 , 3 )
enddo
write ( 6 , * ) ' '
flush ( 6 )
endif
2019-01-13 02:34:03 +05:30
#endif
2019-06-15 21:57:38 +05:30
enddo
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
end associate
2009-10-12 21:31:49 +05:30
2019-05-18 10:53:46 +05:30
end subroutine mech_RGC_partitionDeformation
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-17 02:09:53 +05:30
!> @brief update the internal state of the homogenization scheme and tell whether "done" and
2013-02-21 03:36:15 +05:30
! "happy" with result
!--------------------------------------------------------------------------------------------------
2019-05-18 13:17:20 +05:30
module procedure mech_RGC_updateState
2018-11-03 21:20:43 +05:30
2019-06-15 21:57:38 +05:30
integer , dimension ( 4 ) :: intFaceN , intFaceP , faceID
integer , dimension ( 3 ) :: nGDim , iGr3N , iGr3P
integer :: instance , iNum , i , j , nIntFaceTot , iGrN , iGrP , iMun , iFace , k , l , ipert , iGrain , nGrain , of
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
real ( pReal ) :: residMax , stresMax
logical :: error
real ( pReal ) , dimension ( : , : ) , allocatable :: tract , jmatrix , jnverse , smatrix , pmatrix , rmatrix
real ( pReal ) , dimension ( : ) , allocatable :: resid , relax , p_relax , p_resid , drelax
2019-01-13 14:03:47 +05:30
#ifdef DEBUG
2019-06-15 21:57:38 +05:30
integer , dimension ( 3 ) :: stresLoc
integer , dimension ( 2 ) :: residLoc
2019-01-13 14:03:47 +05:30
#endif
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
zeroTimeStep : if ( dEq0 ( dt ) ) then
mech_RGC_updateState = . true . ! pretend everything is fine and return
2020-03-17 02:09:53 +05:30
return
2019-06-15 21:57:38 +05:30
endif zeroTimeStep
2009-10-12 21:31:49 +05:30
2019-06-15 21:57:38 +05:30
instance = homogenization_typeInstance ( material_homogenizationAt ( el ) )
2020-01-23 17:18:18 +05:30
of = material_homogenizationMemberAt ( ip , el )
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
associate ( stt = > state ( instance ) , st0 = > state0 ( instance ) , dst = > dependentState ( instance ) , prm = > param ( instance ) )
2020-03-17 02:09:53 +05:30
2019-01-13 03:37:35 +05:30
!--------------------------------------------------------------------------------------------------
! get the dimension of the cluster (grains and interfaces)
2019-06-15 21:57:38 +05:30
nGDim = prm % Nconstituents
nGrain = product ( nGDim )
nIntFaceTot = ( nGDim ( 1 ) - 1 ) * nGDim ( 2 ) * nGDim ( 3 ) &
+ nGDim ( 1 ) * ( nGDim ( 2 ) - 1 ) * nGDim ( 3 ) &
+ nGDim ( 1 ) * nGDim ( 2 ) * ( nGDim ( 3 ) - 1 )
2009-10-12 21:31:49 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! allocate the size of the global relaxation arrays/jacobian matrices depending on the size of the cluster
2019-06-15 21:57:38 +05:30
allocate ( resid ( 3 * nIntFaceTot ) , source = 0.0_pReal )
allocate ( tract ( nIntFaceTot , 3 ) , source = 0.0_pReal )
relax = stt % relaxationVector ( : , of )
drelax = stt % relaxationVector ( : , of ) - st0 % relaxationVector ( : , of )
2020-03-17 02:09:53 +05:30
2019-01-13 02:34:03 +05:30
#ifdef DEBUG
2019-06-15 21:57:38 +05:30
if ( iand ( debug_level ( debug_homogenization ) , debug_levelExtensive ) / = 0 ) then
write ( 6 , '(1x,a30)' ) 'Obtained state: '
do i = 1 , size ( stt % relaxationVector ( : , of ) )
write ( 6 , '(1x,2(e15.8,1x))' ) stt % relaxationVector ( i , of )
enddo
write ( 6 , * ) ' '
endif
2019-01-13 02:34:03 +05:30
#endif
2009-10-12 21:31:49 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! computing interface mismatch and stress penalty tensor for all interfaces of all grains
2019-06-15 21:57:38 +05:30
call stressPenalty ( R , NN , avgF , F , ip , el , instance , of )
2009-12-16 21:50:53 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-17 02:09:53 +05:30
! calculating volume discrepancy and stress penalty related to overall volume discrepancy
2019-06-15 21:57:38 +05:30
call volumePenalty ( D , dst % volumeDiscrepancy ( of ) , avgF , F , nGrain , instance , of )
2009-12-16 21:50:53 +05:30
2019-01-13 02:34:03 +05:30
#ifdef DEBUG
2019-06-15 21:57:38 +05:30
if ( iand ( debug_level ( debug_homogenization ) , debug_levelExtensive ) / = 0 ) then
do iGrain = 1 , 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 , 3
write ( 6 , '(1x,3(e15.8,1x),1x,3(e15.8,1x),1x,3(e15.8,1x))' ) ( P ( i , j , iGrain ) , j = 1 , 3 ) , &
( R ( i , j , iGrain ) , j = 1 , 3 ) , &
( D ( i , j , iGrain ) , j = 1 , 3 )
enddo
write ( 6 , * ) ' '
enddo
endif
2019-01-13 02:34:03 +05:30
#endif
2009-10-12 21:31:49 +05:30
2014-09-19 23:29:06 +05:30
!------------------------------------------------------------------------------------------------
2013-02-21 03:36:15 +05:30
! computing the residual stress from the balance of traction at all (interior) interfaces
2019-06-15 21:57:38 +05:30
do iNum = 1 , nIntFaceTot
faceID = interface1to4 ( iNum , param ( instance ) % Nconstituents ) ! identifying the interface ID in local coordinate system (4-dimensional index)
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! identify the left/bottom/back grain (-|N)
2019-06-15 21:57:38 +05:30
iGr3N = faceID ( 2 : 4 ) ! identifying the grain ID in local coordinate system (3-dimensional index)
iGrN = grain3to1 ( iGr3N , param ( instance ) % Nconstituents ) ! translate the local grain ID into global coordinate system (1-dimensional index)
intFaceN = getInterface ( 2 * faceID ( 1 ) , iGr3N )
normN = interfaceNormal ( intFaceN , instance , of )
2020-03-17 02:09:53 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! identify the right/up/front grain (+|P)
2019-06-15 21:57:38 +05:30
iGr3P = iGr3N
iGr3P ( faceID ( 1 ) ) = iGr3N ( faceID ( 1 ) ) + 1 ! identifying the grain ID in local coordinate system (3-dimensional index)
iGrP = grain3to1 ( iGr3P , param ( instance ) % Nconstituents ) ! translate the local grain ID into global coordinate system (1-dimensional index)
intFaceP = getInterface ( 2 * faceID ( 1 ) - 1 , iGr3P )
normP = interfaceNormal ( intFaceP , instance , of )
2010-11-26 17:20:20 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! compute the residual of traction at the interface (in local system, 4-dimensional index)
2019-06-15 21:57:38 +05:30
do i = 1 , 3
2020-03-17 02:09:53 +05:30
tract ( iNum , i ) = sign ( num % viscModus * ( abs ( drelax ( i + 3 * ( iNum - 1 ) ) ) / ( num % refRelaxRate * dt ) ) ** num % viscPower , &
2019-06-15 21:57:38 +05:30
drelax ( i + 3 * ( iNum - 1 ) ) ) ! contribution from the relaxation viscosity
do j = 1 , 3
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 * ( iNum - 1 ) ) = tract ( iNum , i ) ! translate the local residual into global 1-dimensional residual array
enddo
enddo
2020-03-17 02:09:53 +05:30
2019-01-13 02:34:03 +05:30
#ifdef DEBUG
2019-06-15 21:57:38 +05:30
if ( iand ( debug_level ( debug_homogenization ) , debug_levelExtensive ) / = 0 ) then
write ( 6 , '(1x,a30,1x,i3)' ) 'Traction at interface: ' , iNum
write ( 6 , '(1x,3(e15.8,1x))' ) ( tract ( iNum , j ) , j = 1 , 3 )
write ( 6 , * ) ' '
endif
2019-01-13 02:34:03 +05:30
#endif
2019-06-15 21:57:38 +05:30
enddo
2009-10-12 21:31:49 +05:30
2013-12-18 12:58:01 +05:30
!--------------------------------------------------------------------------------------------------
2013-02-21 03:36:15 +05:30
! convergence check for stress residual
2019-06-15 21:57:38 +05:30
stresMax = maxval ( abs ( P ) ) ! get the maximum of first Piola-Kirchhoff (material) stress
residMax = maxval ( abs ( tract ) ) ! get the maximum of the residual
2010-11-26 17:20:20 +05:30
2019-01-13 02:34:03 +05:30
#ifdef DEBUG
2019-06-15 21:57:38 +05:30
if ( iand ( debug_level ( debug_homogenization ) , debug_levelExtensive ) / = 0 . and . prm % of_debug == of ) then
stresLoc = maxloc ( abs ( P ) )
residLoc = maxloc ( abs ( tract ) )
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 )
endif
2019-01-13 02:34:03 +05:30
#endif
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
mech_RGC_updateState = . false .
2020-03-17 02:09:53 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! If convergence reached => done and happy
2020-03-17 02:09:53 +05:30
if ( residMax < num % rtol * stresMax . or . residMax < num % atol ) then
2019-06-15 21:57:38 +05:30
mech_RGC_updateState = . true .
2019-01-13 03:37:35 +05:30
#ifdef DEBUG
2019-06-15 21:57:38 +05:30
if ( iand ( debug_level ( debug_homogenization ) , debug_levelExtensive ) / = 0 . and . prm % of_debug == of ) &
write ( 6 , '(1x,a55,/)' ) '... done and happy' ; flush ( 6 )
2019-01-13 03:37:35 +05:30
#endif
2010-11-26 17:20:20 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! compute/update the state for postResult, i.e., all energy densities computed by time-integration
2019-06-15 21:57:38 +05:30
do iGrain = 1 , product ( prm % Nconstituents )
do i = 1 , 3 ; do j = 1 , 3
stt % work ( of ) = stt % work ( of ) &
+ P ( i , j , iGrain ) * ( F ( i , j , iGrain ) - F0 ( i , j , iGrain ) ) / real ( nGrain , pReal )
stt % penaltyEnergy ( of ) = stt % penaltyEnergy ( of ) &
+ R ( i , j , iGrain ) * ( F ( i , j , iGrain ) - F0 ( i , j , iGrain ) ) / real ( nGrain , pReal )
enddo ; enddo
enddo
dst % mismatch ( 1 : 3 , of ) = sum ( NN , 2 ) / real ( nGrain , pReal )
dst % relaxationRate_avg ( of ) = sum ( abs ( drelax ) ) / dt / real ( 3 * nIntFaceTot , pReal )
dst % relaxationRate_max ( of ) = maxval ( abs ( drelax ) ) / dt
2020-03-17 02:09:53 +05:30
2019-01-13 03:37:35 +05:30
#ifdef DEBUG
2019-06-15 21:57:38 +05:30
if ( iand ( debug_level ( debug_homogenization ) , debug_levelExtensive ) / = 0 . and . prm % of_debug == of ) then
write ( 6 , '(1x,a30,1x,e15.8)' ) 'Constitutive work: ' , stt % work ( of )
write ( 6 , '(1x,a30,3(1x,e15.8))' ) 'Magnitude mismatch: ' , dst % mismatch ( 1 , of ) , &
dst % mismatch ( 2 , of ) , &
dst % mismatch ( 3 , of )
write ( 6 , '(1x,a30,1x,e15.8)' ) 'Penalty energy: ' , stt % penaltyEnergy ( of )
write ( 6 , '(1x,a30,1x,e15.8,/)' ) 'Volume discrepancy: ' , dst % volumeDiscrepancy ( of )
write ( 6 , '(1x,a30,1x,e15.8)' ) 'Maximum relaxation rate: ' , dst % relaxationRate_max ( of )
write ( 6 , '(1x,a30,1x,e15.8,/)' ) 'Average relaxation rate: ' , dst % relaxationRate_avg ( of )
flush ( 6 )
endif
2019-01-13 03:37:35 +05:30
#endif
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
return
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! if residual blows-up => done but unhappy
2020-03-17 02:09:53 +05:30
elseif ( residMax > num % relMax * stresMax . or . residMax > num % absMax ) then ! try to restart when residual blows up exceeding maximum bound
2019-06-15 21:57:38 +05:30
mech_RGC_updateState = [ . true . , . false . ] ! with direct cut-back
2020-03-17 02:09:53 +05:30
2019-01-13 03:37:35 +05:30
#ifdef DEBUG
2019-06-15 21:57:38 +05:30
if ( iand ( debug_level ( debug_homogenization ) , debug_levelExtensive ) / = 0 . and . prm % of_debug == of ) &
write ( 6 , '(1x,a,/)' ) '... broken' ; flush ( 6 )
2019-01-13 03:37:35 +05:30
#endif
return
2020-03-17 02:09:53 +05:30
2013-02-21 03:36:15 +05:30
else ! proceed with computing the Jacobian and state update
2019-01-13 03:37:35 +05:30
#ifdef DEBUG
2019-06-15 21:57:38 +05:30
if ( iand ( debug_level ( debug_homogenization ) , debug_levelExtensive ) / = 0 . and . prm % of_debug == of ) &
write ( 6 , '(1x,a,/)' ) '... not yet done' ; flush ( 6 )
2019-01-13 03:37:35 +05:30
#endif
2020-03-17 02:09:53 +05:30
2009-10-12 21:31:49 +05:30
endif
2014-09-19 23:29:06 +05:30
!---------------------------------------------------------------------------------------------------
2020-03-17 02:09:53 +05:30
! construct the global Jacobian matrix for updating the global relaxation vector array when
2013-02-21 03:36:15 +05:30
! convergence is not yet reached ...
!--------------------------------------------------------------------------------------------------
! ... of the constitutive stress tangent, assembled from dPdF or material constitutive model "smatrix"
2019-06-15 21:57:38 +05:30
allocate ( smatrix ( 3 * nIntFaceTot , 3 * nIntFaceTot ) , source = 0.0_pReal )
do iNum = 1 , nIntFaceTot
faceID = interface1to4 ( iNum , param ( instance ) % Nconstituents ) ! assembling of local dPdF into global Jacobian matrix
2020-03-17 02:09:53 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! identify the left/bottom/back grain (-|N)
2019-06-15 21:57:38 +05:30
iGr3N = faceID ( 2 : 4 ) ! identifying the grain ID in local coordinate sytem
iGrN = grain3to1 ( iGr3N , param ( instance ) % Nconstituents ) ! translate into global grain ID
intFaceN = getInterface ( 2 * faceID ( 1 ) , iGr3N ) ! identifying the connecting interface in local coordinate system
normN = interfaceNormal ( intFaceN , instance , of )
do iFace = 1 , 6
intFaceN = getInterface ( iFace , iGr3N ) ! identifying all interfaces that influence relaxation of the above interface
mornN = interfaceNormal ( intFaceN , instance , of )
iMun = interface4to1 ( intFaceN , param ( instance ) % Nconstituents ) ! 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
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
2014-08-21 23:18:20 +05:30
! projecting the material tangent dPdF into the interface
! to obtain the Jacobian matrix contribution of dPdF
2019-06-15 21:57:38 +05:30
endif
enddo
2020-03-17 02:09:53 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! identify the right/up/front grain (+|P)
2019-06-15 21:57:38 +05:30
iGr3P = iGr3N
iGr3P ( faceID ( 1 ) ) = iGr3N ( faceID ( 1 ) ) + 1 ! identifying the grain ID in local coordinate sytem
iGrP = grain3to1 ( iGr3P , param ( instance ) % Nconstituents ) ! translate into global grain ID
intFaceP = getInterface ( 2 * faceID ( 1 ) - 1 , iGr3P ) ! identifying the connecting interface in local coordinate system
normP = interfaceNormal ( intFaceP , instance , of )
do iFace = 1 , 6
intFaceP = getInterface ( iFace , iGr3P ) ! identifying all interfaces that influence relaxation of the above interface
mornP = interfaceNormal ( intFaceP , instance , of )
iMun = interface4to1 ( intFaceP , param ( instance ) % Nconstituents ) ! 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
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
2020-03-17 02:09:53 +05:30
2019-01-13 03:37:35 +05:30
#ifdef DEBUG
2019-06-15 21:57:38 +05:30
if ( iand ( debug_level ( debug_homogenization ) , debug_levelExtensive ) / = 0 ) then
write ( 6 , '(1x,a30)' ) 'Jacobian matrix of stress'
do i = 1 , 3 * nIntFaceTot
write ( 6 , '(1x,100(e11.4,1x))' ) ( smatrix ( i , j ) , j = 1 , 3 * nIntFaceTot )
enddo
write ( 6 , * ) ' '
flush ( 6 )
endif
2019-01-13 03:37:35 +05:30
#endif
2020-03-17 02:09:53 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-17 02:09:53 +05:30
! ... of the stress penalty tangent (mismatch penalty and volume penalty, computed using numerical
2013-02-21 03:36:15 +05:30
! perturbation method) "pmatrix"
2019-06-15 21:57:38 +05:30
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 , 3 * nIntFaceTot
p_relax = relax
2020-03-17 02:09:53 +05:30
p_relax ( ipert ) = relax ( ipert ) + num % pPert ! perturb the relaxation vector
2019-06-15 21:57:38 +05:30
stt % relaxationVector ( : , of ) = p_relax
call grainDeformation ( pF , avgF , instance , of ) ! rain deformation from perturbed state
call stressPenalty ( pR , DevNull , avgF , pF , ip , el , instance , of ) ! stress penalty due to interface mismatch from perturbed state
call volumePenalty ( pD , devNull ( 1 , 1 ) , avgF , pF , nGrain , instance , of ) ! stress penalty due to volume discrepancy from perturbed state
2010-11-26 17:20:20 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! computing the global stress residual array from the perturbed state
2019-06-15 21:57:38 +05:30
p_resid = 0.0_pReal
do iNum = 1 , nIntFaceTot
faceID = interface1to4 ( iNum , param ( instance ) % Nconstituents ) ! identifying the interface ID in local coordinate system (4-dimensional index)
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! identify the left/bottom/back grain (-|N)
2019-06-15 21:57:38 +05:30
iGr3N = faceID ( 2 : 4 ) ! identify the grain ID in local coordinate system (3-dimensional index)
iGrN = grain3to1 ( iGr3N , param ( instance ) % Nconstituents ) ! 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 , instance , of )
2020-03-17 02:09:53 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! identify the right/up/front grain (+|P)
2019-06-15 21:57:38 +05:30
iGr3P = iGr3N
iGr3P ( faceID ( 1 ) ) = iGr3N ( faceID ( 1 ) ) + 1 ! identify the grain ID in local coordinate system (3-dimensional index)
iGrP = grain3to1 ( iGr3P , param ( instance ) % Nconstituents ) ! 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 , instance , of )
2020-03-17 02:09:53 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-17 02:09:53 +05:30
! compute the residual stress (contribution of mismatch and volume penalties) from perturbed state
2013-02-21 03:36:15 +05:30
! at all interfaces
2019-06-15 21:57:38 +05:30
do i = 1 , 3 ; do j = 1 , 3
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
2020-03-17 02:09:53 +05:30
pmatrix ( : , ipert ) = p_resid / num % pPert
2019-06-15 21:57:38 +05:30
enddo
2020-03-17 02:09:53 +05:30
2019-01-13 03:37:35 +05:30
#ifdef DEBUG
2019-06-15 21:57:38 +05:30
if ( iand ( debug_level ( debug_homogenization ) , debug_levelExtensive ) / = 0 ) then
write ( 6 , '(1x,a30)' ) 'Jacobian matrix of penalty'
do i = 1 , 3 * nIntFaceTot
write ( 6 , '(1x,100(e11.4,1x))' ) ( pmatrix ( i , j ) , j = 1 , 3 * nIntFaceTot )
enddo
write ( 6 , * ) ' '
flush ( 6 )
endif
2019-01-13 03:37:35 +05:30
#endif
2020-03-17 02:09:53 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! ... of the numerical viscosity traction "rmatrix"
2019-06-15 21:57:38 +05:30
allocate ( rmatrix ( 3 * nIntFaceTot , 3 * nIntFaceTot ) , source = 0.0_pReal )
do i = 1 , 3 * nIntFaceTot
2020-03-17 02:09:53 +05:30
rmatrix ( i , i ) = num % viscModus * num % viscPower / ( num % refRelaxRate * dt ) * & ! tangent due to numerical viscosity traction appears
( abs ( drelax ( i ) ) / ( num % refRelaxRate * dt ) ) ** ( num % viscPower - 1.0_pReal ) ! only in the main diagonal term
2019-06-15 21:57:38 +05:30
enddo
2020-03-17 02:09:53 +05:30
2019-01-13 03:37:35 +05:30
#ifdef DEBUG
2019-06-15 21:57:38 +05:30
if ( iand ( debug_level ( debug_homogenization ) , debug_levelExtensive ) / = 0 ) then
write ( 6 , '(1x,a30)' ) 'Jacobian matrix of penalty'
do i = 1 , 3 * nIntFaceTot
write ( 6 , '(1x,100(e11.4,1x))' ) ( rmatrix ( i , j ) , j = 1 , 3 * nIntFaceTot )
enddo
write ( 6 , * ) ' '
flush ( 6 )
endif
2019-01-13 03:37:35 +05:30
#endif
2009-11-17 19:12:38 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! The overall Jacobian matrix summarizing contributions of smatrix, pmatrix, rmatrix
2019-06-15 21:57:38 +05:30
allocate ( jmatrix ( 3 * nIntFaceTot , 3 * nIntFaceTot ) ) ; jmatrix = smatrix + pmatrix + rmatrix
2019-01-13 03:37:35 +05:30
#ifdef DEBUG
2019-06-15 21:57:38 +05:30
if ( iand ( debug_level ( debug_homogenization ) , debug_levelExtensive ) / = 0 ) then
write ( 6 , '(1x,a30)' ) 'Jacobian matrix (total)'
do i = 1 , 3 * nIntFaceTot
write ( 6 , '(1x,100(e11.4,1x))' ) ( jmatrix ( i , j ) , j = 1 , 3 * nIntFaceTot )
enddo
write ( 6 , * ) ' '
flush ( 6 )
endif
2019-01-13 03:37:35 +05:30
#endif
2010-11-26 17:20:20 +05:30
2013-12-18 12:58:01 +05:30
!--------------------------------------------------------------------------------------------------
2013-02-21 03:36:15 +05:30
! computing the update of the state variable (relaxation vectors) using the Jacobian matrix
2019-06-15 21:57:38 +05:30
allocate ( jnverse ( 3 * nIntFaceTot , 3 * nIntFaceTot ) , source = 0.0_pReal )
2019-09-21 06:50:33 +05:30
call math_invert ( jnverse , error , jmatrix )
2020-03-17 02:09:53 +05:30
2019-01-13 03:37:35 +05:30
#ifdef DEBUG
2019-06-15 21:57:38 +05:30
if ( iand ( debug_level ( debug_homogenization ) , debug_levelExtensive ) / = 0 ) then
write ( 6 , '(1x,a30)' ) 'Jacobian inverse'
do i = 1 , 3 * nIntFaceTot
write ( 6 , '(1x,100(e11.4,1x))' ) ( jnverse ( i , j ) , j = 1 , 3 * nIntFaceTot )
enddo
write ( 6 , * ) ' '
flush ( 6 )
endif
2019-01-13 03:37:35 +05:30
#endif
2009-10-12 21:31:49 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! calculate the state update (global relaxation vectors) for the next Newton-Raphson iteration
2019-06-15 21:57:38 +05:30
drelax = 0.0_pReal
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 ( : , of ) = relax + drelax ! Updateing the state variable for the next iteration
2020-03-17 02:09:53 +05:30
if ( any ( abs ( drelax ) > num % maxdRelax ) ) then ! Forcing cutback when the incremental change of relaxation vector becomes too large
2019-06-15 21:57:38 +05:30
mech_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
2010-11-26 17:20:20 +05:30
2019-01-13 03:37:35 +05:30
#ifdef DEBUG
2019-06-15 21:57:38 +05:30
if ( iand ( debug_homogenization , debug_levelExtensive ) > 0 ) then
write ( 6 , '(1x,a30)' ) 'Returned state: '
do i = 1 , size ( stt % relaxationVector ( : , of ) )
write ( 6 , '(1x,2(e15.8,1x))' ) stt % relaxationVector ( i , of )
enddo
write ( 6 , * ) ' '
flush ( 6 )
endif
2019-01-13 03:37:35 +05:30
#endif
2019-06-15 21:57:38 +05:30
end associate
2018-11-04 01:41:43 +05:30
2019-06-15 21:57:38 +05:30
contains
!------------------------------------------------------------------------------------------------
!> @brief calculate stress-like penalty due to deformation mismatch
!------------------------------------------------------------------------------------------------
subroutine stressPenalty ( rPen , nMis , avgF , fDef , ip , el , instance , of )
2019-05-15 02:42:32 +05:30
2019-06-15 21:57:38 +05:30
real ( pReal ) , dimension ( : , : , : ) , intent ( out ) :: rPen !< stress-like penalty
real ( pReal ) , dimension ( : , : ) , intent ( out ) :: nMis !< total amount of mismatch
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
real ( pReal ) , dimension ( : , : , : ) , intent ( in ) :: fDef !< deformation gradients
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: avgF !< initial effective stretch tensor
integer , intent ( in ) :: ip , el , instance , of
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
integer , dimension ( 4 ) :: intFace
integer , dimension ( 3 ) :: iGrain3 , iGNghb3 , nGDim
real ( pReal ) , dimension ( 3 , 3 ) :: gDef , nDef
real ( pReal ) , dimension ( 3 ) :: nVect , surfCorr
real ( pReal ) , dimension ( 2 ) :: Gmoduli
integer :: iGrain , iGNghb , iFace , i , j , k , l
real ( pReal ) :: muGrain , muGNghb , nDefNorm , bgGrain , bgGNghb
real ( pReal ) , parameter :: nDefToler = 1.0e-10_pReal
2019-01-13 03:37:35 +05:30
#ifdef DEBUG
2019-06-15 21:57:38 +05:30
logical :: debugActive
2019-01-13 03:37:35 +05:30
#endif
2018-11-04 03:53:52 +05:30
2019-06-15 21:57:38 +05:30
nGDim = param ( instance ) % Nconstituents
rPen = 0.0_pReal
nMis = 0.0_pReal
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
!----------------------------------------------------------------------------------------------
2020-03-17 02:09:53 +05:30
! get the correction factor the modulus of penalty stress representing the evolution of area of
2019-06-15 21:57:38 +05:30
! the interfaces due to deformations
2018-11-04 03:43:20 +05:30
2019-06-15 21:57:38 +05:30
surfCorr = surfaceCorrection ( avgF , instance , of )
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
associate ( prm = > param ( instance ) )
2018-11-04 03:53:52 +05:30
2019-01-13 03:37:35 +05:30
#ifdef DEBUG
2019-06-15 21:57:38 +05:30
debugActive = iand ( debug_level ( debug_homogenization ) , debug_levelExtensive ) / = 0 . and . prm % of_debug == of
2019-01-13 03:37:35 +05:30
2019-06-15 21:57:38 +05:30
if ( debugActive ) then
write ( 6 , '(1x,a20,2(1x,i3))' ) 'Correction factor: ' , ip , el
write ( 6 , * ) surfCorr
endif
2019-01-13 03:37:35 +05:30
#endif
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
!-----------------------------------------------------------------------------------------------
2020-03-17 02:09:53 +05:30
! computing the mismatch and penalty stress tensor of all grains
2019-06-15 21:57:38 +05:30
grainLoop : do iGrain = 1 , product ( prm % Nconstituents )
Gmoduli = equivalentModuli ( iGrain , ip , el )
muGrain = Gmoduli ( 1 ) ! collecting the equivalent shear modulus of grain
bgGrain = Gmoduli ( 2 ) ! and the lengthh of Burgers vector
iGrain3 = grain1to3 ( iGrain , prm % Nconstituents ) ! get the grain ID in local 3-dimensional index (x,y,z)-position
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
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 , instance , of )
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 ) )
where ( iGNghb3 < 1 ) iGNghb3 = nGDim
where ( iGNghb3 > nGDim ) iGNghb3 = 1
iGNghb = grain3to1 ( iGNghb3 , prm % Nconstituents ) ! get the ID of the neighboring grain
Gmoduli = equivalentModuli ( iGNghb , ip , el ) ! collect 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 ) ) ! difference/jump in deformation gradeint across the neighbor
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
!-------------------------------------------------------------------------------------------
! compute the mismatch tensor of all interfaces
nDefNorm = 0.0_pReal
nDef = 0.0_pReal
do i = 1 , 3 ; do j = 1 , 3
do k = 1 , 3 ; do l = 1 , 3
2020-03-02 03:30:06 +05:30
nDef ( i , j ) = nDef ( i , j ) - nVect ( k ) * gDef ( i , l ) * math_LeviCivita ( j , k , l ) ! compute the interface mismatch tensor from the jump of deformation gradient
2019-06-15 21:57:38 +05:30
enddo ; enddo
nDefNorm = nDefNorm + nDef ( i , j ) ** 2.0_pReal ! 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)
2019-01-13 03:37:35 +05:30
#ifdef DEBUG
2019-06-15 21:57:38 +05:30
if ( debugActive ) then
write ( 6 , '(1x,a20,i2,1x,a20,1x,i3)' ) 'Mismatch to face: ' , intFace ( 1 ) , 'neighbor grain: ' , iGNghb
write ( 6 , * ) transpose ( nDef )
write ( 6 , '(1x,a20,e11.4)' ) 'with magnitude: ' , nDefNorm
endif
2019-01-13 03:37:35 +05:30
#endif
2019-06-15 21:57:38 +05:30
!-------------------------------------------------------------------------------------------
! compute the stress penalty of all interfaces
do i = 1 , 3 ; do j = 1 , 3 ; do k = 1 , 3 ; do l = 1 , 3
rPen ( i , j , iGrain ) = rPen ( i , j , iGrain ) + 0.5_pReal * ( muGrain * bgGrain + muGNghb * bgGNghb ) * prm % xiAlpha &
* surfCorr ( abs ( intFace ( 1 ) ) ) / prm % dAlpha ( abs ( intFace ( 1 ) ) ) &
* cosh ( prm % ciAlpha * nDefNorm ) &
2020-03-02 03:30:06 +05:30
* 0.5_pReal * nVect ( l ) * nDef ( i , k ) / nDefNorm * math_LeviCivita ( k , l , j ) &
2020-03-17 02:09:53 +05:30
* tanh ( nDefNorm / num % xSmoo )
2019-06-15 21:57:38 +05:30
enddo ; enddo ; enddo ; enddo
enddo interfaceLoop
2019-01-13 03:37:35 +05:30
#ifdef DEBUG
2019-06-15 21:57:38 +05:30
if ( debugActive ) then
write ( 6 , '(1x,a20,i2)' ) 'Penalty of grain: ' , iGrain
write ( 6 , * ) transpose ( rPen ( 1 : 3 , 1 : 3 , iGrain ) )
endif
2019-01-13 03:37:35 +05:30
#endif
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
enddo grainLoop
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
end associate
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
end subroutine stressPenalty
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
!------------------------------------------------------------------------------------------------
2020-03-17 02:09:53 +05:30
!> @brief calculate stress-like penalty due to volume discrepancy
2019-06-15 21:57:38 +05:30
!------------------------------------------------------------------------------------------------
subroutine volumePenalty ( vPen , vDiscrep , fAvg , fDef , nGrain , instance , of )
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
real ( pReal ) , dimension ( : , : , : ) , intent ( out ) :: vPen ! stress-like penalty due to volume
real ( pReal ) , intent ( out ) :: vDiscrep ! total volume discrepancy
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
real ( pReal ) , dimension ( : , : , : ) , intent ( in ) :: fDef ! deformation gradients
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: fAvg ! overall deformation gradient
integer , intent ( in ) :: &
Ngrain , &
instance , &
of
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
real ( pReal ) , dimension ( size ( vPen , 3 ) ) :: gVol
integer :: i
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
!----------------------------------------------------------------------------------------------
! compute the volumes of grains and of cluster
vDiscrep = math_det33 ( fAvg ) ! compute the volume of the cluster
do i = 1 , nGrain
gVol ( i ) = math_det33 ( fDef ( 1 : 3 , 1 : 3 , i ) ) ! compute the volume of individual grains
vDiscrep = vDiscrep - gVol ( i ) / real ( nGrain , pReal ) ! calculate the difference/dicrepancy between
2019-01-13 03:37:35 +05:30
! the volume of the cluster and the the total volume of grains
2019-06-15 21:57:38 +05:30
enddo
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
!----------------------------------------------------------------------------------------------
! calculate the stress and penalty due to volume discrepancy
vPen = 0.0_pReal
do i = 1 , nGrain
2020-03-17 02:09:53 +05:30
vPen ( : , : , i ) = - 1.0_pReal / real ( nGrain , pReal ) * num % volDiscrMod * num % volDiscrPow / num % maxVolDiscr * &
sign ( ( abs ( vDiscrep ) / num % maxVolDiscr ) ** ( num % volDiscrPow - 1.0 ) , vDiscrep ) * &
2019-06-15 21:57:38 +05:30
gVol ( i ) * transpose ( math_inv33 ( fDef ( : , : , i ) ) )
2020-03-17 02:09:53 +05:30
2019-01-13 03:37:35 +05:30
#ifdef DEBUG
2019-06-15 21:57:38 +05:30
if ( iand ( debug_level ( debug_homogenization ) , debug_levelExtensive ) / = 0 &
. and . param ( instance ) % of_debug == of ) then
write ( 6 , '(1x,a30,i2)' ) 'Volume penalty of grain: ' , i
write ( 6 , * ) transpose ( vPen ( : , : , i ) )
endif
2019-01-13 03:37:35 +05:30
#endif
2019-06-15 21:57:38 +05:30
enddo
2018-11-04 02:06:49 +05:30
2019-06-15 21:57:38 +05:30
end subroutine volumePenalty
2018-11-04 02:06:49 +05:30
2019-06-15 21:57:38 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-17 02:09:53 +05:30
!> @brief compute the correction factor accouted for surface evolution (area change) due to
2019-06-15 21:57:38 +05:30
! deformation
!--------------------------------------------------------------------------------------------------
function surfaceCorrection ( avgF , instance , of )
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
real ( pReal ) , dimension ( 3 ) :: surfaceCorrection
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: avgF !< average F
integer , intent ( in ) :: &
instance , &
of
real ( pReal ) , dimension ( 3 , 3 ) :: invC
real ( pReal ) , dimension ( 3 ) :: nVect
real ( pReal ) :: detF
integer :: i , j , iBase
logical :: error
2020-03-17 02:09:53 +05:30
2019-09-21 06:58:46 +05:30
call math_invert33 ( invC , detF , error , matmul ( transpose ( avgF ) , avgF ) )
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
surfaceCorrection = 0.0_pReal
do iBase = 1 , 3
nVect = interfaceNormal ( [ iBase , 1 , 1 , 1 ] , instance , of )
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
surfaceCorrection ( iBase ) = sqrt ( surfaceCorrection ( iBase ) ) * detF ! get the surface correction factor (area contraction/enlargement)
enddo
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
end function surfaceCorrection
2018-11-04 02:06:49 +05:30
2019-06-15 21:57:38 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief compute the equivalent shear and bulk moduli from the elasticity tensor
!--------------------------------------------------------------------------------------------------
function equivalentModuli ( grainID , ip , el )
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
real ( pReal ) , dimension ( 2 ) :: equivalentModuli
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
integer , intent ( in ) :: &
grainID , &
ip , & !< integration point number
el !< element number
real ( pReal ) , dimension ( 6 , 6 ) :: elasTens
real ( pReal ) :: &
cEquiv_11 , &
cEquiv_12 , &
cEquiv_44
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
elasTens = constitutive_homogenizedC ( grainID , ip , el )
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
!----------------------------------------------------------------------------------------------
! 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
equivalentModuli ( 1 ) = 0.2_pReal * ( cEquiv_11 - cEquiv_12 ) + 0.6_pReal * cEquiv_44
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
!----------------------------------------------------------------------------------------------
! obtain the length of Burgers vector (could be model dependend)
equivalentModuli ( 2 ) = 2.5e-10_pReal
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
end function equivalentModuli
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-17 02:09:53 +05:30
!> @brief calculating the grain deformation gradient (the same with
2019-06-15 21:57:38 +05:30
! homogenization_RGC_partitionDeformation, but used only for perturbation scheme)
!--------------------------------------------------------------------------------------------------
subroutine grainDeformation ( F , avgF , instance , of )
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
real ( pReal ) , dimension ( : , : , : ) , intent ( out ) :: F !< partioned F per grain
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
real ( pReal ) , dimension ( : , : ) , intent ( in ) :: avgF !< averaged F
integer , intent ( in ) :: &
instance , &
of
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
real ( pReal ) , dimension ( 3 ) :: aVect , nVect
integer , dimension ( 4 ) :: intFace
integer , dimension ( 3 ) :: iGrain3
integer :: iGrain , iFace , i , j
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
!-------------------------------------------------------------------------------------------------
! compute the deformation gradient of individual grains due to relaxations
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
associate ( prm = > param ( instance ) )
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
F = 0.0_pReal
do iGrain = 1 , product ( prm % Nconstituents )
iGrain3 = grain1to3 ( iGrain , prm % Nconstituents )
do iFace = 1 , 6
intFace = getInterface ( iFace , iGrain3 )
aVect = relaxationVector ( intFace , instance , of )
nVect = interfaceNormal ( intFace , instance , of )
forall ( i = 1 : 3 , j = 1 : 3 ) &
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
2018-11-04 11:57:25 +05:30
enddo
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
end associate
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
end subroutine grainDeformation
2020-03-17 02:09:53 +05:30
2019-05-18 13:17:20 +05:30
end procedure mech_RGC_updateState
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-17 02:09:53 +05:30
!> @brief derive average stress and stiffness from constituent quantities
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
2019-05-18 10:53:46 +05:30
module subroutine mech_RGC_averageStressAndItsTangent ( avgP , dAvgPdAvgF , P , dPdF , instance )
2018-11-04 00:38:18 +05:30
2019-06-15 21:57:38 +05:30
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
2019-01-13 13:44:23 +05:30
2019-06-15 21:57:38 +05:30
real ( pReal ) , dimension ( : , : , : ) , intent ( in ) :: P !< partitioned stresses
real ( pReal ) , dimension ( : , : , : , : , : ) , intent ( in ) :: dPdF !< partitioned stiffnesses
integer , intent ( in ) :: instance
2009-10-12 21:31:49 +05:30
2019-06-15 21:57:38 +05:30
avgP = sum ( P , 3 ) / real ( product ( param ( instance ) % Nconstituents ) , pReal )
dAvgPdAvgF = sum ( dPdF , 5 ) / real ( product ( param ( instance ) % Nconstituents ) , pReal )
2009-10-12 21:31:49 +05:30
2019-05-18 10:53:46 +05:30
end subroutine mech_RGC_averageStressAndItsTangent
2009-10-12 21:31:49 +05:30
2019-04-30 22:15:16 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief writes results to HDF5 output file
!--------------------------------------------------------------------------------------------------
2019-05-18 11:09:55 +05:30
module subroutine mech_RGC_results ( instance , group )
2019-04-30 22:15:16 +05:30
2019-05-18 10:53:46 +05:30
integer , intent ( in ) :: instance
character ( len = * ) , intent ( in ) :: group
2020-03-17 02:09:53 +05:30
2019-04-30 22:15:16 +05:30
integer :: o
2020-03-17 02:09:53 +05:30
2019-04-30 22:15:16 +05:30
associate ( stt = > state ( instance ) , dst = > dependentState ( instance ) , prm = > param ( instance ) )
2020-02-15 03:51:58 +05:30
outputsLoop : do o = 1 , size ( prm % output )
select case ( trim ( prm % output ( o ) ) )
case ( 'constitutivework' )
2019-05-01 02:17:49 +05:30
call results_writeDataset ( group , stt % work , 'W' , &
'work density' , 'J/m³' )
2020-02-15 03:51:58 +05:30
case ( 'magnitudemismatch' )
2019-05-12 18:44:29 +05:30
call results_writeDataset ( group , dst % mismatch , 'N' , &
2019-05-13 21:03:45 +05:30
'average mismatch tensor' , '1' )
2020-02-15 03:51:58 +05:30
case ( 'penaltyenergy' )
2019-05-01 02:17:49 +05:30
call results_writeDataset ( group , stt % penaltyEnergy , 'R' , &
2019-05-12 18:44:29 +05:30
'mismatch penalty density' , 'J/m³' )
2020-02-15 03:51:58 +05:30
case ( 'volumediscrepancy' )
2019-05-12 18:44:29 +05:30
call results_writeDataset ( group , dst % volumeDiscrepancy , 'Delta_V' , &
'volume discrepancy' , 'm³' )
2020-02-15 03:51:58 +05:30
case ( 'maximumrelaxrate' )
2019-05-12 18:44:29 +05:30
call results_writeDataset ( group , dst % relaxationrate_max , 'max_alpha_dot' , &
'maximum relaxation rate' , 'm/s' )
2020-02-15 03:51:58 +05:30
case ( 'averagerelaxrate' )
2019-05-12 18:44:29 +05:30
call results_writeDataset ( group , dst % relaxationrate_avg , 'avg_alpha_dot' , &
'average relaxation rate' , 'm/s' )
2019-04-30 22:15:16 +05:30
end select
enddo outputsLoop
end associate
2020-03-17 02:09:53 +05:30
2019-04-30 22:15:16 +05:30
end subroutine mech_RGC_results
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief collect relaxation vectors of an interface
!--------------------------------------------------------------------------------------------------
2018-11-04 12:03:57 +05:30
pure function relaxationVector ( intFace , instance , of )
2013-02-21 03:36:15 +05:30
2019-06-15 21:57:38 +05:30
real ( pReal ) , dimension ( 3 ) :: relaxationVector
2019-04-06 00:15:56 +05:30
2019-06-15 21:57:38 +05:30
integer , intent ( in ) :: instance , of
integer , dimension ( 4 ) , intent ( in ) :: intFace !< set of interface ID in 4D array (normal and position)
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
integer :: iNum
2009-10-12 21:31:49 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! collect the interface relaxation vector from the global state array
2019-01-13 03:37:35 +05:30
2019-06-15 21:57:38 +05:30
iNum = interface4to1 ( intFace , param ( instance ) % Nconstituents ) ! identify the position of the interface in global state array
if ( iNum > 0 ) then
relaxationVector = state ( instance ) % relaxationVector ( ( 3 * iNum - 2 ) : ( 3 * iNum ) , of )
else
relaxationVector = 0.0_pReal
endif
2013-02-21 03:36:15 +05:30
2018-08-30 14:20:52 +05:30
end function relaxationVector
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
2020-03-17 02:09:53 +05:30
!> @brief identify the normal of an interface
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
2018-11-04 12:03:57 +05:30
pure function interfaceNormal ( intFace , instance , of )
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
real ( pReal ) , dimension ( 3 ) :: interfaceNormal
2019-04-06 00:15:56 +05:30
2019-06-15 21:57:38 +05:30
integer , dimension ( 4 ) , intent ( in ) :: intFace !< interface ID in 4D array (normal and position)
integer , intent ( in ) :: &
instance , &
of
2019-04-06 00:15:56 +05:30
2019-06-15 21:57:38 +05:30
integer :: nPos
2009-10-12 21:31:49 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! get the normal of the interface, identified from the value of intFace(1)
2019-06-15 21:57:38 +05:30
interfaceNormal = 0.0_pReal
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
2009-10-12 21:31:49 +05:30
2019-06-15 21:57:38 +05:30
interfaceNormal = matmul ( dependentState ( instance ) % orientation ( 1 : 3 , 1 : 3 , of ) , interfaceNormal ) ! map the normal vector into sample coordinate system (basis)
2010-03-24 18:50:12 +05:30
2018-08-30 14:20:52 +05:30
end function interfaceNormal
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief collect six faces of a grain in 4D (normal and position)
!--------------------------------------------------------------------------------------------------
2018-11-04 12:03:57 +05:30
pure function getInterface ( iFace , iGrain3 )
2013-02-21 03:36:15 +05:30
2019-06-15 21:57:38 +05:30
integer , dimension ( 4 ) :: getInterface
2019-04-06 00:15:56 +05:30
2019-06-15 21:57:38 +05:30
integer , dimension ( 3 ) , intent ( in ) :: iGrain3 !< grain ID in 3D array
integer , intent ( in ) :: iFace !< face index (1..6) mapped like (-e1,-e2,-e3,+e1,+e2,+e3) or iDir = (-1,-2,-3,1,2,3)
2019-04-06 00:15:56 +05:30
2019-06-15 21:57:38 +05:30
integer :: iDir !< direction of interface normal
2020-03-17 02:09:53 +05:30
2019-04-06 00:15:56 +05:30
iDir = ( int ( real ( iFace - 1 , pReal ) / 2.0_pReal ) + 1 ) * ( - 1 ) ** iFace
2018-08-30 14:20:52 +05:30
getInterface ( 1 ) = iDir
2020-03-17 02:09:53 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! identify the interface position by the direction of its normal
2019-06-15 21:57:38 +05:30
getInterface ( 2 : 4 ) = iGrain3
if ( iDir < 0 ) getInterface ( 1 - iDir ) = getInterface ( 1 - iDir ) - 1 ! to have a correlation with coordinate/position in real space
2009-10-12 21:31:49 +05:30
2018-08-30 14:20:52 +05:30
end function getInterface
2009-10-12 21:31:49 +05:30
2018-11-04 12:41:35 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief map grain ID from in 1D (global array) to in 3D (local position)
!--------------------------------------------------------------------------------------------------
2018-11-04 12:41:35 +05:30
pure function grain1to3 ( grain1 , nGDim )
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
integer , dimension ( 3 ) :: grain1to3
2019-04-06 00:15:56 +05:30
2019-06-15 21:57:38 +05:30
integer , intent ( in ) :: grain1 !< grain ID in 1D array
integer , dimension ( 3 ) , intent ( in ) :: nGDim
2009-10-12 21:31:49 +05:30
2019-06-15 21:57:38 +05:30
grain1to3 = 1 + [ mod ( ( grain1 - 1 ) , nGDim ( 1 ) ) , &
mod ( ( grain1 - 1 ) / nGDim ( 1 ) , nGDim ( 2 ) ) , &
2019-04-06 00:15:56 +05:30
( grain1 - 1 ) / ( nGDim ( 1 ) * nGDim ( 2 ) ) ]
2009-10-12 21:31:49 +05:30
2018-08-30 14:20:52 +05:30
end function grain1to3
2009-10-12 21:31:49 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief map grain ID from in 3D (local position) to in 1D (global array)
!--------------------------------------------------------------------------------------------------
2019-04-06 00:15:56 +05:30
integer pure function grain3to1 ( grain3 , nGDim )
2009-10-12 21:31:49 +05:30
2019-06-15 21:57:38 +05:30
integer , dimension ( 3 ) , intent ( in ) :: grain3 !< grain ID in 3D array (pos.x,pos.y,pos.z)
integer , dimension ( 3 ) , intent ( in ) :: nGDim
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
grain3to1 = grain3 ( 1 ) &
+ nGDim ( 1 ) * ( grain3 ( 2 ) - 1 ) &
+ nGDim ( 1 ) * nGDim ( 2 ) * ( grain3 ( 3 ) - 1 )
2009-10-12 21:31:49 +05:30
2018-08-30 14:20:52 +05:30
end function grain3to1
2009-10-12 21:31:49 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief maps interface ID from 4D (normal and local position) into 1D (global array)
!--------------------------------------------------------------------------------------------------
2019-04-06 00:15:56 +05:30
integer pure function interface4to1 ( iFace4D , nGDim )
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
integer , dimension ( 4 ) , intent ( in ) :: iFace4D !< interface ID in 4D array (n.dir,pos.x,pos.y,pos.z)
integer , dimension ( 3 ) , intent ( in ) :: nGDim
2010-11-26 17:20:20 +05:30
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
select case ( abs ( iFace4D ( 1 ) ) )
case ( 1 )
if ( ( iFace4D ( 2 ) == 0 ) . or . ( iFace4D ( 2 ) == nGDim ( 1 ) ) ) then
interface4to1 = 0
else
interface4to1 = iFace4D ( 3 ) + nGDim ( 2 ) * ( iFace4D ( 4 ) - 1 ) &
+ nGDim ( 2 ) * nGDim ( 3 ) * ( iFace4D ( 2 ) - 1 )
endif
2018-11-04 13:49:24 +05:30
2019-06-15 21:57:38 +05:30
case ( 2 )
if ( ( iFace4D ( 3 ) == 0 ) . or . ( iFace4D ( 3 ) == nGDim ( 2 ) ) ) then
interface4to1 = 0
else
interface4to1 = iFace4D ( 4 ) + nGDim ( 3 ) * ( iFace4D ( 2 ) - 1 ) &
+ nGDim ( 3 ) * nGDim ( 1 ) * ( iFace4D ( 3 ) - 1 ) &
+ ( nGDim ( 1 ) - 1 ) * nGDim ( 2 ) * nGDim ( 3 ) ! total # of interfaces normal || e1
endif
2018-11-04 13:49:24 +05:30
2019-06-15 21:57:38 +05:30
case ( 3 )
if ( ( iFace4D ( 4 ) == 0 ) . or . ( iFace4D ( 4 ) == nGDim ( 3 ) ) ) then
interface4to1 = 0
else
interface4to1 = iFace4D ( 2 ) + nGDim ( 1 ) * ( iFace4D ( 3 ) - 1 ) &
+ nGDim ( 1 ) * nGDim ( 2 ) * ( iFace4D ( 4 ) - 1 ) &
+ ( nGDim ( 1 ) - 1 ) * nGDim ( 2 ) * nGDim ( 3 ) & ! total # of interfaces normal || e1
+ nGDim ( 1 ) * ( nGDim ( 2 ) - 1 ) * nGDim ( 3 ) ! total # of interfaces normal || e2
endif
2018-11-04 13:49:24 +05:30
2019-06-15 21:57:38 +05:30
case default
interface4to1 = - 1
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
end select
2009-10-12 21:31:49 +05:30
2018-08-30 14:20:52 +05:30
end function interface4to1
2009-10-12 21:31:49 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief maps interface ID from 1D (global array) into 4D (normal and local position)
!--------------------------------------------------------------------------------------------------
2018-11-04 12:41:35 +05:30
pure function interface1to4 ( iFace1D , nGDim )
2020-03-17 02:09:53 +05:30
2019-06-15 21:57:38 +05:30
integer , dimension ( 4 ) :: interface1to4
2019-04-06 00:15:56 +05:30
2019-06-15 21:57:38 +05:30
integer , intent ( in ) :: iFace1D !< interface ID in 1D array
integer , dimension ( 3 ) , intent ( in ) :: nGDim
integer , dimension ( 3 ) :: nIntFace
2009-10-12 21:31:49 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! compute the total number of interfaces, which ...
2019-06-15 21:57:38 +05:30
nIntFace = [ ( nGDim ( 1 ) - 1 ) * nGDim ( 2 ) * nGDim ( 3 ) , & ! ... normal || e1
nGDim ( 1 ) * ( nGDim ( 2 ) - 1 ) * nGDim ( 3 ) , & ! ... normal || e2
nGDim ( 1 ) * nGDim ( 2 ) * ( nGDim ( 3 ) - 1 ) ] ! ... normal || e3
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! get the corresponding interface ID in 4D (normal and local position)
2019-06-15 21:57:38 +05:30
if ( iFace1D > 0 . and . iFace1D < = nIntFace ( 1 ) ) then ! interface with normal || e1
interface1to4 ( 1 ) = 1
interface1to4 ( 3 ) = mod ( ( iFace1D - 1 ) , nGDim ( 2 ) ) + 1
interface1to4 ( 4 ) = mod ( int ( real ( iFace1D - 1 , pReal ) / real ( nGDim ( 2 ) , pReal ) ) , nGDim ( 3 ) ) + 1
interface1to4 ( 2 ) = int ( real ( iFace1D - 1 , pReal ) / real ( nGDim ( 2 ) , pReal ) / real ( nGDim ( 3 ) , pReal ) ) + 1
elseif ( iFace1D > nIntFace ( 1 ) . and . iFace1D < = ( nIntFace ( 2 ) + nIntFace ( 1 ) ) ) then ! interface with normal || e2
interface1to4 ( 1 ) = 2
interface1to4 ( 4 ) = mod ( ( iFace1D - nIntFace ( 1 ) - 1 ) , nGDim ( 3 ) ) + 1
interface1to4 ( 2 ) = mod ( int ( real ( iFace1D - nIntFace ( 1 ) - 1 , pReal ) / real ( nGDim ( 3 ) , pReal ) ) , nGDim ( 1 ) ) + 1
interface1to4 ( 3 ) = int ( real ( iFace1D - nIntFace ( 1 ) - 1 , pReal ) / real ( nGDim ( 3 ) , pReal ) / real ( nGDim ( 1 ) , pReal ) ) + 1
elseif ( iFace1D > nIntFace ( 2 ) + nIntFace ( 1 ) . and . iFace1D < = ( nIntFace ( 3 ) + nIntFace ( 2 ) + nIntFace ( 1 ) ) ) then ! interface with normal || e3
interface1to4 ( 1 ) = 3
interface1to4 ( 2 ) = mod ( ( iFace1D - nIntFace ( 2 ) - nIntFace ( 1 ) - 1 ) , nGDim ( 1 ) ) + 1
interface1to4 ( 3 ) = mod ( int ( real ( iFace1D - nIntFace ( 2 ) - nIntFace ( 1 ) - 1 , pReal ) / real ( nGDim ( 1 ) , pReal ) ) , nGDim ( 2 ) ) + 1
interface1to4 ( 4 ) = int ( real ( iFace1D - nIntFace ( 2 ) - nIntFace ( 1 ) - 1 , pReal ) / real ( nGDim ( 1 ) , pReal ) / real ( nGDim ( 2 ) , pReal ) ) + 1
endif
2009-10-12 21:31:49 +05:30
2018-08-30 14:20:52 +05:30
end function interface1to4
2013-02-21 03:36:15 +05:30
2019-05-18 10:53:46 +05:30
end submodule homogenization_mech_RGC