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
!--------------------------------------------------------------------------------------------------
module homogenization_RGC
use prec , only : &
pReal , &
pInt
2009-10-12 21:31:49 +05:30
2012-03-09 01:55:28 +05:30
implicit none
2013-02-21 03:36:15 +05:30
private
2013-12-18 12:58:01 +05:30
integer ( pInt ) , dimension ( : ) , allocatable , public :: &
2013-02-21 03:36:15 +05:30
homogenization_RGC_sizePostResults
2013-12-18 12:58:01 +05:30
integer ( pInt ) , dimension ( : , : ) , allocatable , target , public :: &
2013-02-21 03:36:15 +05:30
homogenization_RGC_sizePostResult
2013-12-18 12:58:01 +05:30
character ( len = 64 ) , dimension ( : , : ) , allocatable , target , public :: &
2013-02-21 03:36:15 +05:30
homogenization_RGC_output ! name of each post result output
2014-09-26 16:04:36 +05:30
integer ( pInt ) , dimension ( : ) , allocatable , target , public :: &
2018-10-13 21:51:13 +05:30
homogenization_RGC_Noutput !< number of outputs per homog instance
2018-08-28 16:27:22 +05:30
2018-08-31 04:54:28 +05:30
enum , bind ( c )
enumerator :: undefined_ID , &
constitutivework_ID , &
penaltyenergy_ID , &
volumediscrepancy_ID , &
averagerelaxrate_ID , &
maximumrelaxrate_ID , &
ipcoords_ID , &
magnitudemismatch_ID , &
avgdefgrad_ID , &
avgfirstpiola_ID
end enum
2018-08-28 16:27:22 +05:30
type , private :: tParameters !< container type for internal constitutive parameters
integer ( pInt ) , dimension ( : ) , allocatable :: &
Nconstituents
real ( pReal ) :: &
xiAlpha , &
ciAlpha
real ( pReal ) , dimension ( : ) , allocatable :: &
dAlpha , &
angles
2018-08-31 04:54:28 +05:30
integer ( kind ( undefined_ID ) ) , dimension ( : ) , allocatable :: &
outputID !< ID of each post result output
2018-08-28 16:27:22 +05:30
end type
2018-11-04 02:06:49 +05:30
type , private :: tRGCstate
2018-11-03 21:20:43 +05:30
real ( pReal ) , pointer , dimension ( : ) :: &
2018-10-13 21:51:13 +05:30
work , &
penaltyEnergy , &
volumeDiscrepancy , &
relaxationRate_avg , &
2018-11-04 02:06:49 +05:30
relaxationRate_max
2018-11-03 21:20:43 +05:30
real ( pReal ) , pointer , dimension ( : , : ) :: &
2018-11-04 03:43:20 +05:30
relaxationVector , &
2018-11-03 21:20:43 +05:30
mismatch
2018-11-04 02:06:49 +05:30
end type tRGCstate
type , private :: tRGCdependentState
real ( pReal ) , allocatable , dimension ( : , : , : ) :: &
orientation
end type tRGCdependentState
type ( tparameters ) , dimension ( : ) , allocatable , private :: param !< containers of parameters (len Ninstance)
type ( tRGCstate ) , dimension ( : ) , allocatable , private :: state
type ( tRGCdependentState ) , dimension ( : ) , allocatable , private :: dependentState
2013-02-21 03:36:15 +05:30
public :: &
homogenization_RGC_init , &
homogenization_RGC_partitionDeformation , &
homogenization_RGC_averageStressAndItsTangent , &
homogenization_RGC_updateState , &
homogenization_RGC_postResults
private :: &
2018-08-30 14:20:52 +05:30
relaxationVector , &
interfaceNormal , &
getInterface , &
grain1to3 , &
grain3to1 , &
interface4to1 , &
interface1to4
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
!--------------------------------------------------------------------------------------------------
2018-11-03 21:20:43 +05:30
subroutine homogenization_RGC_init ( )
2018-02-02 17:06:09 +05:30
#if defined(__GFORTRAN__) || __INTEL_COMPILER >= 1800
2017-10-05 20:05:34 +05:30
use , intrinsic :: iso_fortran_env , only : &
compiler_version , &
compiler_options
#endif
2014-08-21 23:18:20 +05:30
use prec , only : &
pReal , &
pInt
2013-02-21 03:36:15 +05:30
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 : &
2018-10-13 21:51:13 +05:30
mesh_NcpElems , &
mesh_NipsPerElem , &
mesh_homogenizationAt , &
mesh_microstructureAt
2009-10-12 21:31:49 +05:30
use IO
use material
2018-06-14 10:09:49 +05:30
use config
2012-03-09 01:55:28 +05:30
implicit none
2018-11-03 21:20:43 +05:30
integer ( pInt ) :: &
2014-08-21 23:18:20 +05:30
NofMyHomog , &
2018-11-04 00:38:18 +05:30
h , &
2018-08-31 04:54:28 +05:30
outputSize , &
2014-08-21 23:18:20 +05:30
instance , &
2018-11-03 21:20:43 +05:30
sizeHState , nIntFaceTot
2018-11-04 03:13:52 +05:30
integer ( pInt ) :: maxNinstance , i , j , e , mySize , of
2018-08-31 04:54:28 +05:30
character ( len = 65536 ) , dimension ( 0 ) , parameter :: emptyStringArray = [ character ( len = 65536 ) :: ]
integer ( kind ( undefined_ID ) ) :: &
outputID !< ID of each post result output
character ( len = 65536 ) , dimension ( : ) , allocatable :: outputs
2018-08-28 16:45:38 +05:30
type ( tParameters ) :: prm
2013-01-09 03:41:59 +05:30
2017-11-21 19:38:45 +05:30
write ( 6 , '(/,a)' ) ' <<<+- homogenization_' / / HOMOGENIZATION_RGC_label / / ' init -+>>>'
2018-04-22 13:37:49 +05:30
write ( 6 , '(/,a)' ) ' Tjahjanto et al., International Journal of Material Forming, 2(1):939– 942, 2009'
2018-08-31 04:54:28 +05:30
write ( 6 , '(a)' ) ' https://doi.org/10.1007/s12289-009-0619-1'
2018-04-22 13:37:49 +05:30
write ( 6 , '(/,a)' ) ' Tjahjanto et al., Modelling and Simulation in Materials Science and Engineering, 18:015006, 2010'
2018-08-31 04:54:28 +05:30
write ( 6 , '(a)' ) ' https://doi.org/10.1088/0965-0393/18/1/015006'
2017-11-21 19:38:45 +05:30
write ( 6 , '(a15,a)' ) ' Current time: ' , IO_timeStamp ( )
2012-02-01 00:48:55 +05:30
#include "compilation_info.f90"
2009-10-12 21:31:49 +05:30
2013-11-27 13:34:05 +05:30
maxNinstance = int ( count ( homogenization_type == HOMOGENIZATION_RGC_ID ) , pInt )
2013-02-21 03:36:15 +05:30
if ( maxNinstance == 0_pInt ) return
2018-08-28 16:27:22 +05:30
if ( iand ( debug_level ( debug_HOMOGENIZATION ) , debug_levelBasic ) / = 0_pInt ) &
2014-09-10 19:44:03 +05:30
write ( 6 , '(a16,1x,i5,/)' ) '# instances:' , maxNinstance
2013-12-18 12:58:01 +05:30
allocate ( homogenization_RGC_sizePostResults ( maxNinstance ) , source = 0_pInt )
2018-08-28 16:27:22 +05:30
2018-08-29 13:22:22 +05:30
allocate ( param ( maxNinstance ) ) ! one container of parameters per instance
2018-11-03 21:20:43 +05:30
allocate ( state ( maxNinstance ) ) ! one container per instance
2018-11-04 02:06:49 +05:30
allocate ( dependentState ( maxNinstance ) ) ! one container per instance
2018-08-29 13:22:22 +05:30
2014-09-10 19:44:03 +05:30
allocate ( homogenization_RGC_Noutput ( maxNinstance ) , source = 0_pInt )
2013-02-21 03:36:15 +05:30
allocate ( homogenization_RGC_output ( maxval ( homogenization_Noutput ) , maxNinstance ) )
2013-12-18 12:58:01 +05:30
homogenization_RGC_output = ''
allocate ( homogenization_RGC_sizePostResult ( maxval ( homogenization_Noutput ) , maxNinstance ) , &
source = 0_pInt )
2018-08-29 13:22:22 +05:30
2018-08-28 16:45:38 +05:30
do h = 1_pInt , size ( homogenization_type )
if ( homogenization_type ( h ) / = HOMOGENIZATION_RGC_ID ) cycle
instance = homogenization_typeInstance ( h )
associate ( prm = > param ( instance ) )
2018-08-31 04:54:28 +05:30
2018-08-28 16:45:38 +05:30
prm % Nconstituents = config_homogenization ( h ) % getInts ( 'clustersize' , requiredShape = [ 3 ] )
2018-08-29 13:22:22 +05:30
if ( homogenization_Ngrains ( h ) / = product ( prm % Nconstituents ) ) &
2018-08-31 04:54:28 +05:30
call IO_error ( 211_pInt , ext_msg = 'clustersize (' / / HOMOGENIZATION_RGC_label / / ')' )
2018-08-28 16:45:38 +05:30
prm % xiAlpha = config_homogenization ( h ) % getFloat ( 'scalingparameter' )
prm % ciAlpha = config_homogenization ( h ) % getFloat ( 'overproportionality' )
prm % dAlpha = config_homogenization ( h ) % getFloats ( 'grainsize' , requiredShape = [ 3 ] )
2018-08-29 17:06:46 +05:30
prm % angles = config_homogenization ( h ) % getFloats ( 'clusterorientation' , requiredShape = [ 3 ] , &
defaultVal = [ 40 0.0_pReal , 40 0.0_pReal , 40 0.0_pReal ] )
2018-08-31 04:54:28 +05:30
outputs = config_homogenization ( h ) % getStrings ( '(output)' , defaultVal = emptyStringArray )
allocate ( prm % outputID ( 0 ) )
do i = 1_pInt , size ( outputs )
outputID = undefined_ID
select case ( outputs ( i ) )
case ( 'constitutivework' )
outputID = constitutivework_ID
outputSize = 1_pInt
case ( 'penaltyenergy' )
outputID = penaltyenergy_ID
outputSize = 1_pInt
case ( 'volumediscrepancy' )
outputID = volumediscrepancy_ID
outputSize = 1_pInt
case ( 'averagerelaxrate' )
outputID = averagerelaxrate_ID
outputSize = 1_pInt
case ( 'maximumrelaxrate' )
outputID = maximumrelaxrate_ID
outputSize = 1_pInt
case ( 'magnitudemismatch' )
outputID = magnitudemismatch_ID
outputSize = 3_pInt
case default
end select
2018-11-03 19:43:11 +05:30
if ( outputID / = undefined_ID ) then
homogenization_RGC_output ( i , instance ) = outputs ( i )
homogenization_RGC_sizePostResult ( i , instance ) = outputSize
prm % outputID = [ prm % outputID , outputID ]
endif
2018-08-31 04:54:28 +05:30
enddo
2018-08-29 17:06:46 +05:30
2018-08-29 17:22:26 +05:30
if ( iand ( debug_level ( debug_homogenization ) , debug_levelExtensive ) / = 0_pInt ) then
write ( 6 , '(a15,1x,i4,/)' ) 'instance: ' , instance
write ( 6 , '(a25,3(1x,i8))' ) 'cluster size: ' , ( prm % Nconstituents ( j ) , j = 1_pInt , 3_pInt )
write ( 6 , '(a25,1x,e10.3)' ) 'scaling parameter: ' , prm % xiAlpha
write ( 6 , '(a25,1x,e10.3)' ) 'over-proportionality: ' , prm % ciAlpha
write ( 6 , '(a25,3(1x,e10.3))' ) 'grain size: ' , ( prm % dAlpha ( j ) , j = 1_pInt , 3_pInt )
write ( 6 , '(a25,3(1x,e10.3))' ) 'cluster orientation: ' , ( prm % angles ( j ) , j = 1_pInt , 3_pInt )
endif
2018-08-31 10:05:30 +05:30
NofMyHomog = count ( material_homog == h )
2018-11-03 21:20:43 +05:30
nIntFaceTot = 3_pInt * ( ( prm % Nconstituents ( 1 ) - 1_pInt ) * prm % Nconstituents ( 2 ) * prm % Nconstituents ( 3 ) &
+ prm % Nconstituents ( 1 ) * ( prm % Nconstituents ( 2 ) - 1_pInt ) * prm % Nconstituents ( 3 ) &
+ prm % Nconstituents ( 1 ) * prm % Nconstituents ( 2 ) * ( prm % Nconstituents ( 3 ) - 1_pInt ) )
sizeHState = nIntFaceTot &
2018-08-31 10:05:30 +05:30
+ 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
homogState ( h ) % sizeState = sizeHState
2018-11-03 19:43:11 +05:30
homogenization_RGC_sizePostResults ( instance ) = sum ( homogenization_RGC_sizePostResult ( : , instance ) )
2018-08-31 10:05:30 +05:30
homogState ( h ) % sizePostResults = homogenization_RGC_sizePostResults ( instance )
allocate ( homogState ( h ) % state0 ( sizeHState , NofMyHomog ) , source = 0.0_pReal )
allocate ( homogState ( h ) % subState0 ( sizeHState , NofMyHomog ) , source = 0.0_pReal )
allocate ( homogState ( h ) % state ( sizeHState , NofMyHomog ) , source = 0.0_pReal )
2018-11-04 03:43:20 +05:30
state ( instance ) % relaxationVector = > homogState ( h ) % state ( 1 : nIntFaceTot , : )
state ( instance ) % work = > homogState ( h ) % state ( nIntFaceTot + 1 , : )
state ( instance ) % mismatch = > homogState ( h ) % state ( nIntFaceTot + 2 : nIntFaceTot + 4 , : )
state ( instance ) % penaltyEnergy = > homogState ( h ) % state ( nIntFaceTot + 5 , : )
state ( instance ) % volumeDiscrepancy = > homogState ( h ) % state ( nIntFaceTot + 6 , : )
state ( instance ) % relaxationRate_avg = > homogState ( h ) % state ( nIntFaceTot + 7 , : )
state ( instance ) % relaxationRate_max = > homogState ( h ) % state ( nIntFaceTot + 8 , : )
2018-11-04 02:06:49 +05:30
allocate ( dependentState ( instance ) % orientation ( 3 , 3 , NofMyHomog ) )
2018-11-04 00:30:02 +05:30
2018-11-04 02:06:49 +05:30
!--------------------------------------------------------------------------------------------------
! * assigning cluster orientations
elementLooping : do e = 1_pInt , mesh_NcpElems
2018-11-04 11:57:25 +05:30
if ( homogenization_typeInstance ( mesh_homogenizationAt ( e ) ) == instance . and . NofMyHomog > 0_pInt ) then
2018-11-04 02:06:49 +05:30
noOrientationGiven : if ( all ( prm % angles > = 39 9.9_pReal ) ) then
2018-11-04 03:43:20 +05:30
of = mappingHomogenization ( 1 , 1 , e )
2018-11-04 03:13:52 +05:30
dependentState ( instance ) % orientation ( 1 : 3 , 1 : 3 , of ) = math_EulerToR ( math_sampleRandomOri ( ) )
2018-11-04 02:06:49 +05:30
do i = 2_pInt , mesh_NipsPerElem
2018-11-04 03:43:20 +05:30
of = mappingHomogenization ( 1 , i , e )
2018-11-04 03:13:52 +05:30
dependentState ( instance ) % orientation ( 1 : 3 , 1 : 3 , of ) = merge ( dependentState ( instance ) % orientation ( 1 : 3 , 1 : 3 , of ) , &
2018-11-04 02:06:49 +05:30
math_EulerToR ( math_sampleRandomOri ( ) ) , &
microstructure_elemhomo ( mesh_microstructureAt ( e ) ) )
enddo
else noOrientationGiven
do i = 1_pInt , mesh_NipsPerElem
2018-11-04 03:43:20 +05:30
of = mappingHomogenization ( 1 , i , e )
2018-11-04 03:13:52 +05:30
dependentState ( instance ) % orientation ( 1 : 3 , 1 : 3 , of ) = math_EulerToR ( prm % angles * inRad )
2018-11-04 02:06:49 +05:30
enddo
endif noOrientationGiven
endif
enddo elementLooping
2018-08-28 16:45:38 +05:30
end associate
enddo
2014-09-10 19:44:03 +05:30
2013-02-21 03:36:15 +05:30
end subroutine homogenization_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
!--------------------------------------------------------------------------------------------------
2018-11-04 03:53:52 +05:30
subroutine homogenization_RGC_partitionDeformation ( F , avgF , instance , of )
2013-02-21 03:36:15 +05:30
use debug , only : &
debug_level , &
debug_homogenization , &
debug_levelExtensive
use material , only : &
2018-11-04 03:53:52 +05:30
homogenization_maxNgrains
2009-10-12 21:31:49 +05:30
implicit none
2013-02-21 03:36:15 +05:30
real ( pReal ) , dimension ( 3 , 3 , homogenization_maxNgrains ) , intent ( out ) :: F !< partioned F per grain
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: avgF !< averaged F
integer ( pInt ) , intent ( in ) :: &
2018-11-04 03:53:52 +05:30
instance , &
of !< element number
2013-02-21 03:36:15 +05:30
real ( pReal ) , dimension ( 3 ) :: aVect , nVect
2009-10-12 21:31:49 +05:30
integer ( pInt ) , dimension ( 4 ) :: intFace
integer ( pInt ) , dimension ( 3 ) :: iGrain3
2018-11-04 03:53:52 +05:30
integer ( pInt ) :: iGrain , iFace , i , j
2018-11-04 03:43:20 +05:30
type ( tParameters ) :: prm
2009-10-12 21:31:49 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! compute the deformation gradient of individual grains due to relaxations
2018-11-04 03:43:20 +05:30
associate ( prm = > param ( instance ) )
2009-10-12 21:31:49 +05:30
F = 0.0_pReal
2018-11-04 03:53:52 +05:30
do iGrain = 1_pInt , product ( prm % Nconstituents )
2018-08-30 14:20:52 +05:30
iGrain3 = grain1to3 ( iGrain , instance )
2018-11-04 02:06:49 +05:30
do iFace = 1_pInt , 6_pInt
2018-11-04 03:43:20 +05:30
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 )
2012-02-13 19:48:07 +05:30
forall ( i = 1_pInt : 3_pInt , j = 1_pInt : 3_pInt ) &
2018-10-13 21:51:13 +05:30
F ( i , j , iGrain ) = F ( i , j , iGrain ) + aVect ( i ) * nVect ( j ) ! calculating deformation relaxations due to interface relaxation
2009-10-12 21:31:49 +05:30
enddo
2013-02-21 03:36:15 +05:30
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
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! debugging the grain deformation gradients
2012-07-05 15:24:50 +05:30
if ( iand ( debug_level ( debug_homogenization ) , debug_levelExtensive ) / = 0_pInt ) then
openmp parallelization working again (at least for j2 and nonlocal constitutive model).
In order to keep it like that, please follow these simple rules:
DON'T use implicit array subscripts:
example: real, dimension(3,3) :: A,B
A(:,2) = B(:,1) <--- DON'T USE
A(1:3,2) = B(1:3,1) <--- BETTER USE
In many cases the use of explicit array subscripts is inevitable for parallelization. Additionally, it is an easy means to prevent memory leaks.
Enclose all write statements with the following:
!$OMP CRITICAL (write2out)
<your write statement>
!$OMP END CRITICAL (write2out)
Whenever you change something in the code and are not sure if it affects parallelization and leads to nonconforming behavior, please ask me and/or Franz to check this.
2011-03-17 16:16:17 +05:30
!$OMP CRITICAL (write2out)
2012-02-01 00:48:55 +05:30
write ( 6 , '(1x,a32,1x,i3)' ) 'Deformation gradient of grain: ' , iGrain
2012-02-13 19:48:07 +05:30
do i = 1_pInt , 3_pInt
2012-02-16 00:28:38 +05:30
write ( 6 , '(1x,3(e15.8,1x))' ) ( F ( i , j , iGrain ) , j = 1_pInt , 3_pInt )
2009-10-12 21:31:49 +05:30
enddo
write ( 6 , * ) ' '
2012-11-06 21:20:20 +05:30
flush ( 6 )
openmp parallelization working again (at least for j2 and nonlocal constitutive model).
In order to keep it like that, please follow these simple rules:
DON'T use implicit array subscripts:
example: real, dimension(3,3) :: A,B
A(:,2) = B(:,1) <--- DON'T USE
A(1:3,2) = B(1:3,1) <--- BETTER USE
In many cases the use of explicit array subscripts is inevitable for parallelization. Additionally, it is an easy means to prevent memory leaks.
Enclose all write statements with the following:
!$OMP CRITICAL (write2out)
<your write statement>
!$OMP END CRITICAL (write2out)
Whenever you change something in the code and are not sure if it affects parallelization and leads to nonconforming behavior, please ask me and/or Franz to check this.
2011-03-17 16:16:17 +05:30
!$OMP END CRITICAL (write2out)
2009-10-12 21:31:49 +05:30
endif
2010-11-26 17:20:20 +05:30
2009-10-12 21:31:49 +05:30
enddo
2018-11-04 03:43:20 +05:30
end associate
2009-10-12 21:31:49 +05:30
2013-02-21 03:36:15 +05:30
end subroutine homogenization_RGC_partitionDeformation
!--------------------------------------------------------------------------------------------------
!> @brief update the internal state of the homogenization scheme and tell whether "done" and
! "happy" with result
!--------------------------------------------------------------------------------------------------
2014-09-19 23:29:06 +05:30
function homogenization_RGC_updateState ( P , F , F0 , avgF , dt , dPdF , ip , el )
2016-05-29 14:15:03 +05:30
use prec , only : &
2016-10-29 13:09:08 +05:30
dEq0
2013-02-21 03:36:15 +05:30
use debug , only : &
debug_level , &
debug_homogenization , &
debug_levelExtensive , &
debug_e , &
debug_i
use math , only : &
math_invert
use mesh , only : &
2018-10-13 21:51:13 +05:30
mesh_homogenizationAt
2013-02-21 03:36:15 +05:30
use material , only : &
homogenization_maxNgrains , &
homogenization_typeInstance , &
2014-08-21 23:18:20 +05:30
homogState , &
2014-09-19 23:29:06 +05:30
mappingHomogenization , &
2013-02-21 03:36:15 +05:30
homogenization_Ngrains
use numerics , only : &
absTol_RGC , &
relTol_RGC , &
absMax_RGC , &
relMax_RGC , &
pPert_RGC , &
maxdRelax_RGC , &
viscPower_RGC , &
viscModus_RGC , &
refRelaxRate_RGC
2009-10-12 21:31:49 +05:30
implicit none
2014-09-19 23:29:06 +05:30
2013-02-21 03:36:15 +05:30
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
2018-11-03 21:20:43 +05:30
2009-10-12 21:31:49 +05:30
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
2018-11-04 00:30:02 +05:30
integer ( pInt ) instance , iNum , i , j , nIntFaceTot , iGrN , iGrP , iMun , iFace , k , l , ipert , iGrain , nGrain , of
2009-12-16 21:50:53 +05:30
real ( pReal ) , dimension ( 3 , 3 , homogenization_maxNgrains ) :: R , pF , pR , D , pD
2010-03-24 18:50:12 +05:30
real ( pReal ) , dimension ( 3 , homogenization_maxNgrains ) :: NN , pNN
2013-12-18 12:58:01 +05:30
real ( pReal ) , dimension ( 3 ) :: normP , normN , mornP , mornN
2018-11-04 00:30:02 +05:30
real ( pReal ) :: residMax , stresMax , volDiscrep
2011-06-06 20:57:35 +05:30
logical error
2013-02-21 03:36:15 +05:30
2009-11-17 19:12:38 +05:30
real ( pReal ) , dimension ( : , : ) , allocatable :: tract , jmatrix , jnverse , smatrix , pmatrix , rmatrix
2009-10-12 21:31:49 +05:30
real ( pReal ) , dimension ( : ) , allocatable :: resid , relax , p_relax , p_resid , drelax
2013-02-13 15:06:06 +05:30
2016-10-29 13:09:08 +05:30
zeroTimeStep : if ( dEq0 ( dt ) ) then
2016-05-29 14:15:03 +05:30
homogenization_RGC_updateState = . true . ! pretend everything is fine and return
2013-02-21 03:36:15 +05:30
return
2016-05-29 14:15:03 +05:30
endif zeroTimeStep
2009-10-12 21:31:49 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! get the dimension of the cluster (grains and interfaces)
2018-10-13 21:51:13 +05:30
instance = homogenization_typeInstance ( mesh_homogenizationAt ( el ) )
2018-11-04 00:30:02 +05:30
of = mappingHomogenization ( 1 , ip , el )
2018-08-30 14:16:30 +05:30
nGDim = param ( instance ) % Nconstituents
2018-10-13 21:51:13 +05:30
nGrain = homogenization_Ngrains ( mesh_homogenizationAt ( el ) )
2012-02-21 22:01:37 +05:30
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 )
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
2013-12-18 12:58:01 +05:30
allocate ( resid ( 3_pInt * nIntFaceTot ) , source = 0.0_pReal )
allocate ( tract ( nIntFaceTot , 3 ) , source = 0.0_pReal )
2018-11-04 00:30:02 +05:30
relax = homogState ( mappingHomogenization ( 2 , ip , el ) ) % state ( 1 : 3_pInt * nIntFaceTot , of )
2018-11-03 21:20:43 +05:30
drelax = relax &
2018-11-04 00:30:02 +05:30
- homogState ( mappingHomogenization ( 2 , ip , el ) ) % state0 ( 1 : 3_pInt * nIntFaceTot , of )
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! debugging the obtained state
2012-07-05 15:24:50 +05:30
if ( iand ( debug_level ( debug_homogenization ) , debug_levelExtensive ) / = 0_pInt ) then
openmp parallelization working again (at least for j2 and nonlocal constitutive model).
In order to keep it like that, please follow these simple rules:
DON'T use implicit array subscripts:
example: real, dimension(3,3) :: A,B
A(:,2) = B(:,1) <--- DON'T USE
A(1:3,2) = B(1:3,1) <--- BETTER USE
In many cases the use of explicit array subscripts is inevitable for parallelization. Additionally, it is an easy means to prevent memory leaks.
Enclose all write statements with the following:
!$OMP CRITICAL (write2out)
<your write statement>
!$OMP END CRITICAL (write2out)
Whenever you change something in the code and are not sure if it affects parallelization and leads to nonconforming behavior, please ask me and/or Franz to check this.
2011-03-17 16:16:17 +05:30
!$OMP CRITICAL (write2out)
2012-02-01 00:48:55 +05:30
write ( 6 , '(1x,a30)' ) 'Obtained state: '
2012-02-21 22:01:37 +05:30
do i = 1_pInt , 3_pInt * nIntFaceTot
2018-11-04 00:30:02 +05:30
write ( 6 , '(1x,2(e15.8,1x))' ) homogState ( mappingHomogenization ( 2 , ip , el ) ) % state ( i , of )
2009-10-12 21:31:49 +05:30
enddo
write ( 6 , * ) ' '
openmp parallelization working again (at least for j2 and nonlocal constitutive model).
In order to keep it like that, please follow these simple rules:
DON'T use implicit array subscripts:
example: real, dimension(3,3) :: A,B
A(:,2) = B(:,1) <--- DON'T USE
A(1:3,2) = B(1:3,1) <--- BETTER USE
In many cases the use of explicit array subscripts is inevitable for parallelization. Additionally, it is an easy means to prevent memory leaks.
Enclose all write statements with the following:
!$OMP CRITICAL (write2out)
<your write statement>
!$OMP END CRITICAL (write2out)
Whenever you change something in the code and are not sure if it affects parallelization and leads to nonconforming behavior, please ask me and/or Franz to check this.
2011-03-17 16:16:17 +05:30
!$OMP END CRITICAL (write2out)
2009-10-12 21:31:49 +05:30
endif
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! computing interface mismatch and stress penalty tensor for all interfaces of all grains
2018-08-30 14:20:52 +05:30
call stressPenalty ( R , NN , avgF , F , ip , el , instance )
2009-12-16 21:50:53 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! calculating volume discrepancy and stress penalty related to overall volume discrepancy
2018-08-30 14:20:52 +05:30
call volumePenalty ( D , volDiscrep , F , avgF , ip , el )
2009-12-16 21:50:53 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! debugging the mismatch, stress and penalties of grains
2012-07-05 15:24:50 +05:30
if ( iand ( debug_level ( debug_homogenization ) , debug_levelExtensive ) / = 0_pInt ) then
openmp parallelization working again (at least for j2 and nonlocal constitutive model).
In order to keep it like that, please follow these simple rules:
DON'T use implicit array subscripts:
example: real, dimension(3,3) :: A,B
A(:,2) = B(:,1) <--- DON'T USE
A(1:3,2) = B(1:3,1) <--- BETTER USE
In many cases the use of explicit array subscripts is inevitable for parallelization. Additionally, it is an easy means to prevent memory leaks.
Enclose all write statements with the following:
!$OMP CRITICAL (write2out)
<your write statement>
!$OMP END CRITICAL (write2out)
Whenever you change something in the code and are not sure if it affects parallelization and leads to nonconforming behavior, please ask me and/or Franz to check this.
2011-03-17 16:16:17 +05:30
!$OMP CRITICAL (write2out)
2012-02-21 22:01:37 +05:30
do iGrain = 1_pInt , nGrain
2013-02-21 03:36:15 +05:30
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
2012-02-21 22:01:37 +05:30
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 ) , &
2018-10-13 21:51:13 +05:30
( R ( i , j , iGrain ) , j = 1_pInt , 3_pInt ) , &
( D ( i , j , iGrain ) , j = 1_pInt , 3_pInt )
2009-10-12 21:31:49 +05:30
enddo
write ( 6 , * ) ' '
enddo
openmp parallelization working again (at least for j2 and nonlocal constitutive model).
In order to keep it like that, please follow these simple rules:
DON'T use implicit array subscripts:
example: real, dimension(3,3) :: A,B
A(:,2) = B(:,1) <--- DON'T USE
A(1:3,2) = B(1:3,1) <--- BETTER USE
In many cases the use of explicit array subscripts is inevitable for parallelization. Additionally, it is an easy means to prevent memory leaks.
Enclose all write statements with the following:
!$OMP CRITICAL (write2out)
<your write statement>
!$OMP END CRITICAL (write2out)
Whenever you change something in the code and are not sure if it affects parallelization and leads to nonconforming behavior, please ask me and/or Franz to check this.
2011-03-17 16:16:17 +05:30
!$OMP END CRITICAL (write2out)
2009-10-12 21:31:49 +05:30
endif
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
2012-02-21 22:01:37 +05:30
do iNum = 1_pInt , nIntFaceTot
2018-08-30 14:20:52 +05:30
faceID = interface1to4 ( iNum , instance ) ! 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)
iGr3N = faceID ( 2 : 4 ) ! identifying the grain ID in local coordinate system (3-dimensional index)
2018-08-30 14:20:52 +05:30
iGrN = grain3to1 ( iGr3N , instance ) ! translate the local grain ID into global coordinate system (1-dimensional index)
intFaceN = getInterface ( 2_pInt * faceID ( 1 ) , iGr3N )
2018-11-04 03:43:20 +05:30
normN = interfaceNormal ( intFaceN , instance , of )
2010-11-26 17:20:20 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! identify the right/up/front grain (+|P)
2009-10-12 21:31:49 +05:30
iGr3P = iGr3N
2013-02-21 03:36:15 +05:30
iGr3P ( faceID ( 1 ) ) = iGr3N ( faceID ( 1 ) ) + 1_pInt ! identifying the grain ID in local coordinate system (3-dimensional index)
2018-08-30 14:20:52 +05:30
iGrP = grain3to1 ( iGr3P , instance ) ! translate the local grain ID into global coordinate system (1-dimensional index)
intFaceP = getInterface ( 2_pInt * faceID ( 1 ) - 1_pInt , iGr3P )
2018-11-04 03:43:20 +05:30
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)
2012-02-21 22:01:37 +05:30
do i = 1_pInt , 3_pInt
2012-02-13 19:48:07 +05:30
tract ( iNum , i ) = sign ( viscModus_RGC * ( abs ( drelax ( i + 3 * ( iNum - 1_pInt ) ) ) / ( refRelaxRate_RGC * dt ) ) ** viscPower_RGC , &
2013-02-21 03:36:15 +05:30
drelax ( i + 3 * ( iNum - 1_pInt ) ) ) ! contribution from the relaxation viscosity
2012-02-21 22:01:37 +05:30
do j = 1_pInt , 3_pInt
2013-02-21 03:36:15 +05:30
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
2009-12-16 21:50:53 +05:30
+ ( P ( i , j , iGrN ) + R ( i , j , iGrN ) + D ( i , j , iGrN ) ) * normN ( j )
2013-02-21 03:36:15 +05:30
resid ( i + 3_pInt * ( iNum - 1_pInt ) ) = tract ( iNum , i ) ! translate the local residual into global 1-dimensional residual array
2009-11-17 19:12:38 +05:30
enddo
2009-10-12 21:31:49 +05:30
enddo
2010-11-26 17:20:20 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! debugging the residual stress
2012-07-05 15:24:50 +05:30
if ( iand ( debug_level ( debug_homogenization ) , debug_levelExtensive ) / = 0_pInt ) then
openmp parallelization working again (at least for j2 and nonlocal constitutive model).
In order to keep it like that, please follow these simple rules:
DON'T use implicit array subscripts:
example: real, dimension(3,3) :: A,B
A(:,2) = B(:,1) <--- DON'T USE
A(1:3,2) = B(1:3,1) <--- BETTER USE
In many cases the use of explicit array subscripts is inevitable for parallelization. Additionally, it is an easy means to prevent memory leaks.
Enclose all write statements with the following:
!$OMP CRITICAL (write2out)
<your write statement>
!$OMP END CRITICAL (write2out)
Whenever you change something in the code and are not sure if it affects parallelization and leads to nonconforming behavior, please ask me and/or Franz to check this.
2011-03-17 16:16:17 +05:30
!$OMP CRITICAL (write2out)
2012-02-01 00:48:55 +05:30
write ( 6 , '(1x,a30,1x,i3)' ) 'Traction at interface: ' , iNum
2012-02-21 22:01:37 +05:30
write ( 6 , '(1x,3(e15.8,1x))' ) ( tract ( iNum , j ) , j = 1_pInt , 3_pInt )
2009-10-12 21:31:49 +05:30
write ( 6 , * ) ' '
openmp parallelization working again (at least for j2 and nonlocal constitutive model).
In order to keep it like that, please follow these simple rules:
DON'T use implicit array subscripts:
example: real, dimension(3,3) :: A,B
A(:,2) = B(:,1) <--- DON'T USE
A(1:3,2) = B(1:3,1) <--- BETTER USE
In many cases the use of explicit array subscripts is inevitable for parallelization. Additionally, it is an easy means to prevent memory leaks.
Enclose all write statements with the following:
!$OMP CRITICAL (write2out)
<your write statement>
!$OMP END CRITICAL (write2out)
Whenever you change something in the code and are not sure if it affects parallelization and leads to nonconforming behavior, please ask me and/or Franz to check this.
2011-03-17 16:16:17 +05:30
!$OMP END CRITICAL (write2out)
2009-10-12 21:31:49 +05:30
endif
enddo
2013-12-18 12:58:01 +05:30
!--------------------------------------------------------------------------------------------------
2013-02-21 03:36:15 +05:30
! 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
2010-11-26 17:20:20 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! Debugging the convergent criteria
2012-07-05 15:24:50 +05:30
if ( iand ( debug_level ( debug_homogenization ) , debug_levelExtensive ) / = 0_pInt &
2012-03-09 01:55:28 +05:30
. and . debug_e == el . and . debug_i == ip ) then
openmp parallelization working again (at least for j2 and nonlocal constitutive model).
In order to keep it like that, please follow these simple rules:
DON'T use implicit array subscripts:
example: real, dimension(3,3) :: A,B
A(:,2) = B(:,1) <--- DON'T USE
A(1:3,2) = B(1:3,1) <--- BETTER USE
In many cases the use of explicit array subscripts is inevitable for parallelization. Additionally, it is an easy means to prevent memory leaks.
Enclose all write statements with the following:
!$OMP CRITICAL (write2out)
<your write statement>
!$OMP END CRITICAL (write2out)
Whenever you change something in the code and are not sure if it affects parallelization and leads to nonconforming behavior, please ask me and/or Franz to check this.
2011-03-17 16:16:17 +05:30
!$OMP CRITICAL (write2out)
2012-02-01 00:48:55 +05:30
write ( 6 , '(1x,a)' ) ' '
write ( 6 , '(1x,a,1x,i2,1x,i4)' ) 'RGC residual check ...' , ip , el
2012-02-16 00:28:38 +05:30
write ( 6 , '(1x,a15,1x,e15.8,1x,a7,i3,1x,a12,i2,i2)' ) 'Max stress: ' , stresMax , &
2009-10-12 21:31:49 +05:30
'@ grain' , stresLoc ( 3 ) , 'in component' , stresLoc ( 1 ) , stresLoc ( 2 )
2012-02-16 00:28:38 +05:30
write ( 6 , '(1x,a15,1x,e15.8,1x,a7,i3,1x,a12,i2)' ) 'Max residual: ' , residMax , &
2009-10-12 21:31:49 +05:30
'@ iface' , residLoc ( 1 ) , 'in direction' , residLoc ( 2 )
2012-11-06 21:20:20 +05:30
flush ( 6 )
openmp parallelization working again (at least for j2 and nonlocal constitutive model).
In order to keep it like that, please follow these simple rules:
DON'T use implicit array subscripts:
example: real, dimension(3,3) :: A,B
A(:,2) = B(:,1) <--- DON'T USE
A(1:3,2) = B(1:3,1) <--- BETTER USE
In many cases the use of explicit array subscripts is inevitable for parallelization. Additionally, it is an easy means to prevent memory leaks.
Enclose all write statements with the following:
!$OMP CRITICAL (write2out)
<your write statement>
!$OMP END CRITICAL (write2out)
Whenever you change something in the code and are not sure if it affects parallelization and leads to nonconforming behavior, please ask me and/or Franz to check this.
2011-03-17 16:16:17 +05:30
!$OMP END CRITICAL (write2out)
2009-10-12 21:31:49 +05:30
endif
2010-11-26 17:20:20 +05:30
2009-10-12 21:31:49 +05:30
homogenization_RGC_updateState = . false .
2013-10-16 18:34:59 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! If convergence reached => done and happy
2009-10-12 21:31:49 +05:30
if ( residMax < relTol_RGC * stresMax . or . residMax < absTol_RGC ) then
homogenization_RGC_updateState = . true .
2010-11-26 17:20:20 +05:30
2012-07-05 15:24:50 +05:30
if ( iand ( debug_level ( debug_homogenization ) , debug_levelExtensive ) / = 0_pInt &
2012-03-09 01:55:28 +05:30
. and . debug_e == el . and . debug_i == ip ) then
openmp parallelization working again (at least for j2 and nonlocal constitutive model).
In order to keep it like that, please follow these simple rules:
DON'T use implicit array subscripts:
example: real, dimension(3,3) :: A,B
A(:,2) = B(:,1) <--- DON'T USE
A(1:3,2) = B(1:3,1) <--- BETTER USE
In many cases the use of explicit array subscripts is inevitable for parallelization. Additionally, it is an easy means to prevent memory leaks.
Enclose all write statements with the following:
!$OMP CRITICAL (write2out)
<your write statement>
!$OMP END CRITICAL (write2out)
Whenever you change something in the code and are not sure if it affects parallelization and leads to nonconforming behavior, please ask me and/or Franz to check this.
2011-03-17 16:16:17 +05:30
!$OMP CRITICAL (write2out)
2013-10-16 18:34:59 +05:30
write ( 6 , '(1x,a55,/)' ) '... done and happy'
2012-11-06 21:20:20 +05:30
flush ( 6 )
openmp parallelization working again (at least for j2 and nonlocal constitutive model).
In order to keep it like that, please follow these simple rules:
DON'T use implicit array subscripts:
example: real, dimension(3,3) :: A,B
A(:,2) = B(:,1) <--- DON'T USE
A(1:3,2) = B(1:3,1) <--- BETTER USE
In many cases the use of explicit array subscripts is inevitable for parallelization. Additionally, it is an easy means to prevent memory leaks.
Enclose all write statements with the following:
!$OMP CRITICAL (write2out)
<your write statement>
!$OMP END CRITICAL (write2out)
Whenever you change something in the code and are not sure if it affects parallelization and leads to nonconforming behavior, please ask me and/or Franz to check this.
2011-03-17 16:16:17 +05:30
!$OMP END CRITICAL (write2out)
2009-10-12 21:31:49 +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
2018-11-04 03:53:52 +05:30
do iGrain = 1_pInt , homogenization_Ngrains ( mesh_homogenizationAt ( el ) ) ! time-integration loop for work and energy
do i = 1_pInt , 3_pInt ; do j = 1_pInt , 3_pInt
2018-11-04 00:30:02 +05:30
state ( instance ) % work ( of ) = state ( instance ) % work ( of ) &
+ P ( i , j , iGrain ) * ( F ( i , j , iGrain ) - F0 ( i , j , iGrain ) ) / real ( nGrain , pReal )
state ( instance ) % penaltyEnergy ( of ) = state ( instance ) % penaltyEnergy ( of ) &
+ R ( i , j , iGrain ) * ( F ( i , j , iGrain ) - F0 ( i , j , iGrain ) ) / real ( nGrain , pReal )
2018-11-04 03:53:52 +05:30
enddo ; enddo
2009-10-12 21:31:49 +05:30
enddo
2018-11-04 02:06:49 +05:30
2018-11-04 03:53:52 +05:30
state ( instance ) % mismatch ( 1 : 3 , of ) = sum ( NN , 2 ) / real ( nGrain , pReal ) ! the overall mismatch of all interface normals
2018-11-04 02:06:49 +05:30
state ( instance ) % volumeDiscrepancy ( of ) = volDiscrep
state ( instance ) % relaxationRate_avg ( of ) = sum ( abs ( drelax ) ) / dt / real ( 3_pInt * nIntFaceTot , pReal )
state ( instance ) % relaxationRate_max ( of ) = maxval ( abs ( drelax ) ) / dt
2014-08-21 23:18:20 +05:30
2012-07-05 15:24:50 +05:30
if ( iand ( debug_level ( debug_homogenization ) , debug_levelExtensive ) / = 0_pInt &
2012-03-09 01:55:28 +05:30
. and . debug_e == el . and . debug_i == ip ) then
openmp parallelization working again (at least for j2 and nonlocal constitutive model).
In order to keep it like that, please follow these simple rules:
DON'T use implicit array subscripts:
example: real, dimension(3,3) :: A,B
A(:,2) = B(:,1) <--- DON'T USE
A(1:3,2) = B(1:3,1) <--- BETTER USE
In many cases the use of explicit array subscripts is inevitable for parallelization. Additionally, it is an easy means to prevent memory leaks.
Enclose all write statements with the following:
!$OMP CRITICAL (write2out)
<your write statement>
!$OMP END CRITICAL (write2out)
Whenever you change something in the code and are not sure if it affects parallelization and leads to nonconforming behavior, please ask me and/or Franz to check this.
2011-03-17 16:16:17 +05:30
!$OMP CRITICAL (write2out)
2018-11-04 00:30:02 +05:30
write ( 6 , '(1x,a30,1x,e15.8)' ) 'Constitutive work: ' , state ( instance ) % work ( of )
2018-11-04 02:06:49 +05:30
write ( 6 , '(1x,a30,3(1x,e15.8))' ) 'Magnitude mismatch: ' , state ( instance ) % mismatch ( 1 , of ) , &
state ( instance ) % mismatch ( 2 , of ) , &
state ( instance ) % mismatch ( 3 , of )
write ( 6 , '(1x,a30,1x,e15.8)' ) 'Penalty energy: ' , state ( instance ) % penaltyEnergy ( of )
write ( 6 , '(1x,a30,1x,e15.8,/)' ) 'Volume discrepancy: ' , state ( instance ) % volumeDiscrepancy ( of )
write ( 6 , '(1x,a30,1x,e15.8)' ) 'Maximum relaxation rate: ' , state ( instance ) % relaxationRate_max ( of )
write ( 6 , '(1x,a30,1x,e15.8,/)' ) 'Average relaxation rate: ' , state ( instance ) % relaxationRate_avg ( of )
2012-11-06 21:20:20 +05:30
flush ( 6 )
openmp parallelization working again (at least for j2 and nonlocal constitutive model).
In order to keep it like that, please follow these simple rules:
DON'T use implicit array subscripts:
example: real, dimension(3,3) :: A,B
A(:,2) = B(:,1) <--- DON'T USE
A(1:3,2) = B(1:3,1) <--- BETTER USE
In many cases the use of explicit array subscripts is inevitable for parallelization. Additionally, it is an easy means to prevent memory leaks.
Enclose all write statements with the following:
!$OMP CRITICAL (write2out)
<your write statement>
!$OMP END CRITICAL (write2out)
Whenever you change something in the code and are not sure if it affects parallelization and leads to nonconforming behavior, please ask me and/or Franz to check this.
2011-03-17 16:16:17 +05:30
!$OMP END CRITICAL (write2out)
2009-10-12 21:31:49 +05:30
endif
2010-11-26 17:20:20 +05:30
2009-10-12 21:31:49 +05:30
return
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! 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
2010-11-26 17:20:20 +05:30
2012-07-05 15:24:50 +05:30
if ( iand ( debug_level ( debug_homogenization ) , debug_levelExtensive ) / = 0_pInt &
2012-03-09 01:55:28 +05:30
. and . debug_e == el . and . debug_i == ip ) then
openmp parallelization working again (at least for j2 and nonlocal constitutive model).
In order to keep it like that, please follow these simple rules:
DON'T use implicit array subscripts:
example: real, dimension(3,3) :: A,B
A(:,2) = B(:,1) <--- DON'T USE
A(1:3,2) = B(1:3,1) <--- BETTER USE
In many cases the use of explicit array subscripts is inevitable for parallelization. Additionally, it is an easy means to prevent memory leaks.
Enclose all write statements with the following:
!$OMP CRITICAL (write2out)
<your write statement>
!$OMP END CRITICAL (write2out)
Whenever you change something in the code and are not sure if it affects parallelization and leads to nonconforming behavior, please ask me and/or Franz to check this.
2011-03-17 16:16:17 +05:30
!$OMP CRITICAL (write2out)
2013-10-16 18:34:59 +05:30
write ( 6 , '(1x,a55,/)' ) '... broken'
2012-11-06 21:20:20 +05:30
flush ( 6 )
openmp parallelization working again (at least for j2 and nonlocal constitutive model).
In order to keep it like that, please follow these simple rules:
DON'T use implicit array subscripts:
example: real, dimension(3,3) :: A,B
A(:,2) = B(:,1) <--- DON'T USE
A(1:3,2) = B(1:3,1) <--- BETTER USE
In many cases the use of explicit array subscripts is inevitable for parallelization. Additionally, it is an easy means to prevent memory leaks.
Enclose all write statements with the following:
!$OMP CRITICAL (write2out)
<your write statement>
!$OMP END CRITICAL (write2out)
Whenever you change something in the code and are not sure if it affects parallelization and leads to nonconforming behavior, please ask me and/or Franz to check this.
2011-03-17 16:16:17 +05:30
!$OMP END CRITICAL (write2out)
2009-10-12 21:31:49 +05:30
endif
2010-11-26 17:20:20 +05:30
2009-10-12 21:31:49 +05:30
return
2013-02-21 03:36:15 +05:30
else ! proceed with computing the Jacobian and state update
2012-07-05 15:24:50 +05:30
if ( iand ( debug_level ( debug_homogenization ) , debug_levelExtensive ) / = 0_pInt &
2012-03-09 01:55:28 +05:30
. and . debug_e == el . and . debug_i == ip ) then
openmp parallelization working again (at least for j2 and nonlocal constitutive model).
In order to keep it like that, please follow these simple rules:
DON'T use implicit array subscripts:
example: real, dimension(3,3) :: A,B
A(:,2) = B(:,1) <--- DON'T USE
A(1:3,2) = B(1:3,1) <--- BETTER USE
In many cases the use of explicit array subscripts is inevitable for parallelization. Additionally, it is an easy means to prevent memory leaks.
Enclose all write statements with the following:
!$OMP CRITICAL (write2out)
<your write statement>
!$OMP END CRITICAL (write2out)
Whenever you change something in the code and are not sure if it affects parallelization and leads to nonconforming behavior, please ask me and/or Franz to check this.
2011-03-17 16:16:17 +05:30
!$OMP CRITICAL (write2out)
2013-10-16 18:34:59 +05:30
write ( 6 , '(1x,a55,/)' ) '... not yet done'
2012-11-06 21:20:20 +05:30
flush ( 6 )
openmp parallelization working again (at least for j2 and nonlocal constitutive model).
In order to keep it like that, please follow these simple rules:
DON'T use implicit array subscripts:
example: real, dimension(3,3) :: A,B
A(:,2) = B(:,1) <--- DON'T USE
A(1:3,2) = B(1:3,1) <--- BETTER USE
In many cases the use of explicit array subscripts is inevitable for parallelization. Additionally, it is an easy means to prevent memory leaks.
Enclose all write statements with the following:
!$OMP CRITICAL (write2out)
<your write statement>
!$OMP END CRITICAL (write2out)
Whenever you change something in the code and are not sure if it affects parallelization and leads to nonconforming behavior, please ask me and/or Franz to check this.
2011-03-17 16:16:17 +05:30
!$OMP END CRITICAL (write2out)
2009-10-12 21:31:49 +05:30
endif
2010-11-26 17:20:20 +05:30
2009-10-12 21:31:49 +05:30
endif
2014-09-19 23:29:06 +05:30
!---------------------------------------------------------------------------------------------------
2013-02-21 03:36:15 +05:30
! 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"
2013-12-18 12:58:01 +05:30
allocate ( smatrix ( 3 * nIntFaceTot , 3 * nIntFaceTot ) , source = 0.0_pReal )
2012-02-21 22:01:37 +05:30
do iNum = 1_pInt , nIntFaceTot
2018-08-30 14:20:52 +05:30
faceID = interface1to4 ( iNum , instance ) ! assembling of local dPdF into global Jacobian matrix
2010-11-26 17:20:20 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! identify the left/bottom/back grain (-|N)
iGr3N = faceID ( 2 : 4 ) ! identifying the grain ID in local coordinate sytem
2018-08-30 14:20:52 +05:30
iGrN = grain3to1 ( iGr3N , instance ) ! translate into global grain ID
intFaceN = getInterface ( 2_pInt * faceID ( 1 ) , iGr3N ) ! identifying the connecting interface in local coordinate system
2018-11-04 03:43:20 +05:30
normN = interfaceNormal ( intFaceN , instance , of )
2018-11-04 02:06:49 +05:30
do iFace = 1_pInt , 6_pInt
2018-08-30 14:20:52 +05:30
intFaceN = getInterface ( iFace , iGr3N ) ! identifying all interfaces that influence relaxation of the above interface
2018-11-04 03:43:20 +05:30
mornN = interfaceNormal ( intFaceN , instance , of )
2018-08-30 14:20:52 +05:30
iMun = interface4to1 ( intFaceN , instance ) ! translate the interfaces ID into local 4-dimensional index
2013-12-18 12:58:01 +05:30
if ( iMun > 0 ) then ! get the corresponding tangent
2012-02-21 22:01:37 +05:30
do i = 1_pInt , 3_pInt ; do j = 1_pInt , 3_pInt ; do k = 1_pInt , 3_pInt ; do l = 1_pInt , 3_pInt
2018-11-04 03:53:52 +05:30
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 )
2011-08-03 23:28:16 +05:30
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
2009-10-12 21:31:49 +05:30
endif
enddo
2010-11-26 17:20:20 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! identify the right/up/front grain (+|P)
2009-10-12 21:31:49 +05:30
iGr3P = iGr3N
2013-02-21 03:36:15 +05:30
iGr3P ( faceID ( 1 ) ) = iGr3N ( faceID ( 1 ) ) + 1_pInt ! identifying the grain ID in local coordinate sytem
2018-08-30 14:20:52 +05:30
iGrP = grain3to1 ( iGr3P , instance ) ! translate into global grain ID
intFaceP = getInterface ( 2_pInt * faceID ( 1 ) - 1_pInt , iGr3P ) ! identifying the connecting interface in local coordinate system
2018-11-04 03:43:20 +05:30
normP = interfaceNormal ( intFaceP , instance , of )
2018-11-04 02:06:49 +05:30
do iFace = 1_pInt , 6_pInt
2018-08-30 14:20:52 +05:30
intFaceP = getInterface ( iFace , iGr3P ) ! identifying all interfaces that influence relaxation of the above interface
2018-11-04 03:43:20 +05:30
mornP = interfaceNormal ( intFaceP , instance , of )
2018-08-30 14:20:52 +05:30
iMun = interface4to1 ( intFaceP , instance ) ! translate the interfaces ID into local 4-dimensional index
2013-02-21 03:36:15 +05:30
if ( iMun > 0_pInt ) then ! get the corresponding tangent
2012-02-21 22:01:37 +05:30
do i = 1_pInt , 3_pInt ; do j = 1_pInt , 3_pInt ; do k = 1_pInt , 3_pInt ; do l = 1_pInt , 3_pInt
2011-08-03 23:28:16 +05:30
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
2009-10-12 21:31:49 +05:30
endif
enddo
enddo
2010-11-26 17:20:20 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! debugging the global Jacobian matrix of stress tangent
2012-07-05 15:24:50 +05:30
if ( iand ( debug_level ( debug_homogenization ) , debug_levelExtensive ) / = 0_pInt ) then
openmp parallelization working again (at least for j2 and nonlocal constitutive model).
In order to keep it like that, please follow these simple rules:
DON'T use implicit array subscripts:
example: real, dimension(3,3) :: A,B
A(:,2) = B(:,1) <--- DON'T USE
A(1:3,2) = B(1:3,1) <--- BETTER USE
In many cases the use of explicit array subscripts is inevitable for parallelization. Additionally, it is an easy means to prevent memory leaks.
Enclose all write statements with the following:
!$OMP CRITICAL (write2out)
<your write statement>
!$OMP END CRITICAL (write2out)
Whenever you change something in the code and are not sure if it affects parallelization and leads to nonconforming behavior, please ask me and/or Franz to check this.
2011-03-17 16:16:17 +05:30
!$OMP CRITICAL (write2out)
2012-02-01 00:48:55 +05:30
write ( 6 , '(1x,a30)' ) 'Jacobian matrix of stress'
2012-02-21 22:01:37 +05:30
do i = 1_pInt , 3_pInt * nIntFaceTot
write ( 6 , '(1x,100(e11.4,1x))' ) ( smatrix ( i , j ) , j = 1_pInt , 3_pInt * nIntFaceTot )
2009-10-12 21:31:49 +05:30
enddo
write ( 6 , * ) ' '
2012-11-06 21:20:20 +05:30
flush ( 6 )
openmp parallelization working again (at least for j2 and nonlocal constitutive model).
In order to keep it like that, please follow these simple rules:
DON'T use implicit array subscripts:
example: real, dimension(3,3) :: A,B
A(:,2) = B(:,1) <--- DON'T USE
A(1:3,2) = B(1:3,1) <--- BETTER USE
In many cases the use of explicit array subscripts is inevitable for parallelization. Additionally, it is an easy means to prevent memory leaks.
Enclose all write statements with the following:
!$OMP CRITICAL (write2out)
<your write statement>
!$OMP END CRITICAL (write2out)
Whenever you change something in the code and are not sure if it affects parallelization and leads to nonconforming behavior, please ask me and/or Franz to check this.
2011-03-17 16:16:17 +05:30
!$OMP END CRITICAL (write2out)
2009-10-12 21:31:49 +05:30
endif
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! ... of the stress penalty tangent (mismatch penalty and volume penalty, computed using numerical
! perturbation method) "pmatrix"
2013-12-18 12:58:01 +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 )
2012-02-21 22:01:37 +05:30
do ipert = 1_pInt , 3_pInt * nIntFaceTot
2009-10-12 21:31:49 +05:30
p_relax = relax
2013-02-21 03:36:15 +05:30
p_relax ( ipert ) = relax ( ipert ) + pPert_RGC ! perturb the relaxation vector
2018-11-04 00:30:02 +05:30
homogState ( mappingHomogenization ( 2 , ip , el ) ) % state ( 1 : 3 * nIntFaceTot , of ) = p_relax
2018-11-04 11:57:25 +05:30
call grainDeformation ( pF , avgF , instance , of ) ! compute the grains deformation from perturbed state
2018-08-30 14:20:52 +05:30
call stressPenalty ( pR , pNN , avgF , pF , ip , el , instance ) ! compute stress penalty due to interface mismatch from perturbed state
call volumePenalty ( pD , volDiscrep , pF , avgF , ip , el ) ! compute 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
2009-10-12 21:31:49 +05:30
p_resid = 0.0_pReal
2012-02-21 22:01:37 +05:30
do iNum = 1_pInt , nIntFaceTot
2018-08-30 14:20:52 +05:30
faceID = interface1to4 ( iNum , instance ) ! 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)
iGr3N = faceID ( 2 : 4 ) ! identifying the grain ID in local coordinate system (3-dimensional index)
2018-08-30 14:20:52 +05:30
iGrN = grain3to1 ( iGr3N , instance ) ! translate the local grain ID into global coordinate system (1-dimensional index)
intFaceN = getInterface ( 2_pInt * faceID ( 1 ) , iGr3N ) ! identifying the interface ID of the grain
2018-11-04 03:43:20 +05:30
normN = interfaceNormal ( intFaceN , instance , of )
2010-11-26 17:20:20 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! identify the right/up/front grain (+|P)
2009-10-12 21:31:49 +05:30
iGr3P = iGr3N
2013-02-21 03:36:15 +05:30
iGr3P ( faceID ( 1 ) ) = iGr3N ( faceID ( 1 ) ) + 1_pInt ! identifying the grain ID in local coordinate system (3-dimensional index)
2018-08-30 14:20:52 +05:30
iGrP = grain3to1 ( iGr3P , instance ) ! translate the local grain ID into global coordinate system (1-dimensional index)
intFaceP = getInterface ( 2_pInt * faceID ( 1 ) - 1_pInt , iGr3P ) ! identifying the interface ID of the grain
2018-11-04 03:43:20 +05:30
normP = interfaceNormal ( intFaceP , instance , of )
2010-11-26 17:20:20 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! 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
2009-10-12 21:31:49 +05:30
p_resid ( i + 3 * ( iNum - 1 ) ) = p_resid ( i + 3 * ( iNum - 1 ) ) + ( pR ( i , j , iGrP ) - R ( i , j , iGrP ) ) * normP ( j ) &
2009-12-16 21:50:53 +05:30
+ ( 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 )
2013-02-21 03:36:15 +05:30
enddo ; enddo
2009-10-12 21:31:49 +05:30
enddo
pmatrix ( : , ipert ) = p_resid / pPert_RGC
enddo
2010-11-26 17:20:20 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! debugging the global Jacobian matrix of penalty tangent
2012-07-05 15:24:50 +05:30
if ( iand ( debug_level ( debug_homogenization ) , debug_levelExtensive ) / = 0_pInt ) then
openmp parallelization working again (at least for j2 and nonlocal constitutive model).
In order to keep it like that, please follow these simple rules:
DON'T use implicit array subscripts:
example: real, dimension(3,3) :: A,B
A(:,2) = B(:,1) <--- DON'T USE
A(1:3,2) = B(1:3,1) <--- BETTER USE
In many cases the use of explicit array subscripts is inevitable for parallelization. Additionally, it is an easy means to prevent memory leaks.
Enclose all write statements with the following:
!$OMP CRITICAL (write2out)
<your write statement>
!$OMP END CRITICAL (write2out)
Whenever you change something in the code and are not sure if it affects parallelization and leads to nonconforming behavior, please ask me and/or Franz to check this.
2011-03-17 16:16:17 +05:30
!$OMP CRITICAL (write2out)
2012-02-01 00:48:55 +05:30
write ( 6 , '(1x,a30)' ) 'Jacobian matrix of penalty'
2012-02-21 22:01:37 +05:30
do i = 1_pInt , 3_pInt * nIntFaceTot
write ( 6 , '(1x,100(e11.4,1x))' ) ( pmatrix ( i , j ) , j = 1_pInt , 3_pInt * nIntFaceTot )
2009-10-12 21:31:49 +05:30
enddo
write ( 6 , * ) ' '
2012-11-06 21:20:20 +05:30
flush ( 6 )
openmp parallelization working again (at least for j2 and nonlocal constitutive model).
In order to keep it like that, please follow these simple rules:
DON'T use implicit array subscripts:
example: real, dimension(3,3) :: A,B
A(:,2) = B(:,1) <--- DON'T USE
A(1:3,2) = B(1:3,1) <--- BETTER USE
In many cases the use of explicit array subscripts is inevitable for parallelization. Additionally, it is an easy means to prevent memory leaks.
Enclose all write statements with the following:
!$OMP CRITICAL (write2out)
<your write statement>
!$OMP END CRITICAL (write2out)
Whenever you change something in the code and are not sure if it affects parallelization and leads to nonconforming behavior, please ask me and/or Franz to check this.
2011-03-17 16:16:17 +05:30
!$OMP END CRITICAL (write2out)
2009-10-12 21:31:49 +05:30
endif
2009-11-17 19:12:38 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! ... of the numerical viscosity traction "rmatrix"
2013-12-18 12:58:01 +05:30
allocate ( rmatrix ( 3 * nIntFaceTot , 3 * nIntFaceTot ) , source = 0.0_pReal )
2012-02-21 22:01:37 +05:30
forall ( i = 1_pInt : 3_pInt * nIntFaceTot ) &
2013-02-21 03:36:15 +05:30
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
2010-11-26 17:20:20 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! debugging the global Jacobian matrix of numerical viscosity tangent
2012-07-05 15:24:50 +05:30
if ( iand ( debug_level ( debug_homogenization ) , debug_levelExtensive ) / = 0_pInt ) then
openmp parallelization working again (at least for j2 and nonlocal constitutive model).
In order to keep it like that, please follow these simple rules:
DON'T use implicit array subscripts:
example: real, dimension(3,3) :: A,B
A(:,2) = B(:,1) <--- DON'T USE
A(1:3,2) = B(1:3,1) <--- BETTER USE
In many cases the use of explicit array subscripts is inevitable for parallelization. Additionally, it is an easy means to prevent memory leaks.
Enclose all write statements with the following:
!$OMP CRITICAL (write2out)
<your write statement>
!$OMP END CRITICAL (write2out)
Whenever you change something in the code and are not sure if it affects parallelization and leads to nonconforming behavior, please ask me and/or Franz to check this.
2011-03-17 16:16:17 +05:30
!$OMP CRITICAL (write2out)
2012-02-01 00:48:55 +05:30
write ( 6 , '(1x,a30)' ) 'Jacobian matrix of penalty'
2012-02-21 22:01:37 +05:30
do i = 1_pInt , 3_pInt * nIntFaceTot
write ( 6 , '(1x,100(e11.4,1x))' ) ( rmatrix ( i , j ) , j = 1_pInt , 3_pInt * nIntFaceTot )
2009-11-17 19:12:38 +05:30
enddo
write ( 6 , * ) ' '
2012-11-06 21:20:20 +05:30
flush ( 6 )
openmp parallelization working again (at least for j2 and nonlocal constitutive model).
In order to keep it like that, please follow these simple rules:
DON'T use implicit array subscripts:
example: real, dimension(3,3) :: A,B
A(:,2) = B(:,1) <--- DON'T USE
A(1:3,2) = B(1:3,1) <--- BETTER USE
In many cases the use of explicit array subscripts is inevitable for parallelization. Additionally, it is an easy means to prevent memory leaks.
Enclose all write statements with the following:
!$OMP CRITICAL (write2out)
<your write statement>
!$OMP END CRITICAL (write2out)
Whenever you change something in the code and are not sure if it affects parallelization and leads to nonconforming behavior, please ask me and/or Franz to check this.
2011-03-17 16:16:17 +05:30
!$OMP END CRITICAL (write2out)
2009-11-17 19:12:38 +05:30
endif
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! The overall Jacobian matrix summarizing contributions of smatrix, pmatrix, rmatrix
2009-11-17 19:12:38 +05:30
allocate ( jmatrix ( 3 * nIntFaceTot , 3 * nIntFaceTot ) ) ; jmatrix = smatrix + pmatrix + rmatrix
2010-11-26 17:20:20 +05:30
2012-07-05 15:24:50 +05:30
if ( iand ( debug_level ( debug_homogenization ) , debug_levelExtensive ) / = 0_pInt ) then
openmp parallelization working again (at least for j2 and nonlocal constitutive model).
In order to keep it like that, please follow these simple rules:
DON'T use implicit array subscripts:
example: real, dimension(3,3) :: A,B
A(:,2) = B(:,1) <--- DON'T USE
A(1:3,2) = B(1:3,1) <--- BETTER USE
In many cases the use of explicit array subscripts is inevitable for parallelization. Additionally, it is an easy means to prevent memory leaks.
Enclose all write statements with the following:
!$OMP CRITICAL (write2out)
<your write statement>
!$OMP END CRITICAL (write2out)
Whenever you change something in the code and are not sure if it affects parallelization and leads to nonconforming behavior, please ask me and/or Franz to check this.
2011-03-17 16:16:17 +05:30
!$OMP CRITICAL (write2out)
2012-02-01 00:48:55 +05:30
write ( 6 , '(1x,a30)' ) 'Jacobian matrix (total)'
2012-02-21 22:01:37 +05:30
do i = 1_pInt , 3_pInt * nIntFaceTot
write ( 6 , '(1x,100(e11.4,1x))' ) ( jmatrix ( i , j ) , j = 1_pInt , 3_pInt * nIntFaceTot )
2009-10-12 21:31:49 +05:30
enddo
write ( 6 , * ) ' '
2012-11-06 21:20:20 +05:30
flush ( 6 )
openmp parallelization working again (at least for j2 and nonlocal constitutive model).
In order to keep it like that, please follow these simple rules:
DON'T use implicit array subscripts:
example: real, dimension(3,3) :: A,B
A(:,2) = B(:,1) <--- DON'T USE
A(1:3,2) = B(1:3,1) <--- BETTER USE
In many cases the use of explicit array subscripts is inevitable for parallelization. Additionally, it is an easy means to prevent memory leaks.
Enclose all write statements with the following:
!$OMP CRITICAL (write2out)
<your write statement>
!$OMP END CRITICAL (write2out)
Whenever you change something in the code and are not sure if it affects parallelization and leads to nonconforming behavior, please ask me and/or Franz to check this.
2011-03-17 16:16:17 +05:30
!$OMP END CRITICAL (write2out)
2009-10-12 21:31:49 +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
2013-12-18 12:58:01 +05:30
allocate ( jnverse ( 3_pInt * nIntFaceTot , 3_pInt * nIntFaceTot ) , source = 0.0_pReal )
2013-02-21 03:36:15 +05:30
call math_invert ( size ( jmatrix , 1 ) , jmatrix , jnverse , error ) ! Compute the inverse of the overall Jacobian matrix
2010-11-26 17:20:20 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! debugging the inverse Jacobian matrix
2012-07-05 15:24:50 +05:30
if ( iand ( debug_level ( debug_homogenization ) , debug_levelExtensive ) / = 0_pInt ) then
openmp parallelization working again (at least for j2 and nonlocal constitutive model).
In order to keep it like that, please follow these simple rules:
DON'T use implicit array subscripts:
example: real, dimension(3,3) :: A,B
A(:,2) = B(:,1) <--- DON'T USE
A(1:3,2) = B(1:3,1) <--- BETTER USE
In many cases the use of explicit array subscripts is inevitable for parallelization. Additionally, it is an easy means to prevent memory leaks.
Enclose all write statements with the following:
!$OMP CRITICAL (write2out)
<your write statement>
!$OMP END CRITICAL (write2out)
Whenever you change something in the code and are not sure if it affects parallelization and leads to nonconforming behavior, please ask me and/or Franz to check this.
2011-03-17 16:16:17 +05:30
!$OMP CRITICAL (write2out)
2012-02-01 00:48:55 +05:30
write ( 6 , '(1x,a30)' ) 'Jacobian inverse'
2012-02-21 22:01:37 +05:30
do i = 1_pInt , 3_pInt * nIntFaceTot
2012-02-16 00:28:38 +05:30
write ( 6 , '(1x,100(e11.4,1x))' ) ( jnverse ( i , j ) , j = 1_pInt , 3_pInt * nIntFaceTot )
2009-10-12 21:31:49 +05:30
enddo
write ( 6 , * ) ' '
2012-11-06 21:20:20 +05:30
flush ( 6 )
openmp parallelization working again (at least for j2 and nonlocal constitutive model).
In order to keep it like that, please follow these simple rules:
DON'T use implicit array subscripts:
example: real, dimension(3,3) :: A,B
A(:,2) = B(:,1) <--- DON'T USE
A(1:3,2) = B(1:3,1) <--- BETTER USE
In many cases the use of explicit array subscripts is inevitable for parallelization. Additionally, it is an easy means to prevent memory leaks.
Enclose all write statements with the following:
!$OMP CRITICAL (write2out)
<your write statement>
!$OMP END CRITICAL (write2out)
Whenever you change something in the code and are not sure if it affects parallelization and leads to nonconforming behavior, please ask me and/or Franz to check this.
2011-03-17 16:16:17 +05:30
!$OMP END CRITICAL (write2out)
2009-10-12 21:31:49 +05:30
endif
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! calculate the state update (global relaxation vectors) for the next Newton-Raphson iteration
2009-11-17 19:12:38 +05:30
drelax = 0.0_pReal
2012-02-21 22:01:37 +05:30
do i = 1_pInt , 3_pInt * nIntFaceTot
do j = 1_pInt , 3_pInt * nIntFaceTot
2013-02-21 03:36:15 +05:30
drelax ( i ) = drelax ( i ) - jnverse ( i , j ) * resid ( j ) ! Calculate the correction for the state variable
2009-10-12 21:31:49 +05:30
enddo
enddo
2013-02-21 03:36:15 +05:30
relax = relax + drelax ! Updateing the state variable for the next iteration
2018-11-04 00:30:02 +05:30
homogState ( mappingHomogenization ( 2 , ip , el ) ) % state ( 1 : 3 * nIntFaceTot , of ) = relax
2013-02-21 03:36:15 +05:30
if ( any ( abs ( drelax ) > maxdRelax_RGC ) ) then ! Forcing cutback when the incremental change of relaxation vector becomes too large
homogenization_RGC_updateState = [ . true . , . false . ]
openmp parallelization working again (at least for j2 and nonlocal constitutive model).
In order to keep it like that, please follow these simple rules:
DON'T use implicit array subscripts:
example: real, dimension(3,3) :: A,B
A(:,2) = B(:,1) <--- DON'T USE
A(1:3,2) = B(1:3,1) <--- BETTER USE
In many cases the use of explicit array subscripts is inevitable for parallelization. Additionally, it is an easy means to prevent memory leaks.
Enclose all write statements with the following:
!$OMP CRITICAL (write2out)
<your write statement>
!$OMP END CRITICAL (write2out)
Whenever you change something in the code and are not sure if it affects parallelization and leads to nonconforming behavior, please ask me and/or Franz to check this.
2011-03-17 16:16:17 +05:30
!$OMP CRITICAL (write2out)
2012-02-01 00:48:55 +05:30
write ( 6 , '(1x,a,1x,i3,1x,a,1x,i3,1x,a)' ) 'RGC_updateState: ip' , ip , '| el' , el , 'enforces cutback'
2012-02-16 00:28:38 +05:30
write ( 6 , '(1x,a,1x,e15.8)' ) 'due to large relaxation change =' , maxval ( abs ( drelax ) )
2012-11-06 21:20:20 +05:30
flush ( 6 )
openmp parallelization working again (at least for j2 and nonlocal constitutive model).
In order to keep it like that, please follow these simple rules:
DON'T use implicit array subscripts:
example: real, dimension(3,3) :: A,B
A(:,2) = B(:,1) <--- DON'T USE
A(1:3,2) = B(1:3,1) <--- BETTER USE
In many cases the use of explicit array subscripts is inevitable for parallelization. Additionally, it is an easy means to prevent memory leaks.
Enclose all write statements with the following:
!$OMP CRITICAL (write2out)
<your write statement>
!$OMP END CRITICAL (write2out)
Whenever you change something in the code and are not sure if it affects parallelization and leads to nonconforming behavior, please ask me and/or Franz to check this.
2011-03-17 16:16:17 +05:30
!$OMP END CRITICAL (write2out)
2009-10-12 21:31:49 +05:30
endif
2010-11-26 17:20:20 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! debugging the return state
2012-03-09 01:55:28 +05:30
if ( iand ( debug_homogenization , debug_levelExtensive ) > 0_pInt ) then
openmp parallelization working again (at least for j2 and nonlocal constitutive model).
In order to keep it like that, please follow these simple rules:
DON'T use implicit array subscripts:
example: real, dimension(3,3) :: A,B
A(:,2) = B(:,1) <--- DON'T USE
A(1:3,2) = B(1:3,1) <--- BETTER USE
In many cases the use of explicit array subscripts is inevitable for parallelization. Additionally, it is an easy means to prevent memory leaks.
Enclose all write statements with the following:
!$OMP CRITICAL (write2out)
<your write statement>
!$OMP END CRITICAL (write2out)
Whenever you change something in the code and are not sure if it affects parallelization and leads to nonconforming behavior, please ask me and/or Franz to check this.
2011-03-17 16:16:17 +05:30
!$OMP CRITICAL (write2out)
2012-02-01 00:48:55 +05:30
write ( 6 , '(1x,a30)' ) 'Returned state: '
2012-02-21 22:01:37 +05:30
do i = 1_pInt , 3_pInt * nIntFaceTot
2018-11-04 00:30:02 +05:30
write ( 6 , '(1x,2(e15.8,1x))' ) homogState ( mappingHomogenization ( 2 , ip , el ) ) % state ( i , of )
2009-10-12 21:31:49 +05:30
enddo
write ( 6 , * ) ' '
2012-11-06 21:20:20 +05:30
flush ( 6 )
openmp parallelization working again (at least for j2 and nonlocal constitutive model).
In order to keep it like that, please follow these simple rules:
DON'T use implicit array subscripts:
example: real, dimension(3,3) :: A,B
A(:,2) = B(:,1) <--- DON'T USE
A(1:3,2) = B(1:3,1) <--- BETTER USE
In many cases the use of explicit array subscripts is inevitable for parallelization. Additionally, it is an easy means to prevent memory leaks.
Enclose all write statements with the following:
!$OMP CRITICAL (write2out)
<your write statement>
!$OMP END CRITICAL (write2out)
Whenever you change something in the code and are not sure if it affects parallelization and leads to nonconforming behavior, please ask me and/or Franz to check this.
2011-03-17 16:16:17 +05:30
!$OMP END CRITICAL (write2out)
2009-10-12 21:31:49 +05:30
endif
2018-11-04 01:41:43 +05:30
contains
!--------------------------------------------------------------------------------------------------
!> @brief calculate stress-like penalty due to deformation mismatch
!--------------------------------------------------------------------------------------------------
subroutine stressPenalty ( rPen , nMis , avgF , fDef , ip , el , instance )
use debug , only : &
debug_level , &
debug_homogenization , &
debug_levelExtensive , &
debug_e , &
debug_i
use mesh , only : &
mesh_homogenizationAt
use constitutive , only : &
constitutive_homogenizedC
use math , only : &
math_civita
use material , only : &
2018-11-04 03:43:20 +05:30
homogenization_maxNgrains
2018-11-04 01:41:43 +05:30
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 , instance
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
2018-11-04 03:43:20 +05:30
integer ( pInt ) :: iGrain , iGNghb , iFace , i , j , k , l , of
2018-11-04 01:41:43 +05:30
real ( pReal ) :: muGrain , muGNghb , nDefNorm , bgGrain , bgGNghb
type ( tParameters ) :: prm
real ( pReal ) , parameter :: nDefToler = 1.0e-10_pReal
2018-11-04 03:53:52 +05:30
logical :: debugActive
debugActive = iand ( debug_level ( debug_homogenization ) , debug_levelExtensive ) / = 0_pInt &
. and . debug_e == el . and . debug_i == ip
2018-11-04 01:41:43 +05:30
nGDim = param ( instance ) % Nconstituents
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
2018-11-04 03:43:20 +05:30
of = mappingHomogenization ( 1 , ip , el )
surfCorr = surfaceCorrection ( avgF , instance , of )
2018-11-04 01:41:43 +05:30
associate ( prm = > param ( instance ) )
2018-11-04 03:53:52 +05:30
if ( debugActive ) then
write ( 6 , '(1x,a20,2(1x,i3))' ) 'Correction factor: ' , ip , el
write ( 6 , * ) surfCorr
2018-11-04 01:41:43 +05:30
endif
!--------------------------------------------------------------------------------------------------
! computing the mismatch and penalty stress tensor of all grains
2018-11-04 03:43:20 +05:30
do iGrain = 1_pInt , product ( param ( instance ) % Nconstituents )
2018-11-04 01:41:43 +05:30
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 , instance ) ! get the grain ID in local 3-dimensional index (x,y,z)-position
!* Looping over all six interfaces of each grain
2018-11-04 02:06:49 +05:30
do iFace = 1_pInt , 6_pInt
2018-11-04 01:41:43 +05:30
intFace = getInterface ( iFace , iGrain3 ) ! get the 4-dimensional index of the interface in local numbering system of the grain
2018-11-04 03:43:20 +05:30
nVect = interfaceNormal ( intFace , instance , of )
2018-11-04 01:41:43 +05:30
iGNghb3 = iGrain3 ! identify the neighboring grain across the interface
2018-11-04 02:06:49 +05:30
iGNghb3 ( abs ( intFace ( 1 ) ) ) = iGNghb3 ( abs ( intFace ( 1 ) ) ) &
+ int ( real ( intFace ( 1 ) , pReal ) / real ( abs ( intFace ( 1 ) ) , pReal ) , pInt )
where ( iGNghb3 < 1 ) iGNghb3 = nGDim
where ( iGNghb3 > nGDim ) iGNghb3 = 1_pInt
2018-11-04 01:41:43 +05:30
iGNghb = grain3to1 ( iGNghb3 , instance ) ! get the ID of the neighboring grain
Gmoduli = 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
2018-11-04 03:53:52 +05:30
nDefNorm = nDefNorm + nDef ( i , j ) ** 2.0_pReal ! compute the norm of the mismatch tensor
2018-11-04 01:41:43 +05:30
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)
2018-11-04 03:53:52 +05:30
if ( debugActive ) then
2018-11-04 01:41:43 +05:30
write ( 6 , '(1x,a20,i2,1x,a20,1x,i3)' ) 'Mismatch to face: ' , intFace ( 1 ) , 'neighbor grain: ' , iGNghb
2018-11-04 03:53:52 +05:30
write ( 6 , * ) transpose ( nDef )
write ( 6 , '(1x,a20,e11.4)' ) 'with magnitude: ' , nDefNorm
2018-11-04 01:41:43 +05:30
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 ) * prm % xiAlpha &
* surfCorr ( abs ( intFace ( 1 ) ) ) / prm % dAlpha ( abs ( intFace ( 1 ) ) ) &
* cosh ( prm % ciAlpha * nDefNorm ) &
* 0.5_pReal * nVect ( l ) * nDef ( i , k ) / nDefNorm * math_civita ( k , l , j ) &
* tanh ( nDefNorm / xSmoo_RGC )
enddo ; enddo
enddo ; enddo
enddo
2018-11-04 03:53:52 +05:30
if ( debugActive ) then
write ( 6 , '(1x,a20,i2)' ) 'Penalty of grain: ' , iGrain
write ( 6 , * ) transpose ( rPen ( 1 : 3 , 1 : 3 , iGrain ) )
2018-11-04 01:41:43 +05:30
endif
enddo
end associate
end subroutine stressPenalty
!--------------------------------------------------------------------------------------------------
!> @brief calculate stress-like penalty due to volume discrepancy
!--------------------------------------------------------------------------------------------------
subroutine volumePenalty ( vPen , vDiscrep , fDef , fAvg , ip , el )
use debug , only : &
debug_level , &
debug_homogenization , &
debug_levelExtensive , &
debug_e , &
debug_i
use mesh , only : &
mesh_homogenizationAt
use math , only : &
math_det33 , &
math_inv33
use material , only : &
homogenization_maxNgrains , &
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
2018-11-04 03:53:52 +05:30
integer ( pInt ) :: iGrain , nGrain
logical :: debugActive
debugActive = iand ( debug_level ( debug_homogenization ) , debug_levelExtensive ) / = 0_pInt &
. and . debug_e == el . and . debug_i == ip
2018-11-04 01:41:43 +05:30
nGrain = homogenization_Ngrains ( mesh_homogenizationAt ( el ) )
2018-11-04 03:53:52 +05:30
2018-11-04 01:41:43 +05:30
!--------------------------------------------------------------------------------------------------
! 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 ) ) )
2018-11-04 03:53:52 +05:30
if ( debugActive ) then
2018-11-04 01:41:43 +05:30
write ( 6 , '(1x,a30,i2)' ) 'Volume penalty of grain: ' , iGrain
2018-11-04 03:53:52 +05:30
write ( 6 , * ) transpose ( vPen ( : , : , iGrain ) )
2018-11-04 01:41:43 +05:30
endif
enddo
2018-11-04 02:06:49 +05:30
2018-11-04 01:41:43 +05:30
end subroutine volumePenalty
2018-11-04 02:06:49 +05:30
2018-11-04 01:41:43 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief compute the correction factor accouted for surface evolution (area change) due to
! deformation
!--------------------------------------------------------------------------------------------------
2018-11-04 03:43:20 +05:30
function surfaceCorrection ( avgF , instance , of )
2018-11-04 01:41:43 +05:30
use math , only : &
math_invert33 , &
math_mul33x33
implicit none
real ( pReal ) , dimension ( 3 ) :: surfaceCorrection
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: avgF !< average F
2018-11-04 03:43:20 +05:30
integer ( pInt ) , intent ( in ) :: &
instance , &
of
2018-11-04 01:41:43 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: invC
real ( pReal ) , dimension ( 3 ) :: nVect
real ( pReal ) :: detF
integer ( pInt ) :: i , j , iBase
logical :: error
call math_invert33 ( math_mul33x33 ( transpose ( avgF ) , avgF ) , invC , detF , error )
surfaceCorrection = 0.0_pReal
do iBase = 1_pInt , 3_pInt
2018-11-04 03:43:20 +05:30
nVect = interfaceNormal ( [ iBase , 1_pInt , 1_pInt , 1_pInt ] , instance , of )
2018-11-04 01:41:43 +05:30
do i = 1_pInt , 3_pInt ; do j = 1_pInt , 3_pInt
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
end function surfaceCorrection
2018-11-04 02:06:49 +05:30
2018-11-04 01:41:43 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief compute the equivalent shear and bulk moduli from the elasticity tensor
!--------------------------------------------------------------------------------------------------
function 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 ) :: 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
equivalentModuli ( 1 ) = 0.2_pReal * ( cEquiv_11 - cEquiv_12 ) + 0.6_pReal * cEquiv_44
!--------------------------------------------------------------------------------------------------
! obtain the length of Burgers vector (could be model dependend)
equivalentModuli ( 2 ) = 2.5e-10_pReal
end function equivalentModuli
!--------------------------------------------------------------------------------------------------
!> @brief calculating the grain deformation gradient (the same with
! homogenization_RGC_partionDeformation, but used only for perturbation scheme)
!--------------------------------------------------------------------------------------------------
2018-11-04 11:57:25 +05:30
subroutine grainDeformation ( F , avgF , instance , of )
2018-11-04 01:41:43 +05:30
use material , only : &
2018-11-04 11:57:25 +05:30
homogenization_maxNgrains
2018-11-04 01:41:43 +05:30
implicit none
real ( pReal ) , dimension ( 3 , 3 , homogenization_maxNgrains ) , intent ( out ) :: F !< partioned F per grain
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: avgF !<
integer ( pInt ) , intent ( in ) :: &
2018-11-04 11:57:25 +05:30
instance , &
of
2018-11-04 01:41:43 +05:30
real ( pReal ) , dimension ( 3 ) :: aVect , nVect
integer ( pInt ) , dimension ( 4 ) :: intFace
integer ( pInt ) , dimension ( 3 ) :: iGrain3
2018-11-04 11:57:25 +05:30
integer ( pInt ) :: iGrain , iFace , i , j
2018-11-04 01:41:43 +05:30
2018-11-04 03:43:20 +05:30
!--------------------------------------------------------------------------------------------------
! compute the deformation gradient of individual grains due to relaxations
2018-11-04 11:57:25 +05:30
F = 0.0_pReal
do iGrain = 1_pInt , product ( param ( instance ) % Nconstituents )
iGrain3 = grain1to3 ( iGrain , instance )
do iFace = 1_pInt , 6_pInt
intFace = getInterface ( iFace , iGrain3 )
aVect = relaxationVector ( intFace , instance , of )
nVect = interfaceNormal ( intFace , instance , of )
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
2018-11-04 03:43:20 +05:30
2018-11-04 01:41:43 +05:30
end subroutine grainDeformation
2009-10-12 21:31:49 +05:30
2013-02-21 03:36:15 +05:30
end function homogenization_RGC_updateState
!--------------------------------------------------------------------------------------------------
!> @brief derive average stress and stiffness from constituent quantities
!--------------------------------------------------------------------------------------------------
2018-11-04 00:38:18 +05:30
subroutine homogenization_RGC_averageStressAndItsTangent ( avgP , dAvgPdAvgF , P , dPdF , instance )
2014-08-21 23:18:20 +05:30
use material , only : &
2018-11-04 00:38:18 +05:30
homogenization_maxNgrains
2009-10-12 21:31:49 +05:30
implicit none
2013-10-16 18:34:59 +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
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
2018-11-04 00:38:18 +05:30
integer ( pInt ) , intent ( in ) :: instance
2009-10-12 21:31:49 +05:30
2018-11-04 00:38:18 +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
2013-10-11 21:31:53 +05:30
end subroutine homogenization_RGC_averageStressAndItsTangent
2009-10-12 21:31:49 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief return array of homogenization results for post file inclusion
!--------------------------------------------------------------------------------------------------
2018-10-13 21:51:13 +05:30
pure function homogenization_RGC_postResults ( ip , el ) result ( postResults )
2013-02-21 03:36:15 +05:30
use mesh , only : &
2018-10-13 21:51:13 +05:30
mesh_homogenizationAt
2013-02-21 03:36:15 +05:30
use material , only : &
homogenization_typeInstance , &
2014-08-21 23:18:20 +05:30
homogState , &
2014-09-19 23:29:06 +05:30
mappingHomogenization , &
2013-02-21 03:36:15 +05:30
homogenization_Noutput
2012-03-09 01:55:28 +05:30
2009-10-12 21:31:49 +05:30
implicit none
2013-10-19 00:27:28 +05:30
integer ( pInt ) , intent ( in ) :: &
ip , & !< integration point number
el !< element number
2014-09-19 23:29:06 +05:30
2018-08-30 14:16:30 +05:30
integer ( pInt ) instance , o , c , nIntFaceTot
type ( tParameters ) :: prm
2018-10-13 21:51:13 +05:30
real ( pReal ) , dimension ( homogenization_RGC_sizePostResults ( homogenization_typeInstance ( mesh_homogenizationAt ( el ) ) ) ) :: &
2018-08-30 14:20:52 +05:30
postResults
2009-10-12 21:31:49 +05:30
2018-10-13 21:51:13 +05:30
instance = homogenization_typeInstance ( mesh_homogenizationAt ( el ) )
2018-08-30 14:16:30 +05:30
associate ( prm = > param ( instance ) )
nIntFaceTot = ( prm % Nconstituents ( 1 ) - 1_pInt ) * prm % Nconstituents ( 2 ) * prm % Nconstituents ( 3 ) &
+ prm % Nconstituents ( 1 ) * ( prm % Nconstituents ( 2 ) - 1_pInt ) * prm % Nconstituents ( 3 ) &
+ prm % Nconstituents ( 1 ) * prm % Nconstituents ( 2 ) * ( prm % Nconstituents ( 3 ) - 1_pInt )
2009-10-12 21:31:49 +05:30
c = 0_pInt
2018-08-30 14:20:52 +05:30
postResults = 0.0_pReal
2018-08-31 10:05:30 +05:30
outputsLoop : do o = 1_pInt , size ( prm % outputID )
select case ( prm % outputID ( o ) )
2013-12-18 12:58:01 +05:30
case ( constitutivework_ID )
2018-08-30 14:20:52 +05:30
postResults ( c + 1 ) = homogState ( mappingHomogenization ( 2 , ip , el ) ) % &
2014-08-21 23:18:20 +05:30
state ( 3 * nIntFaceTot + 1 , mappingHomogenization ( 1 , ip , el ) )
2012-02-21 22:01:37 +05:30
c = c + 1_pInt
2013-12-18 12:58:01 +05:30
case ( magnitudemismatch_ID )
2018-08-30 14:20:52 +05:30
postResults ( c + 1 ) = homogState ( mappingHomogenization ( 2 , ip , el ) ) % &
2014-08-21 23:18:20 +05:30
state ( 3 * nIntFaceTot + 2 , mappingHomogenization ( 1 , ip , el ) )
2018-08-30 14:20:52 +05:30
postResults ( c + 2 ) = homogState ( mappingHomogenization ( 2 , ip , el ) ) % &
2014-08-21 23:18:20 +05:30
state ( 3 * nIntFaceTot + 3 , mappingHomogenization ( 1 , ip , el ) )
2018-08-30 14:20:52 +05:30
postResults ( c + 3 ) = homogState ( mappingHomogenization ( 2 , ip , el ) ) % &
2014-08-21 23:18:20 +05:30
state ( 3 * nIntFaceTot + 4 , mappingHomogenization ( 1 , ip , el ) )
2012-02-21 22:01:37 +05:30
c = c + 3_pInt
2013-12-18 12:58:01 +05:30
case ( penaltyenergy_ID )
2018-08-30 14:20:52 +05:30
postResults ( c + 1 ) = homogState ( mappingHomogenization ( 2 , ip , el ) ) % &
2014-08-21 23:18:20 +05:30
state ( 3 * nIntFaceTot + 5 , mappingHomogenization ( 1 , ip , el ) )
c = c + 1_pInt
2013-12-18 12:58:01 +05:30
case ( volumediscrepancy_ID )
2018-08-30 14:20:52 +05:30
postResults ( c + 1 ) = homogState ( mappingHomogenization ( 2 , ip , el ) ) % &
2014-08-21 23:18:20 +05:30
state ( 3 * nIntFaceTot + 6 , mappingHomogenization ( 1 , ip , el ) )
c = c + 1_pInt
2013-12-18 12:58:01 +05:30
case ( averagerelaxrate_ID )
2018-08-30 14:20:52 +05:30
postResults ( c + 1 ) = homogState ( mappingHomogenization ( 2 , ip , el ) ) % &
2014-08-21 23:18:20 +05:30
state ( 3 * nIntFaceTot + 7 , mappingHomogenization ( 1 , ip , el ) )
c = c + 1_pInt
2013-12-18 12:58:01 +05:30
case ( maximumrelaxrate_ID )
2018-08-30 14:20:52 +05:30
postResults ( c + 1 ) = homogState ( mappingHomogenization ( 2 , ip , el ) ) % &
2014-08-21 23:18:20 +05:30
state ( 3 * nIntFaceTot + 8 , mappingHomogenization ( 1 , ip , el ) )
c = c + 1_pInt
2009-10-12 21:31:49 +05:30
end select
2018-08-31 10:05:30 +05:30
enddo outputsLoop
2018-08-30 14:16:30 +05:30
end associate
2013-02-21 03:36:15 +05:30
end function homogenization_RGC_postResults
!--------------------------------------------------------------------------------------------------
!> @brief collect relaxation vectors of an interface
!--------------------------------------------------------------------------------------------------
2018-11-04 03:43:20 +05:30
function relaxationVector ( intFace , instance , of )
2013-02-21 03:36:15 +05:30
2009-10-12 21:31:49 +05:30
implicit none
2018-11-04 03:43:20 +05:30
integer ( pInt ) , intent ( in ) :: instance , of
2018-08-30 14:20:52 +05:30
real ( pReal ) , dimension ( 3 ) :: relaxationVector
2013-02-21 03:36:15 +05:30
integer ( pInt ) , dimension ( 4 ) , intent ( in ) :: intFace !< set of interface ID in 4D array (normal and position)
2014-02-08 16:18:09 +05:30
integer ( pInt ) :: &
2018-11-04 03:43:20 +05:30
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
2018-08-30 14:20:52 +05:30
relaxationVector = 0.0_pReal
iNum = interface4to1 ( intFace , instance ) ! identify the position of the interface in global state array
2018-11-04 03:43:20 +05:30
if ( iNum > 0_pInt ) relaxationVector = state ( instance ) % relaxationVector ( ( 3 * iNum - 2 ) : ( 3 * iNum ) , of )
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
!--------------------------------------------------------------------------------------------------
!> @brief identify the normal of an interface
!--------------------------------------------------------------------------------------------------
2018-11-04 03:43:20 +05:30
function interfaceNormal ( intFace , instance , of )
2013-02-21 03:36:15 +05:30
use math , only : &
math_mul33x3
2012-03-09 01:55:28 +05:30
2009-10-12 21:31:49 +05:30
implicit none
2018-08-30 14:20:52 +05:30
real ( pReal ) , dimension ( 3 ) :: interfaceNormal
2013-02-21 03:36:15 +05:30
integer ( pInt ) , dimension ( 4 ) , intent ( in ) :: intFace !< interface ID in 4D array (normal and position)
2014-02-08 16:18:09 +05:30
integer ( pInt ) , intent ( in ) :: &
2018-11-04 03:43:20 +05:30
instance , &
of
integer ( pInt ) :: 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)
2018-08-30 14:20:52 +05:30
interfaceNormal = 0.0_pReal
2012-02-16 00:28:38 +05:30
nPos = abs ( intFace ( 1 ) ) ! identify the position of the interface in global state array
2018-08-30 21:21:31 +05:30
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
2018-11-04 03:53:52 +05:30
interfaceNormal = math_mul33x3 ( 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-08-30 14:20:52 +05:30
function getInterface ( iFace , iGrain3 )
2013-02-21 03:36:15 +05:30
2009-10-12 21:31:49 +05:30
implicit none
2018-08-30 14:20:52 +05:30
integer ( pInt ) , dimension ( 4 ) :: getInterface
2014-02-08 16:18:09 +05:30
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
2009-10-12 21:31:49 +05:30
!* Direction of interface normal
2012-02-21 22:01:37 +05:30
iDir = ( int ( real ( iFace - 1_pInt , pReal ) / 2.0_pReal , pInt ) + 1_pInt ) * ( - 1_pInt ) ** iFace
2018-08-30 14:20:52 +05:30
getInterface ( 1 ) = iDir
2009-10-12 21:31:49 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! identify the interface position by the direction of its normal
2018-08-30 14:20:52 +05:30
getInterface ( 2 : 4 ) = iGrain3
2013-02-21 03:36:15 +05:30
if ( iDir < 0_pInt ) & ! to have a correlation with coordinate/position in real space
2018-08-30 14:20:52 +05:30
getInterface ( 1_pInt - iDir ) = getInterface ( 1_pInt - iDir ) - 1_pInt
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
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief map grain ID from in 1D (global array) to in 3D (local position)
!--------------------------------------------------------------------------------------------------
2018-08-30 14:20:52 +05:30
function grain1to3 ( grain1 , instance )
2012-03-09 01:55:28 +05:30
2009-10-12 21:31:49 +05:30
implicit none
2018-08-30 14:20:52 +05:30
integer ( pInt ) , dimension ( 3 ) :: grain1to3
2014-02-08 16:18:09 +05:30
integer ( pInt ) , intent ( in ) :: &
2013-02-21 03:36:15 +05:30
grain1 , & !< grain ID in 1D array
2018-08-30 14:16:30 +05:30
instance
2014-02-08 16:18:09 +05:30
integer ( pInt ) , dimension ( 3 ) :: nGDim
2009-10-12 21:31:49 +05:30
2018-08-30 14:16:30 +05:30
nGDim = param ( instance ) % Nconstituents
2018-08-30 21:21:31 +05:30
grain1to3 = 1_pInt + [ mod ( ( grain1 - 1_pInt ) , nGDim ( 1 ) ) , &
mod ( ( grain1 - 1_pInt ) / nGDim ( 1 ) , nGDim ( 2 ) ) , &
( grain1 - 1_pInt ) / ( 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)
!--------------------------------------------------------------------------------------------------
2018-08-30 14:20:52 +05:30
pure function grain3to1 ( grain3 , instance )
2009-10-12 21:31:49 +05:30
2012-03-09 01:55:28 +05:30
implicit none
2013-02-21 03:36:15 +05:30
integer ( pInt ) , dimension ( 3 ) , intent ( in ) :: grain3 !< grain ID in 3D array (pos.x,pos.y,pos.z)
2018-11-03 21:20:43 +05:30
integer ( pInt ) , intent ( in ) :: instance ! homogenization ID
2018-08-30 21:21:31 +05:30
integer ( pInt ) :: grain3to1
integer ( pInt ) , dimension ( 3 ) :: nGDim
2009-10-12 21:31:49 +05:30
2018-08-29 17:22:26 +05:30
nGDim = param ( instance ) % Nconstituents
2018-08-30 14:20:52 +05:30
grain3to1 = grain3 ( 1 ) + nGDim ( 1 ) * ( grain3 ( 2 ) - 1_pInt ) + nGDim ( 1 ) * nGDim ( 2 ) * ( grain3 ( 3 ) - 1_pInt )
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)
!--------------------------------------------------------------------------------------------------
2018-08-30 14:20:52 +05:30
integer ( pInt ) pure function interface4to1 ( iFace4D , instance )
2012-03-09 01:55:28 +05:30
2009-10-12 21:31:49 +05:30
implicit none
2013-02-21 03:36:15 +05:30
integer ( pInt ) , dimension ( 4 ) , intent ( in ) :: iFace4D !< interface ID in 4D array (n.dir,pos.x,pos.y,pos.z)
2018-08-30 21:21:31 +05:30
integer ( pInt ) , intent ( in ) :: instance
2009-10-12 21:31:49 +05:30
integer ( pInt ) , dimension ( 3 ) :: nGDim , nIntFace
2018-08-29 17:22:26 +05:30
nGDim = param ( instance ) % Nconstituents
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! 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
2010-11-26 17:20:20 +05:30
2018-08-30 14:20:52 +05:30
interface4to1 = - 1_pInt
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! get the corresponding interface ID in 1D global array
if ( abs ( iFace4D ( 1 ) ) == 1_pInt ) then ! interface with normal //e1
2018-08-30 14:20:52 +05:30
interface4to1 = iFace4D ( 3 ) + nGDim ( 2 ) * ( iFace4D ( 4 ) - 1_pInt ) &
2018-10-13 21:51:13 +05:30
+ nGDim ( 2 ) * nGDim ( 3 ) * ( iFace4D ( 2 ) - 1_pInt )
2018-08-30 14:20:52 +05:30
if ( ( iFace4D ( 2 ) == 0_pInt ) . or . ( iFace4D ( 2 ) == nGDim ( 1 ) ) ) interface4to1 = 0_pInt
2013-02-21 03:36:15 +05:30
elseif ( abs ( iFace4D ( 1 ) ) == 2_pInt ) then ! interface with normal //e2
2018-08-30 14:20:52 +05:30
interface4to1 = iFace4D ( 4 ) + nGDim ( 3 ) * ( iFace4D ( 2 ) - 1_pInt ) &
2018-10-13 21:51:13 +05:30
+ nGDim ( 3 ) * nGDim ( 1 ) * ( iFace4D ( 3 ) - 1_pInt ) + nIntFace ( 1 )
2018-08-30 14:20:52 +05:30
if ( ( iFace4D ( 3 ) == 0_pInt ) . or . ( iFace4D ( 3 ) == nGDim ( 2 ) ) ) interface4to1 = 0_pInt
2013-02-21 03:36:15 +05:30
elseif ( abs ( iFace4D ( 1 ) ) == 3_pInt ) then ! interface with normal //e3
2018-08-30 14:20:52 +05:30
interface4to1 = iFace4D ( 2 ) + nGDim ( 1 ) * ( iFace4D ( 3 ) - 1_pInt ) &
2018-10-13 21:51:13 +05:30
+ nGDim ( 1 ) * nGDim ( 2 ) * ( iFace4D ( 4 ) - 1_pInt ) + nIntFace ( 1 ) + nIntFace ( 2 )
2018-08-30 14:20:52 +05:30
if ( ( iFace4D ( 4 ) == 0_pInt ) . or . ( iFace4D ( 4 ) == nGDim ( 3 ) ) ) interface4to1 = 0_pInt
2009-10-12 21:31:49 +05:30
endif
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-08-30 14:20:52 +05:30
pure function interface1to4 ( iFace1D , instance )
2012-03-09 01:55:28 +05:30
2009-10-12 21:31:49 +05:30
implicit none
2018-08-30 21:21:31 +05:30
integer ( pInt ) , dimension ( 4 ) :: interface1to4
integer ( pInt ) , intent ( in ) :: iFace1D !< interface ID in 1D array
integer ( pInt ) , intent ( in ) :: instance
integer ( pInt ) , dimension ( 3 ) :: nGDim , nIntFace
2009-10-12 21:31:49 +05:30
2018-08-29 17:22:26 +05:30
nGDim = param ( instance ) % Nconstituents
2009-10-12 21:31:49 +05:30
2013-02-21 03:36:15 +05:30
!--------------------------------------------------------------------------------------------------
! 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
2018-08-30 14:20:52 +05:30
interface1to4 ( 1 ) = 1_pInt
interface1to4 ( 3 ) = mod ( ( iFace1D - 1_pInt ) , nGDim ( 2 ) ) + 1_pInt
interface1to4 ( 4 ) = mod ( &
2012-02-21 22:01:37 +05:30
int ( &
real ( iFace1D - 1_pInt , pReal ) / &
real ( nGDim ( 2 ) , pReal ) &
, pInt ) &
, nGDim ( 3 ) ) + 1_pInt
2018-08-30 14:20:52 +05:30
interface1to4 ( 2 ) = int ( &
2012-02-21 22:01:37 +05:30
real ( iFace1D - 1_pInt , pReal ) / &
real ( nGDim ( 2 ) , pReal ) / &
real ( nGDim ( 3 ) , pReal ) &
, pInt ) + 1_pInt
2013-02-21 03:36:15 +05:30
elseif ( iFace1D > nIntFace ( 1 ) . and . iFace1D < = ( nIntFace ( 2 ) + nIntFace ( 1 ) ) ) then ! interface with normal //e2
2018-08-30 14:20:52 +05:30
interface1to4 ( 1 ) = 2_pInt
interface1to4 ( 4 ) = mod ( ( iFace1D - nIntFace ( 1 ) - 1_pInt ) , nGDim ( 3 ) ) + 1_pInt
interface1to4 ( 2 ) = mod ( &
2012-02-21 22:01:37 +05:30
int ( &
real ( iFace1D - nIntFace ( 1 ) - 1_pInt , pReal ) / &
real ( nGDim ( 3 ) , pReal ) &
, pInt ) &
, nGDim ( 1 ) ) + 1_pInt
2018-08-30 14:20:52 +05:30
interface1to4 ( 3 ) = int ( &
2012-02-21 22:01:37 +05:30
real ( iFace1D - nIntFace ( 1 ) - 1_pInt , pReal ) / &
real ( nGDim ( 3 ) , pReal ) / &
real ( nGDim ( 1 ) , pReal ) &
, pInt ) + 1_pInt
2013-02-21 03:36:15 +05:30
elseif ( iFace1D > nIntFace ( 2 ) + nIntFace ( 1 ) . and . iFace1D < = ( nIntFace ( 3 ) + nIntFace ( 2 ) + nIntFace ( 1 ) ) ) then ! interface with normal //e3
2018-08-30 14:20:52 +05:30
interface1to4 ( 1 ) = 3_pInt
interface1to4 ( 2 ) = mod ( ( iFace1D - nIntFace ( 2 ) - nIntFace ( 1 ) - 1_pInt ) , nGDim ( 1 ) ) + 1_pInt
interface1to4 ( 3 ) = mod ( &
2012-02-21 22:01:37 +05:30
int ( &
real ( iFace1D - nIntFace ( 2 ) - nIntFace ( 1 ) - 1_pInt , pReal ) / &
real ( nGDim ( 1 ) , pReal ) &
, pInt ) &
, nGDim ( 2 ) ) + 1_pInt
2018-08-30 14:20:52 +05:30
interface1to4 ( 4 ) = int ( &
2012-02-21 22:01:37 +05:30
real ( iFace1D - nIntFace ( 2 ) - nIntFace ( 1 ) - 1_pInt , pReal ) / &
real ( nGDim ( 1 ) , pReal ) / &
real ( nGDim ( 2 ) , pReal ) &
, pInt ) + 1_pInt
2009-10-12 21:31:49 +05:30
endif
2018-08-30 14:20:52 +05:30
end function interface1to4
2013-02-21 03:36:15 +05:30
end module homogenization_RGC