2020-07-15 18:05:21 +05:30
!----------------------------------------------------------------------------------------------------
2020-09-14 00:30:34 +05:30
!> @brief internal microstructure state for all plasticity constitutive models
2020-07-15 18:05:21 +05:30
!----------------------------------------------------------------------------------------------------
2020-11-03 05:19:17 +05:30
submodule ( constitutive ) constitutive_mech
2020-07-09 04:31:08 +05:30
2020-12-15 22:15:11 +05:30
integer ( kind ( ELASTICITY_undefined_ID ) ) , dimension ( : ) , allocatable :: &
phase_elasticity !< elasticity of each phase
integer ( kind ( SOURCE_undefined_ID ) ) , dimension ( : , : ) , allocatable :: &
phase_stiffnessDegradation !< active stiffness degradation mechanisms of each phase
2020-07-09 04:31:08 +05:30
interface
2020-08-15 19:32:10 +05:30
module function plastic_none_init ( ) result ( myPlasticity )
logical , dimension ( : ) , allocatable :: &
myPlasticity
end function plastic_none_init
module function plastic_isotropic_init ( ) result ( myPlasticity )
logical , dimension ( : ) , allocatable :: &
myPlasticity
end function plastic_isotropic_init
module function plastic_phenopowerlaw_init ( ) result ( myPlasticity )
logical , dimension ( : ) , allocatable :: &
myPlasticity
end function plastic_phenopowerlaw_init
module function plastic_kinehardening_init ( ) result ( myPlasticity )
logical , dimension ( : ) , allocatable :: &
myPlasticity
end function plastic_kinehardening_init
module function plastic_dislotwin_init ( ) result ( myPlasticity )
logical , dimension ( : ) , allocatable :: &
myPlasticity
end function plastic_dislotwin_init
2020-11-29 16:28:39 +05:30
module function plastic_dislotungsten_init ( ) result ( myPlasticity )
2020-08-15 19:32:10 +05:30
logical , dimension ( : ) , allocatable :: &
myPlasticity
2020-11-29 16:28:39 +05:30
end function plastic_dislotungsten_init
2020-08-15 19:32:10 +05:30
module function plastic_nonlocal_init ( ) result ( myPlasticity )
logical , dimension ( : ) , allocatable :: &
myPlasticity
end function plastic_nonlocal_init
2020-07-09 04:31:08 +05:30
module subroutine plastic_isotropic_LpAndItsTangent ( Lp , dLp_dMp , Mp , instance , of )
real ( pReal ) , dimension ( 3 , 3 ) , intent ( out ) :: &
Lp !< plastic velocity gradient
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) , intent ( out ) :: &
dLp_dMp !< derivative of Lp with respect to the Mandel stress
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: &
Mp !< Mandel stress
integer , intent ( in ) :: &
instance , &
of
end subroutine plastic_isotropic_LpAndItsTangent
pure module subroutine plastic_phenopowerlaw_LpAndItsTangent ( Lp , dLp_dMp , Mp , instance , of )
real ( pReal ) , dimension ( 3 , 3 ) , intent ( out ) :: &
Lp !< plastic velocity gradient
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) , intent ( out ) :: &
dLp_dMp !< derivative of Lp with respect to the Mandel stress
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: &
Mp !< Mandel stress
integer , intent ( in ) :: &
instance , &
of
end subroutine plastic_phenopowerlaw_LpAndItsTangent
pure module subroutine plastic_kinehardening_LpAndItsTangent ( Lp , dLp_dMp , Mp , instance , of )
real ( pReal ) , dimension ( 3 , 3 ) , intent ( out ) :: &
Lp !< plastic velocity gradient
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) , intent ( out ) :: &
dLp_dMp !< derivative of Lp with respect to the Mandel stress
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: &
Mp !< Mandel stress
integer , intent ( in ) :: &
instance , &
of
end subroutine plastic_kinehardening_LpAndItsTangent
module subroutine plastic_dislotwin_LpAndItsTangent ( Lp , dLp_dMp , Mp , T , instance , of )
real ( pReal ) , dimension ( 3 , 3 ) , intent ( out ) :: &
Lp !< plastic velocity gradient
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) , intent ( out ) :: &
dLp_dMp !< derivative of Lp with respect to the Mandel stress
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: &
Mp !< Mandel stress
real ( pReal ) , intent ( in ) :: &
T
integer , intent ( in ) :: &
instance , &
of
end subroutine plastic_dislotwin_LpAndItsTangent
2020-11-29 16:28:39 +05:30
pure module subroutine plastic_dislotungsten_LpAndItsTangent ( Lp , dLp_dMp , Mp , T , instance , of )
2020-07-09 04:31:08 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( out ) :: &
Lp !< plastic velocity gradient
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) , intent ( out ) :: &
dLp_dMp !< derivative of Lp with respect to the Mandel stress
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: &
Mp !< Mandel stress
real ( pReal ) , intent ( in ) :: &
T
integer , intent ( in ) :: &
instance , &
of
2020-11-29 16:28:39 +05:30
end subroutine plastic_dislotungsten_LpAndItsTangent
2020-07-09 04:31:08 +05:30
module subroutine plastic_nonlocal_LpAndItsTangent ( Lp , dLp_dMp , &
Mp , Temperature , instance , of , ip , el )
real ( pReal ) , dimension ( 3 , 3 ) , intent ( out ) :: &
Lp !< plastic velocity gradient
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) , intent ( out ) :: &
dLp_dMp !< derivative of Lp with respect to the Mandel stress
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: &
Mp !< Mandel stress
real ( pReal ) , intent ( in ) :: &
Temperature
integer , intent ( in ) :: &
instance , &
of , &
ip , & !< current integration point
el !< current element number
end subroutine plastic_nonlocal_LpAndItsTangent
2020-12-22 23:54:00 +05:30
module subroutine plastic_isotropic_dotState ( Mp , instance , of )
2020-12-20 21:02:33 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: &
Mp !< Mandel stress
integer , intent ( in ) :: &
instance , &
of
end subroutine plastic_isotropic_dotState
module subroutine plastic_phenopowerlaw_dotState ( Mp , instance , of )
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: &
Mp !< Mandel stress
integer , intent ( in ) :: &
instance , &
of
end subroutine plastic_phenopowerlaw_dotState
module subroutine plastic_kinehardening_dotState ( Mp , instance , of )
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: &
Mp !< Mandel stress
integer , intent ( in ) :: &
instance , &
of
end subroutine plastic_kinehardening_dotState
module subroutine plastic_dislotwin_dotState ( Mp , T , instance , of )
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: &
Mp !< Mandel stress
real ( pReal ) , intent ( in ) :: &
T
integer , intent ( in ) :: &
instance , &
of
end subroutine plastic_dislotwin_dotState
module subroutine plastic_disloTungsten_dotState ( Mp , T , instance , of )
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: &
Mp !< Mandel stress
real ( pReal ) , intent ( in ) :: &
T
integer , intent ( in ) :: &
instance , &
of
end subroutine plastic_disloTungsten_dotState
2020-12-23 01:15:27 +05:30
module subroutine plastic_nonlocal_dotState ( Mp , F , Temperature , timestep , &
2020-12-20 21:02:33 +05:30
instance , of , ip , el )
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: &
Mp !< MandelStress
real ( pReal ) , dimension ( 3 , 3 , homogenization_maxNconstituents , discretization_nIPs , discretization_Nelems ) , intent ( in ) :: &
2020-12-23 01:15:27 +05:30
F !< deformation gradient
2020-12-20 21:02:33 +05:30
real ( pReal ) , intent ( in ) :: &
Temperature , & !< temperature
timestep !< substepped crystallite time increment
integer , intent ( in ) :: &
instance , &
of , &
ip , & !< current integration point
el !< current element number
end subroutine plastic_nonlocal_dotState
2020-07-09 04:31:08 +05:30
module subroutine plastic_dislotwin_dependentState ( T , instance , of )
integer , intent ( in ) :: &
instance , &
of
real ( pReal ) , intent ( in ) :: &
T
end subroutine plastic_dislotwin_dependentState
2020-11-29 16:28:39 +05:30
module subroutine plastic_dislotungsten_dependentState ( instance , of )
2020-07-09 04:31:08 +05:30
integer , intent ( in ) :: &
instance , &
of
2020-11-29 16:28:39 +05:30
end subroutine plastic_dislotungsten_dependentState
2020-07-09 04:31:08 +05:30
2020-12-23 01:15:27 +05:30
module subroutine plastic_nonlocal_dependentState ( F , instance , of , ip , el )
2020-07-09 04:31:08 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: &
2020-12-23 01:15:27 +05:30
F !< deformation gradient
2020-07-09 04:31:08 +05:30
integer , intent ( in ) :: &
instance , &
of , &
2020-07-15 18:05:21 +05:30
ip , & !< current integration point
el !< current element number
2020-07-09 04:31:08 +05:30
end subroutine plastic_nonlocal_dependentState
2020-12-22 23:54:00 +05:30
module subroutine plastic_kinehardening_deltaState ( Mp , instance , of )
2020-12-20 21:41:43 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: &
Mp !< Mandel stress
integer , intent ( in ) :: &
instance , &
of
end subroutine plastic_kinehardening_deltaState
module subroutine plastic_nonlocal_deltaState ( Mp , instance , of , ip , el )
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: &
Mp
integer , intent ( in ) :: &
instance , &
of , &
ip , &
el
end subroutine plastic_nonlocal_deltaState
2020-07-12 18:52:40 +05:30
module subroutine plastic_isotropic_results ( instance , group )
integer , intent ( in ) :: instance
character ( len = * ) , intent ( in ) :: group
end subroutine plastic_isotropic_results
module subroutine plastic_phenopowerlaw_results ( instance , group )
integer , intent ( in ) :: instance
character ( len = * ) , intent ( in ) :: group
end subroutine plastic_phenopowerlaw_results
module subroutine plastic_kinehardening_results ( instance , group )
integer , intent ( in ) :: instance
character ( len = * ) , intent ( in ) :: group
end subroutine plastic_kinehardening_results
module subroutine plastic_dislotwin_results ( instance , group )
integer , intent ( in ) :: instance
character ( len = * ) , intent ( in ) :: group
end subroutine plastic_dislotwin_results
2020-11-29 16:28:39 +05:30
module subroutine plastic_dislotungsten_results ( instance , group )
2020-07-12 18:52:40 +05:30
integer , intent ( in ) :: instance
character ( len = * ) , intent ( in ) :: group
2020-11-29 16:28:39 +05:30
end subroutine plastic_dislotungsten_results
2020-07-12 18:52:40 +05:30
module subroutine plastic_nonlocal_results ( instance , group )
integer , intent ( in ) :: instance
character ( len = * ) , intent ( in ) :: group
end subroutine plastic_nonlocal_results
2020-07-09 04:31:08 +05:30
end interface
2020-12-22 15:14:43 +05:30
type :: tOutput !< new requested output (per phase)
character ( len = pStringLen ) , allocatable , dimension ( : ) :: &
label
end type tOutput
type ( tOutput ) , allocatable , dimension ( : ) :: output_constituent
2020-07-09 04:31:08 +05:30
contains
!--------------------------------------------------------------------------------------------------
2020-11-05 21:19:59 +05:30
!> @brief Initialize mechanical field related constitutive models
!> @details Initialize elasticity, plasticity and stiffness degradation models.
2020-07-09 04:31:08 +05:30
!--------------------------------------------------------------------------------------------------
2020-11-03 02:09:33 +05:30
module subroutine mech_init
2020-07-09 04:31:08 +05:30
2020-11-03 02:09:33 +05:30
integer :: &
2020-12-23 11:44:07 +05:30
ph , &
2020-11-03 02:09:33 +05:30
stiffDegradationCtr
class ( tNode ) , pointer :: &
2020-12-22 23:54:00 +05:30
num_crystallite , &
2020-11-03 02:09:33 +05:30
phases , &
phase , &
2020-11-03 03:16:46 +05:30
mech , &
2020-11-03 02:09:33 +05:30
elastic , &
stiffDegradation
2020-08-15 19:32:10 +05:30
2020-11-03 05:19:17 +05:30
print '(/,a)' , ' <<<+- constitutive_mech init -+>>>'
2020-09-14 00:30:34 +05:30
2020-11-03 02:09:33 +05:30
!-------------------------------------------------------------------------------------------------
! initialize elasticity (hooke) !ToDO: Maybe move to elastic submodule along with function homogenizedC?
2020-09-13 14:09:17 +05:30
phases = > config_material % get ( 'phase' )
2020-11-03 02:09:33 +05:30
allocate ( phase_elasticity ( phases % length ) , source = ELASTICITY_undefined_ID )
allocate ( phase_elasticityInstance ( phases % length ) , source = 0 )
allocate ( phase_NstiffnessDegradations ( phases % length ) , source = 0 )
2020-12-22 13:24:32 +05:30
allocate ( output_constituent ( phases % length ) )
2020-11-03 02:09:33 +05:30
2020-12-23 11:44:07 +05:30
do ph = 1 , phases % length
phase = > phases % get ( ph )
2020-11-18 01:54:40 +05:30
mech = > phase % get ( 'mechanics' )
2020-12-22 13:24:32 +05:30
#if defined(__GFORTRAN__)
2020-12-23 11:44:07 +05:30
output_constituent ( ph ) % label = output_asStrings ( mech )
2020-12-22 13:24:32 +05:30
#else
2020-12-23 11:44:07 +05:30
output_constituent ( ph ) % label = mech % get_asStrings ( 'output' , defaultVal = emptyStringArray )
2020-12-22 13:24:32 +05:30
#endif
2020-11-03 03:16:46 +05:30
elastic = > mech % get ( 'elasticity' )
2020-11-03 02:09:33 +05:30
if ( elastic % get_asString ( 'type' ) == 'hooke' ) then
2020-12-23 11:44:07 +05:30
phase_elasticity ( ph ) = ELASTICITY_HOOKE_ID
2020-11-03 02:09:33 +05:30
else
call IO_error ( 200 , ext_msg = elastic % get_asString ( 'type' ) )
endif
2020-11-05 21:19:59 +05:30
stiffDegradation = > mech % get ( 'stiffness_degradation' , defaultVal = emptyList ) ! check for stiffness degradation mechanisms
2020-12-23 11:44:07 +05:30
phase_NstiffnessDegradations ( ph ) = stiffDegradation % length
2020-11-03 02:09:33 +05:30
enddo
allocate ( phase_stiffnessDegradation ( maxval ( phase_NstiffnessDegradations ) , phases % length ) , &
source = STIFFNESS_DEGRADATION_undefined_ID )
if ( maxVal ( phase_NstiffnessDegradations ) / = 0 ) then
2020-12-23 11:44:07 +05:30
do ph = 1 , phases % length
phase = > phases % get ( ph )
2020-11-18 01:54:40 +05:30
mech = > phase % get ( 'mechanics' )
2020-11-03 03:16:46 +05:30
stiffDegradation = > mech % get ( 'stiffness_degradation' , defaultVal = emptyList )
2020-11-03 02:09:33 +05:30
do stiffDegradationCtr = 1 , stiffDegradation % length
if ( stiffDegradation % get_asString ( stiffDegradationCtr ) == 'damage' ) &
2020-12-23 11:44:07 +05:30
phase_stiffnessDegradation ( stiffDegradationCtr , ph ) = STIFFNESS_DEGRADATION_damage_ID
2020-11-03 02:09:33 +05:30
enddo
enddo
endif
2020-08-15 19:32:10 +05:30
2020-11-03 05:19:17 +05:30
! initialize plasticity
2020-08-15 19:32:10 +05:30
allocate ( plasticState ( phases % length ) )
allocate ( phase_plasticity ( phases % length ) , source = PLASTICITY_undefined_ID )
allocate ( phase_plasticityInstance ( phases % length ) , source = 0 )
allocate ( phase_localPlasticity ( phases % length ) , source = . true . )
where ( plastic_none_init ( ) ) phase_plasticity = PLASTICITY_NONE_ID
where ( plastic_isotropic_init ( ) ) phase_plasticity = PLASTICITY_ISOTROPIC_ID
where ( plastic_phenopowerlaw_init ( ) ) phase_plasticity = PLASTICITY_PHENOPOWERLAW_ID
where ( plastic_kinehardening_init ( ) ) phase_plasticity = PLASTICITY_KINEHARDENING_ID
where ( plastic_dislotwin_init ( ) ) phase_plasticity = PLASTICITY_DISLOTWIN_ID
2020-11-29 16:28:39 +05:30
where ( plastic_dislotungsten_init ( ) ) phase_plasticity = PLASTICITY_DISLOTUNGSTEN_ID
2020-08-15 19:32:10 +05:30
where ( plastic_nonlocal_init ( ) ) phase_plasticity = PLASTICITY_NONLOCAL_ID
2020-12-23 11:44:07 +05:30
do ph = 1 , phases % length
phase_elasticityInstance ( ph ) = count ( phase_elasticity ( 1 : ph ) == phase_elasticity ( ph ) )
phase_plasticityInstance ( ph ) = count ( phase_plasticity ( 1 : ph ) == phase_plasticity ( ph ) )
2020-09-14 00:30:34 +05:30
enddo
2020-08-15 19:32:10 +05:30
2020-12-22 23:54:00 +05:30
num_crystallite = > config_numerics % get ( 'crystallite' , defaultVal = emptyDict )
select case ( num_crystallite % get_asString ( 'integrator' , defaultVal = 'FPI' ) )
case ( 'FPI' )
integrateState = > integrateStateFPI
case ( 'Euler' )
integrateState = > integrateStateEuler
case ( 'AdaptiveEuler' )
integrateState = > integrateStateAdaptiveEuler
case ( 'RK4' )
integrateState = > integrateStateRK4
case ( 'RKCK45' )
integrateState = > integrateStateRKCK45
case default
call IO_error ( 301 , ext_msg = 'integrator' )
end select
2020-11-03 02:09:33 +05:30
end subroutine mech_init
2020-07-09 04:31:08 +05:30
2020-07-10 18:29:07 +05:30
2020-08-15 19:32:10 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief checks if a plastic module is active or not
!--------------------------------------------------------------------------------------------------
2020-12-22 14:33:19 +05:30
function plastic_active ( plastic_label ) result ( active_plastic )
2020-08-15 19:32:10 +05:30
character ( len = * ) , intent ( in ) :: plastic_label !< type of plasticity model
logical , dimension ( : ) , allocatable :: active_plastic
class ( tNode ) , pointer :: &
phases , &
phase , &
2020-11-03 03:16:46 +05:30
mech , &
2020-08-15 19:32:10 +05:30
pl
2020-12-23 11:44:07 +05:30
integer :: ph
2020-08-15 19:32:10 +05:30
2020-09-13 14:09:17 +05:30
phases = > config_material % get ( 'phase' )
2020-08-15 19:32:10 +05:30
allocate ( active_plastic ( phases % length ) , source = . false . )
2020-12-23 11:44:07 +05:30
do ph = 1 , phases % length
phase = > phases % get ( ph )
2020-11-18 01:54:40 +05:30
mech = > phase % get ( 'mechanics' )
2020-11-03 03:16:46 +05:30
pl = > mech % get ( 'plasticity' )
2020-12-23 11:44:07 +05:30
if ( pl % get_asString ( 'type' ) == plastic_label ) active_plastic ( ph ) = . true .
2020-08-15 19:32:10 +05:30
enddo
end function plastic_active
2020-11-03 02:09:33 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief returns the 2nd Piola-Kirchhoff stress tensor and its tangent with respect to
!> the elastic and intermediate deformation gradients using Hooke's law
!--------------------------------------------------------------------------------------------------
2020-11-04 17:13:52 +05:30
module subroutine constitutive_hooke_SandItsTangents ( S , dS_dFe , dS_dFi , &
2020-12-23 11:37:18 +05:30
Fe , Fi , co , ip , el )
2020-11-03 02:09:33 +05:30
integer , intent ( in ) :: &
2020-12-23 11:37:18 +05:30
co , & !< component-ID of integration point
2020-11-03 02:09:33 +05:30
ip , & !< integration point
el !< element
real ( pReal ) , intent ( in ) , dimension ( 3 , 3 ) :: &
Fe , & !< elastic deformation gradient
Fi !< intermediate deformation gradient
real ( pReal ) , intent ( out ) , dimension ( 3 , 3 ) :: &
S !< 2nd Piola-Kirchhoff stress tensor in lattice configuration
real ( pReal ) , intent ( out ) , dimension ( 3 , 3 , 3 , 3 ) :: &
dS_dFe , & !< derivative of 2nd P-K stress with respect to elastic deformation gradient
dS_dFi !< derivative of 2nd P-K stress with respect to intermediate deformation gradient
real ( pReal ) , dimension ( 3 , 3 ) :: E
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) :: C
integer :: &
ho , & !< homogenization
d !< counter in degradation loop
integer :: &
i , j
ho = material_homogenizationAt ( el )
2020-12-23 11:37:18 +05:30
C = math_66toSym3333 ( constitutive_homogenizedC ( co , ip , el ) )
2020-11-03 02:09:33 +05:30
2020-12-23 11:37:18 +05:30
DegradationLoop : do d = 1 , phase_NstiffnessDegradations ( material_phaseAt ( co , el ) )
degradationType : select case ( phase_stiffnessDegradation ( d , material_phaseAt ( co , el ) ) )
2020-11-03 02:09:33 +05:30
case ( STIFFNESS_DEGRADATION_damage_ID ) degradationType
2020-12-16 00:25:55 +05:30
C = C * damage ( ho ) % p ( material_homogenizationMemberAt ( ip , el ) ) ** 2
2020-11-03 02:09:33 +05:30
end select degradationType
enddo DegradationLoop
E = 0.5_pReal * ( matmul ( transpose ( Fe ) , Fe ) - math_I3 ) !< Green-Lagrange strain in unloaded configuration
S = math_mul3333xx33 ( C , matmul ( matmul ( transpose ( Fi ) , E ) , Fi ) ) !< 2PK stress in lattice configuration in work conjugate with GL strain pulled back to lattice configuration
do i = 1 , 3 ; do j = 1 , 3
dS_dFe ( i , j , 1 : 3 , 1 : 3 ) = matmul ( Fe , matmul ( matmul ( Fi , C ( i , j , 1 : 3 , 1 : 3 ) ) , transpose ( Fi ) ) ) !< dS_ij/dFe_kl = C_ijmn * Fi_lm * Fi_on * Fe_ko
dS_dFi ( i , j , 1 : 3 , 1 : 3 ) = 2.0_pReal * matmul ( matmul ( E , Fi ) , C ( i , j , 1 : 3 , 1 : 3 ) ) !< dS_ij/dFi_kl = C_ijln * E_km * Fe_mn
enddo ; enddo
end subroutine constitutive_hooke_SandItsTangents
2020-07-09 04:31:08 +05:30
!--------------------------------------------------------------------------------------------------
2020-07-12 20:14:26 +05:30
!> @brief calls microstructure function of the different plasticity constitutive models
2020-07-09 04:31:08 +05:30
!--------------------------------------------------------------------------------------------------
2020-12-23 11:37:18 +05:30
module subroutine constitutive_plastic_dependentState ( F , co , ip , el )
2020-07-12 20:14:26 +05:30
integer , intent ( in ) :: &
2020-12-23 14:29:47 +05:30
co , & !< component-ID of integration point
2020-07-12 20:14:26 +05:30
ip , & !< integration point
el !< element
real ( pReal ) , intent ( in ) , dimension ( 3 , 3 ) :: &
2020-12-23 14:29:47 +05:30
F !< deformation gradient
2020-07-10 18:43:56 +05:30
2020-07-09 04:31:08 +05:30
integer :: &
ho , & !< homogenization
tme , & !< thermal member position
instance , of
ho = material_homogenizationAt ( el )
2020-12-15 22:15:11 +05:30
tme = material_homogenizationMemberAt ( ip , el )
2020-12-23 11:37:18 +05:30
of = material_phasememberAt ( co , ip , el )
instance = phase_plasticityInstance ( material_phaseAt ( co , el ) )
2020-07-09 04:31:08 +05:30
2020-12-23 11:37:18 +05:30
plasticityType : select case ( phase_plasticity ( material_phaseAt ( co , el ) ) )
2020-07-09 04:31:08 +05:30
case ( PLASTICITY_DISLOTWIN_ID ) plasticityType
call plastic_dislotwin_dependentState ( temperature ( ho ) % p ( tme ) , instance , of )
2020-08-15 19:32:10 +05:30
case ( PLASTICITY_DISLOTUNGSTEN_ID ) plasticityType
2020-11-29 16:28:39 +05:30
call plastic_dislotungsten_dependentState ( instance , of )
2020-07-09 04:31:08 +05:30
case ( PLASTICITY_NONLOCAL_ID ) plasticityType
2020-12-23 01:15:27 +05:30
call plastic_nonlocal_dependentState ( F , instance , of , ip , el )
2020-07-09 04:31:08 +05:30
end select plasticityType
2020-07-12 20:14:26 +05:30
end subroutine constitutive_plastic_dependentState
2020-07-09 04:31:08 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the velocity gradient
! ToDo: Discuss whether it makes sense if crystallite handles the configuration conversion, i.e.
! Mp in, dLp_dMp out
!--------------------------------------------------------------------------------------------------
2020-07-12 20:14:26 +05:30
module subroutine constitutive_plastic_LpAndItsTangents ( Lp , dLp_dS , dLp_dFi , &
2020-12-23 11:37:18 +05:30
S , Fi , co , ip , el )
2020-07-12 20:14:26 +05:30
integer , intent ( in ) :: &
2020-12-23 14:29:47 +05:30
co , & !< component-ID of integration point
2020-07-12 20:14:26 +05:30
ip , & !< integration point
el !< element
real ( pReal ) , intent ( in ) , dimension ( 3 , 3 ) :: &
S , & !< 2nd Piola-Kirchhoff stress
Fi !< intermediate deformation gradient
real ( pReal ) , intent ( out ) , dimension ( 3 , 3 ) :: &
Lp !< plastic velocity gradient
real ( pReal ) , intent ( out ) , dimension ( 3 , 3 , 3 , 3 ) :: &
dLp_dS , &
dLp_dFi !< derivative of Lp with respect to Fi
2020-07-10 18:43:56 +05:30
2020-07-09 04:31:08 +05:30
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) :: &
dLp_dMp !< derivative of Lp with respect to Mandel stress
real ( pReal ) , dimension ( 3 , 3 ) :: &
Mp !< Mandel stress work conjugate with Lp
integer :: &
ho , & !< homogenization
tme !< thermal member position
integer :: &
i , j , instance , of
ho = material_homogenizationAt ( el )
2020-12-15 22:15:11 +05:30
tme = material_homogenizationMemberAt ( ip , el )
2020-07-09 04:31:08 +05:30
Mp = matmul ( matmul ( transpose ( Fi ) , Fi ) , S )
2020-12-23 11:37:18 +05:30
of = material_phasememberAt ( co , ip , el )
instance = phase_plasticityInstance ( material_phaseAt ( co , el ) )
2020-07-09 04:31:08 +05:30
2020-12-23 11:37:18 +05:30
plasticityType : select case ( phase_plasticity ( material_phaseAt ( co , el ) ) )
2020-07-09 04:31:08 +05:30
case ( PLASTICITY_NONE_ID ) plasticityType
Lp = 0.0_pReal
dLp_dMp = 0.0_pReal
case ( PLASTICITY_ISOTROPIC_ID ) plasticityType
2020-08-15 19:32:10 +05:30
call plastic_isotropic_LpAndItsTangent ( Lp , dLp_dMp , Mp , instance , of )
2020-07-09 04:31:08 +05:30
case ( PLASTICITY_PHENOPOWERLAW_ID ) plasticityType
call plastic_phenopowerlaw_LpAndItsTangent ( Lp , dLp_dMp , Mp , instance , of )
case ( PLASTICITY_KINEHARDENING_ID ) plasticityType
call plastic_kinehardening_LpAndItsTangent ( Lp , dLp_dMp , Mp , instance , of )
case ( PLASTICITY_NONLOCAL_ID ) plasticityType
2020-08-15 19:32:10 +05:30
call plastic_nonlocal_LpAndItsTangent ( Lp , dLp_dMp , Mp , temperature ( ho ) % p ( tme ) , instance , of , ip , el )
2020-07-09 04:31:08 +05:30
case ( PLASTICITY_DISLOTWIN_ID ) plasticityType
2020-08-15 19:32:10 +05:30
call plastic_dislotwin_LpAndItsTangent ( Lp , dLp_dMp , Mp , temperature ( ho ) % p ( tme ) , instance , of )
2020-07-09 04:31:08 +05:30
2020-08-15 19:32:10 +05:30
case ( PLASTICITY_DISLOTUNGSTEN_ID ) plasticityType
2020-11-29 16:28:39 +05:30
call plastic_dislotungsten_LpAndItsTangent ( Lp , dLp_dMp , Mp , temperature ( ho ) % p ( tme ) , instance , of )
2020-07-09 04:31:08 +05:30
end select plasticityType
do i = 1 , 3 ; do j = 1 , 3
dLp_dFi ( i , j , 1 : 3 , 1 : 3 ) = matmul ( matmul ( Fi , S ) , transpose ( dLp_dMp ( i , j , 1 : 3 , 1 : 3 ) ) ) + &
matmul ( matmul ( Fi , dLp_dMp ( i , j , 1 : 3 , 1 : 3 ) ) , S )
dLp_dS ( i , j , 1 : 3 , 1 : 3 ) = matmul ( matmul ( transpose ( Fi ) , Fi ) , dLp_dMp ( i , j , 1 : 3 , 1 : 3 ) ) ! ToDo: @PS: why not: dLp_dMp:(FiT Fi)
enddo ; enddo
2020-07-12 20:14:26 +05:30
end subroutine constitutive_plastic_LpAndItsTangents
2020-07-09 04:31:08 +05:30
2020-12-20 21:02:33 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief contains the constitutive equation for calculating the rate of change of microstructure
!--------------------------------------------------------------------------------------------------
2020-12-23 11:44:07 +05:30
function mech_collectDotState ( subdt , co , ip , el , ph , of ) result ( broken )
2020-12-20 21:02:33 +05:30
integer , intent ( in ) :: &
2020-12-23 11:37:18 +05:30
co , & !< component-ID of integration point
2020-12-20 21:02:33 +05:30
ip , & !< integration point
el , & !< element
2020-12-23 11:44:07 +05:30
ph , &
2020-12-20 21:02:33 +05:30
of
real ( pReal ) , intent ( in ) :: &
subdt !< timestep
real ( pReal ) , dimension ( 3 , 3 ) :: &
Mp
integer :: &
ho , & !< homogenization
tme , & !< thermal member position
i , & !< counter in source loop
instance
logical :: broken
ho = material_homogenizationAt ( el )
tme = material_homogenizationMemberAt ( ip , el )
2020-12-23 11:44:07 +05:30
instance = phase_plasticityInstance ( ph )
2020-12-20 21:02:33 +05:30
2020-12-23 11:44:07 +05:30
Mp = matmul ( matmul ( transpose ( constitutive_mech_Fi ( ph ) % data ( 1 : 3 , 1 : 3 , of ) ) , &
constitutive_mech_Fi ( ph ) % data ( 1 : 3 , 1 : 3 , of ) ) , crystallite_S ( 1 : 3 , 1 : 3 , co , ip , el ) )
2020-12-20 21:02:33 +05:30
2020-12-23 11:44:07 +05:30
plasticityType : select case ( phase_plasticity ( ph ) )
2020-12-20 21:02:33 +05:30
case ( PLASTICITY_ISOTROPIC_ID ) plasticityType
call plastic_isotropic_dotState ( Mp , instance , of )
case ( PLASTICITY_PHENOPOWERLAW_ID ) plasticityType
call plastic_phenopowerlaw_dotState ( Mp , instance , of )
case ( PLASTICITY_KINEHARDENING_ID ) plasticityType
call plastic_kinehardening_dotState ( Mp , instance , of )
case ( PLASTICITY_DISLOTWIN_ID ) plasticityType
call plastic_dislotwin_dotState ( Mp , temperature ( ho ) % p ( tme ) , instance , of )
case ( PLASTICITY_DISLOTUNGSTEN_ID ) plasticityType
call plastic_disloTungsten_dotState ( Mp , temperature ( ho ) % p ( tme ) , instance , of )
case ( PLASTICITY_NONLOCAL_ID ) plasticityType
2020-12-23 01:15:27 +05:30
call plastic_nonlocal_dotState ( Mp , crystallite_partitionedF0 , temperature ( ho ) % p ( tme ) , subdt , &
2020-12-20 21:02:33 +05:30
instance , of , ip , el )
end select plasticityType
2020-12-23 11:44:07 +05:30
broken = any ( IEEE_is_NaN ( plasticState ( ph ) % dotState ( : , of ) ) )
2020-12-20 21:02:33 +05:30
2020-12-22 16:08:29 +05:30
end function mech_collectDotState
2020-12-20 21:02:33 +05:30
2020-12-20 21:41:43 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief for constitutive models having an instantaneous change of state
!> will return false if delta state is not needed/supported by the constitutive model
!--------------------------------------------------------------------------------------------------
2020-12-23 11:44:07 +05:30
function constitutive_deltaState ( S , Fi , co , ip , el , ph , of ) result ( broken )
2020-12-20 21:41:43 +05:30
integer , intent ( in ) :: &
2020-12-23 11:37:18 +05:30
co , & !< component-ID of integration point
2020-12-20 21:41:43 +05:30
ip , & !< integration point
el , & !< element
2020-12-23 11:44:07 +05:30
ph , &
2020-12-20 21:41:43 +05:30
of
real ( pReal ) , intent ( in ) , dimension ( 3 , 3 ) :: &
S , & !< 2nd Piola Kirchhoff stress
Fi !< intermediate deformation gradient
real ( pReal ) , dimension ( 3 , 3 ) :: &
Mp
integer :: &
instance , &
myOffset , &
mySize
logical :: &
broken
Mp = matmul ( matmul ( transpose ( Fi ) , Fi ) , S )
2020-12-23 11:44:07 +05:30
instance = phase_plasticityInstance ( ph )
2020-12-20 21:41:43 +05:30
2020-12-23 11:44:07 +05:30
plasticityType : select case ( phase_plasticity ( ph ) )
2020-12-20 21:41:43 +05:30
case ( PLASTICITY_KINEHARDENING_ID ) plasticityType
call plastic_kinehardening_deltaState ( Mp , instance , of )
2020-12-23 11:44:07 +05:30
broken = any ( IEEE_is_NaN ( plasticState ( ph ) % deltaState ( : , of ) ) )
2020-12-20 21:41:43 +05:30
case ( PLASTICITY_NONLOCAL_ID ) plasticityType
call plastic_nonlocal_deltaState ( Mp , instance , of , ip , el )
2020-12-23 11:44:07 +05:30
broken = any ( IEEE_is_NaN ( plasticState ( ph ) % deltaState ( : , of ) ) )
2020-12-20 21:41:43 +05:30
case default
broken = . false .
end select plasticityType
if ( . not . broken ) then
2020-12-23 11:44:07 +05:30
select case ( phase_plasticity ( ph ) )
2020-12-20 21:41:43 +05:30
case ( PLASTICITY_NONLOCAL_ID , PLASTICITY_KINEHARDENING_ID )
2020-12-23 11:44:07 +05:30
myOffset = plasticState ( ph ) % offsetDeltaState
mySize = plasticState ( ph ) % sizeDeltaState
plasticState ( ph ) % state ( myOffset + 1 : myOffset + mySize , of ) = &
plasticState ( ph ) % state ( myOffset + 1 : myOffset + mySize , of ) + plasticState ( ph ) % deltaState ( 1 : mySize , of )
2020-12-20 21:41:43 +05:30
end select
endif
end function constitutive_deltaState
2020-12-21 16:44:09 +05:30
module subroutine mech_results ( group , ph )
character ( len = * ) , intent ( in ) :: group
integer , intent ( in ) :: ph
2020-12-22 13:15:01 +05:30
if ( phase_plasticity ( ph ) / = PLASTICITY_NONE_ID ) &
2020-12-22 13:24:32 +05:30
call results_closeGroup ( results_addGroup ( group / / 'plastic/' ) )
2020-12-22 13:15:01 +05:30
2020-12-21 16:44:09 +05:30
select case ( phase_plasticity ( ph ) )
case ( PLASTICITY_ISOTROPIC_ID )
2020-12-22 13:24:32 +05:30
call plastic_isotropic_results ( phase_plasticityInstance ( ph ) , group / / 'plastic/' )
2020-12-21 16:44:09 +05:30
case ( PLASTICITY_PHENOPOWERLAW_ID )
2020-12-22 13:24:32 +05:30
call plastic_phenopowerlaw_results ( phase_plasticityInstance ( ph ) , group / / 'plastic/' )
2020-12-21 16:44:09 +05:30
case ( PLASTICITY_KINEHARDENING_ID )
2020-12-22 13:24:32 +05:30
call plastic_kinehardening_results ( phase_plasticityInstance ( ph ) , group / / 'plastic/' )
2020-12-21 16:44:09 +05:30
case ( PLASTICITY_DISLOTWIN_ID )
2020-12-22 13:24:32 +05:30
call plastic_dislotwin_results ( phase_plasticityInstance ( ph ) , group / / 'plastic/' )
2020-12-21 16:44:09 +05:30
case ( PLASTICITY_DISLOTUNGSTEN_ID )
2020-12-22 13:24:32 +05:30
call plastic_dislotungsten_results ( phase_plasticityInstance ( ph ) , group / / 'plastic/' )
2020-12-21 16:44:09 +05:30
case ( PLASTICITY_NONLOCAL_ID )
2020-12-22 13:24:32 +05:30
call plastic_nonlocal_results ( phase_plasticityInstance ( ph ) , group / / 'plastic/' )
2020-12-22 13:15:01 +05:30
2020-12-21 16:44:09 +05:30
end select
2020-12-22 13:24:32 +05:30
call crystallite_results ( group , ph )
2020-12-21 16:44:09 +05:30
end subroutine mech_results
module subroutine mech_restart_read ( fileHandle )
integer ( HID_T ) , intent ( in ) :: fileHandle
end subroutine mech_restart_read
2020-12-21 15:27:18 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief calculation of stress (P) with time integration based on a residuum in Lp and
!> intermediate acceleration of the Newton-Raphson correction
!--------------------------------------------------------------------------------------------------
2020-12-23 11:37:18 +05:30
function integrateStress ( co , ip , el , timeFraction ) result ( broken )
2020-12-21 15:27:18 +05:30
integer , intent ( in ) :: el , & ! element index
ip , & ! integration point index
2020-12-23 11:37:18 +05:30
co ! grain index
2020-12-21 15:27:18 +05:30
real ( pReal ) , optional , intent ( in ) :: timeFraction ! fraction of timestep
real ( pReal ) , dimension ( 3 , 3 ) :: F , & ! deformation gradient at end of timestep
Fp_new , & ! plastic deformation gradient at end of timestep
invFp_new , & ! inverse of Fp_new
invFp_current , & ! inverse of Fp_current
Lpguess , & ! current guess for plastic velocity gradient
Lpguess_old , & ! known last good guess for plastic velocity gradient
Lp_constitutive , & ! plastic velocity gradient resulting from constitutive law
residuumLp , & ! current residuum of plastic velocity gradient
residuumLp_old , & ! last residuum of plastic velocity gradient
deltaLp , & ! direction of next guess
Fi_new , & ! gradient of intermediate deformation stages
invFi_new , &
invFi_current , & ! inverse of Fi_current
Liguess , & ! current guess for intermediate velocity gradient
Liguess_old , & ! known last good guess for intermediate velocity gradient
Li_constitutive , & ! intermediate velocity gradient resulting from constitutive law
residuumLi , & ! current residuum of intermediate velocity gradient
residuumLi_old , & ! last residuum of intermediate velocity gradient
deltaLi , & ! direction of next guess
Fe , & ! elastic deformation gradient
S , & ! 2nd Piola-Kirchhoff Stress in plastic (lattice) configuration
A , &
B , &
temp_33
real ( pReal ) , dimension ( 9 ) :: temp_9 ! needed for matrix inversion by LAPACK
integer , dimension ( 9 ) :: devNull_9 ! needed for matrix inversion by LAPACK
real ( pReal ) , dimension ( 9 , 9 ) :: dRLp_dLp , & ! partial derivative of residuum (Jacobian for Newton-Raphson scheme)
dRLi_dLi ! partial derivative of residuumI (Jacobian for Newton-Raphson scheme)
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) :: dS_dFe , & ! partial derivative of 2nd Piola-Kirchhoff stress
dS_dFi , &
dFe_dLp , & ! partial derivative of elastic deformation gradient
dFe_dLi , &
dFi_dLi , &
dLp_dFi , &
dLi_dFi , &
dLp_dS , &
dLi_dS
real ( pReal ) steplengthLp , &
steplengthLi , &
dt , & ! time increment
atol_Lp , &
atol_Li , &
devNull
integer NiterationStressLp , & ! number of stress integrations
NiterationStressLi , & ! number of inner stress integrations
ierr , & ! error indicator for LAPACK
o , &
p , &
m , &
jacoCounterLp , &
jacoCounterLi ! counters to check for Jacobian update
logical :: error , broken
broken = . true .
if ( present ( timeFraction ) ) then
2020-12-23 11:37:18 +05:30
dt = crystallite_subdt ( co , ip , el ) * timeFraction
F = crystallite_subF0 ( 1 : 3 , 1 : 3 , co , ip , el ) &
+ ( crystallite_subF ( 1 : 3 , 1 : 3 , co , ip , el ) - crystallite_subF0 ( 1 : 3 , 1 : 3 , co , ip , el ) ) * timeFraction
2020-12-21 15:27:18 +05:30
else
2020-12-23 11:37:18 +05:30
dt = crystallite_subdt ( co , ip , el )
F = crystallite_subF ( 1 : 3 , 1 : 3 , co , ip , el )
2020-12-21 15:27:18 +05:30
endif
2020-12-23 11:37:18 +05:30
call constitutive_plastic_dependentState ( crystallite_partitionedF ( 1 : 3 , 1 : 3 , co , ip , el ) , co , ip , el )
2020-12-21 15:27:18 +05:30
2020-12-23 11:37:18 +05:30
p = material_phaseAt ( co , el )
m = material_phaseMemberAt ( co , ip , el )
2020-12-21 15:27:18 +05:30
2020-12-23 11:37:18 +05:30
Lpguess = crystallite_Lp ( 1 : 3 , 1 : 3 , co , ip , el ) ! take as first guess
2020-12-23 14:29:47 +05:30
Liguess = constitutive_mech_Li ( p ) % data ( 1 : 3 , 1 : 3 , m ) ! take as first guess
2020-12-21 15:27:18 +05:30
2020-12-23 11:37:18 +05:30
call math_invert33 ( invFp_current , devNull , error , crystallite_subFp0 ( 1 : 3 , 1 : 3 , co , ip , el ) )
2020-12-21 15:27:18 +05:30
if ( error ) return ! error
2020-12-23 11:37:18 +05:30
call math_invert33 ( invFi_current , devNull , error , crystallite_subFi0 ( 1 : 3 , 1 : 3 , co , ip , el ) )
2020-12-21 15:27:18 +05:30
if ( error ) return ! error
A = matmul ( F , invFp_current ) ! intermediate tensor needed later to calculate dFe_dLp
jacoCounterLi = 0
steplengthLi = 1.0_pReal
residuumLi_old = 0.0_pReal
Liguess_old = Liguess
NiterationStressLi = 0
LiLoop : do
NiterationStressLi = NiterationStressLi + 1
if ( NiterationStressLi > num % nStress ) return ! error
invFi_new = matmul ( invFi_current , math_I3 - dt * Liguess )
Fi_new = math_inv33 ( invFi_new )
jacoCounterLp = 0
steplengthLp = 1.0_pReal
residuumLp_old = 0.0_pReal
Lpguess_old = Lpguess
NiterationStressLp = 0
LpLoop : do
NiterationStressLp = NiterationStressLp + 1
if ( NiterationStressLp > num % nStress ) return ! error
B = math_I3 - dt * Lpguess
Fe = matmul ( matmul ( A , B ) , invFi_new )
call constitutive_hooke_SandItsTangents ( S , dS_dFe , dS_dFi , &
2020-12-23 11:37:18 +05:30
Fe , Fi_new , co , ip , el )
2020-12-21 15:27:18 +05:30
call constitutive_plastic_LpAndItsTangents ( Lp_constitutive , dLp_dS , dLp_dFi , &
2020-12-23 11:37:18 +05:30
S , Fi_new , co , ip , el )
2020-12-21 15:27:18 +05:30
!* update current residuum and check for convergence of loop
atol_Lp = max ( num % rtol_crystalliteStress * max ( norm2 ( Lpguess ) , norm2 ( Lp_constitutive ) ) , & ! absolute tolerance from largest acceptable relative error
num % atol_crystalliteStress ) ! minimum lower cutoff
residuumLp = Lpguess - Lp_constitutive
if ( any ( IEEE_is_NaN ( residuumLp ) ) ) then
return ! error
elseif ( norm2 ( residuumLp ) < atol_Lp ) then ! converged if below absolute tolerance
exit LpLoop
elseif ( NiterationStressLp == 1 . or . norm2 ( residuumLp ) < norm2 ( residuumLp_old ) ) then ! not converged, but improved norm of residuum (always proceed in first iteration)...
residuumLp_old = residuumLp ! ...remember old values and...
Lpguess_old = Lpguess
steplengthLp = 1.0_pReal ! ...proceed with normal step length (calculate new search direction)
else ! not converged and residuum not improved...
steplengthLp = num % subStepSizeLp * steplengthLp ! ...try with smaller step length in same direction
Lpguess = Lpguess_old &
+ deltaLp * stepLengthLp
cycle LpLoop
endif
calculateJacobiLi : if ( mod ( jacoCounterLp , num % iJacoLpresiduum ) == 0 ) then
jacoCounterLp = jacoCounterLp + 1
do o = 1 , 3 ; do p = 1 , 3
dFe_dLp ( o , 1 : 3 , p , 1 : 3 ) = - dt * A ( o , p ) * transpose ( invFi_new ) ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j)
enddo ; enddo
dRLp_dLp = math_eye ( 9 ) &
- math_3333to99 ( math_mul3333xx3333 ( math_mul3333xx3333 ( dLp_dS , dS_dFe ) , dFe_dLp ) )
temp_9 = math_33to9 ( residuumLp )
call dgesv ( 9 , 1 , dRLp_dLp , 9 , devNull_9 , temp_9 , 9 , ierr ) ! solve dRLp/dLp * delta Lp = -res for delta Lp
if ( ierr / = 0 ) return ! error
deltaLp = - math_9to33 ( temp_9 )
endif calculateJacobiLi
Lpguess = Lpguess &
+ deltaLp * steplengthLp
enddo LpLoop
call constitutive_LiAndItsTangents ( Li_constitutive , dLi_dS , dLi_dFi , &
2020-12-23 11:37:18 +05:30
S , Fi_new , co , ip , el )
2020-12-21 15:27:18 +05:30
!* update current residuum and check for convergence of loop
atol_Li = max ( num % rtol_crystalliteStress * max ( norm2 ( Liguess ) , norm2 ( Li_constitutive ) ) , & ! absolute tolerance from largest acceptable relative error
num % atol_crystalliteStress ) ! minimum lower cutoff
residuumLi = Liguess - Li_constitutive
if ( any ( IEEE_is_NaN ( residuumLi ) ) ) then
return ! error
elseif ( norm2 ( residuumLi ) < atol_Li ) then ! converged if below absolute tolerance
exit LiLoop
elseif ( NiterationStressLi == 1 . or . norm2 ( residuumLi ) < norm2 ( residuumLi_old ) ) then ! not converged, but improved norm of residuum (always proceed in first iteration)...
residuumLi_old = residuumLi ! ...remember old values and...
Liguess_old = Liguess
steplengthLi = 1.0_pReal ! ...proceed with normal step length (calculate new search direction)
else ! not converged and residuum not improved...
steplengthLi = num % subStepSizeLi * steplengthLi ! ...try with smaller step length in same direction
Liguess = Liguess_old &
+ deltaLi * steplengthLi
cycle LiLoop
endif
calculateJacobiLp : if ( mod ( jacoCounterLi , num % iJacoLpresiduum ) == 0 ) then
jacoCounterLi = jacoCounterLi + 1
temp_33 = matmul ( matmul ( A , B ) , invFi_current )
do o = 1 , 3 ; do p = 1 , 3
dFe_dLi ( 1 : 3 , o , 1 : 3 , p ) = - dt * math_I3 ( o , p ) * temp_33 ! dFe_dLp(i,j,k,l) = -dt * A(i,k) invFi(l,j)
dFi_dLi ( 1 : 3 , o , 1 : 3 , p ) = - dt * math_I3 ( o , p ) * invFi_current
enddo ; enddo
do o = 1 , 3 ; do p = 1 , 3
dFi_dLi ( 1 : 3 , 1 : 3 , o , p ) = matmul ( matmul ( Fi_new , dFi_dLi ( 1 : 3 , 1 : 3 , o , p ) ) , Fi_new )
enddo ; enddo
dRLi_dLi = math_eye ( 9 ) &
- math_3333to99 ( math_mul3333xx3333 ( dLi_dS , math_mul3333xx3333 ( dS_dFe , dFe_dLi ) &
+ math_mul3333xx3333 ( dS_dFi , dFi_dLi ) ) ) &
- math_3333to99 ( math_mul3333xx3333 ( dLi_dFi , dFi_dLi ) )
temp_9 = math_33to9 ( residuumLi )
call dgesv ( 9 , 1 , dRLi_dLi , 9 , devNull_9 , temp_9 , 9 , ierr ) ! solve dRLi/dLp * delta Li = -res for delta Li
if ( ierr / = 0 ) return ! error
deltaLi = - math_9to33 ( temp_9 )
endif calculateJacobiLp
Liguess = Liguess &
+ deltaLi * steplengthLi
enddo LiLoop
invFp_new = matmul ( invFp_current , B )
call math_invert33 ( Fp_new , devNull , error , invFp_new )
if ( error ) return ! error
2020-12-23 11:37:18 +05:30
p = material_phaseAt ( co , el )
m = material_phaseMemberAt ( co , ip , el )
2020-12-21 15:27:18 +05:30
2020-12-23 11:37:18 +05:30
crystallite_P ( 1 : 3 , 1 : 3 , co , ip , el ) = matmul ( matmul ( F , invFp_new ) , matmul ( S , transpose ( invFp_new ) ) )
crystallite_S ( 1 : 3 , 1 : 3 , co , ip , el ) = S
crystallite_Lp ( 1 : 3 , 1 : 3 , co , ip , el ) = Lpguess
2020-12-21 15:27:18 +05:30
constitutive_mech_Li ( p ) % data ( 1 : 3 , 1 : 3 , m ) = Liguess
2020-12-23 02:51:11 +05:30
constitutive_mech_Fp ( p ) % data ( 1 : 3 , 1 : 3 , m ) = Fp_new / math_det33 ( Fp_new ) ** ( 1.0_pReal / 3.0_pReal ) ! regularize
2020-12-21 15:27:18 +05:30
constitutive_mech_Fi ( p ) % data ( 1 : 3 , 1 : 3 , m ) = Fi_new
2020-12-23 11:37:18 +05:30
crystallite_Fe ( 1 : 3 , 1 : 3 , co , ip , el ) = matmul ( matmul ( F , invFp_new ) , invFi_new )
2020-12-21 15:27:18 +05:30
broken = . false .
end function integrateStress
!--------------------------------------------------------------------------------------------------
!> @brief integrate stress, state with adaptive 1st order explicit Euler method
!> using Fixed Point Iteration to adapt the stepsize
!--------------------------------------------------------------------------------------------------
2020-12-23 11:44:07 +05:30
module subroutine integrateStateFPI ( co , ip , el )
2020-12-21 15:27:18 +05:30
integer , intent ( in ) :: &
2020-12-23 11:44:07 +05:30
el , & !< element index in element loop
ip , & !< integration point index in ip loop
co !< grain index in grain loop
2020-12-21 15:27:18 +05:30
integer :: &
NiterationState , & !< number of iterations in state loop
2020-12-23 11:44:07 +05:30
ph , &
me , &
2020-12-21 15:27:18 +05:30
s , &
size_pl
integer , dimension ( maxval ( phase_Nsources ) ) :: &
size_so
real ( pReal ) :: &
zeta
real ( pReal ) , dimension ( max ( constitutive_plasticity_maxSizeDotState , constitutive_source_maxSizeDotState ) ) :: &
r ! state residuum
real ( pReal ) , dimension ( constitutive_plasticity_maxSizeDotState , 2 ) :: &
plastic_dotState
real ( pReal ) , dimension ( constitutive_source_maxSizeDotState , 2 , maxval ( phase_Nsources ) ) :: source_dotState
logical :: &
broken
2020-12-23 11:44:07 +05:30
ph = material_phaseAt ( co , el )
me = material_phaseMemberAt ( co , ip , el )
2020-12-21 15:27:18 +05:30
2020-12-23 11:44:07 +05:30
broken = mech_collectDotState ( crystallite_subdt ( co , ip , el ) , co , ip , el , ph , me )
2020-12-21 15:27:18 +05:30
if ( broken ) return
2020-12-23 11:44:07 +05:30
size_pl = plasticState ( ph ) % sizeDotState
plasticState ( ph ) % state ( 1 : size_pl , me ) = plasticState ( ph ) % subState0 ( 1 : size_pl , me ) &
+ plasticState ( ph ) % dotState ( 1 : size_pl , me ) &
* crystallite_subdt ( co , ip , el )
2020-12-21 15:27:18 +05:30
plastic_dotState ( 1 : size_pl , 2 ) = 0.0_pReal
iteration : do NiterationState = 1 , num % nState
if ( nIterationState > 1 ) plastic_dotState ( 1 : size_pl , 2 ) = plastic_dotState ( 1 : size_pl , 1 )
2020-12-23 11:44:07 +05:30
plastic_dotState ( 1 : size_pl , 1 ) = plasticState ( ph ) % dotState ( : , me )
2020-12-21 15:27:18 +05:30
2020-12-23 11:44:07 +05:30
broken = integrateStress ( co , ip , el )
2020-12-21 15:27:18 +05:30
if ( broken ) exit iteration
2020-12-23 11:44:07 +05:30
broken = mech_collectDotState ( crystallite_subdt ( co , ip , el ) , co , ip , el , ph , me )
2020-12-21 15:27:18 +05:30
if ( broken ) exit iteration
2020-12-23 11:44:07 +05:30
zeta = damper ( plasticState ( ph ) % dotState ( : , me ) , plastic_dotState ( 1 : size_pl , 1 ) , &
2020-12-21 15:27:18 +05:30
plastic_dotState ( 1 : size_pl , 2 ) )
2020-12-23 11:44:07 +05:30
plasticState ( ph ) % dotState ( : , me ) = plasticState ( ph ) % dotState ( : , me ) * zeta &
2020-12-21 15:27:18 +05:30
+ plastic_dotState ( 1 : size_pl , 1 ) * ( 1.0_pReal - zeta )
2020-12-23 11:44:07 +05:30
r ( 1 : size_pl ) = plasticState ( ph ) % state ( 1 : size_pl , me ) &
- plasticState ( ph ) % subState0 ( 1 : size_pl , me ) &
- plasticState ( ph ) % dotState ( 1 : size_pl , me ) * crystallite_subdt ( co , ip , el )
plasticState ( ph ) % state ( 1 : size_pl , me ) = plasticState ( ph ) % state ( 1 : size_pl , me ) &
2020-12-21 15:27:18 +05:30
- r ( 1 : size_pl )
2020-12-23 11:44:07 +05:30
crystallite_converged ( co , ip , el ) = converged ( r ( 1 : size_pl ) , &
plasticState ( ph ) % state ( 1 : size_pl , me ) , &
plasticState ( ph ) % atol ( 1 : size_pl ) )
2020-12-21 15:27:18 +05:30
2020-12-23 11:44:07 +05:30
if ( crystallite_converged ( co , ip , el ) ) then
broken = constitutive_deltaState ( crystallite_S ( 1 : 3 , 1 : 3 , co , ip , el ) , &
constitutive_mech_Fi ( ph ) % data ( 1 : 3 , 1 : 3 , me ) , co , ip , el , ph , me )
2020-12-21 15:27:18 +05:30
exit iteration
endif
enddo iteration
contains
!--------------------------------------------------------------------------------------------------
!> @brief calculate the damping for correction of state and dot state
!--------------------------------------------------------------------------------------------------
real ( pReal ) pure function damper ( current , previous , previous2 )
real ( pReal ) , dimension ( : ) , intent ( in ) :: &
current , previous , previous2
real ( pReal ) :: dot_prod12 , dot_prod22
dot_prod12 = dot_product ( current - previous , previous - previous2 )
dot_prod22 = dot_product ( previous - previous2 , previous - previous2 )
if ( ( dot_product ( current , previous ) < 0.0_pReal . or . dot_prod12 < 0.0_pReal ) . and . dot_prod22 > 0.0_pReal ) then
damper = 0.75_pReal + 0.25_pReal * tanh ( 2.0_pReal + 4.0_pReal * dot_prod12 / dot_prod22 )
else
damper = 1.0_pReal
endif
end function damper
end subroutine integrateStateFPI
!--------------------------------------------------------------------------------------------------
!> @brief integrate state with 1st order explicit Euler method
!--------------------------------------------------------------------------------------------------
2020-12-23 11:44:07 +05:30
subroutine integrateStateEuler ( co , ip , el )
2020-12-21 15:27:18 +05:30
integer , intent ( in ) :: &
2020-12-23 11:44:07 +05:30
el , & !< element index in element loop
ip , & !< integration point index in ip loop
co !< grain index in grain loop
2020-12-21 15:27:18 +05:30
integer :: &
2020-12-23 11:44:07 +05:30
ph , &
me , &
2020-12-21 15:27:18 +05:30
sizeDotState
logical :: &
broken
2020-12-23 11:44:07 +05:30
ph = material_phaseAt ( co , el )
me = material_phaseMemberAt ( co , ip , el )
2020-12-21 15:27:18 +05:30
2020-12-23 11:44:07 +05:30
broken = mech_collectDotState ( crystallite_subdt ( co , ip , el ) , co , ip , el , ph , me )
2020-12-21 15:27:18 +05:30
if ( broken ) return
2020-12-23 11:44:07 +05:30
sizeDotState = plasticState ( ph ) % sizeDotState
plasticState ( ph ) % state ( 1 : sizeDotState , me ) = plasticState ( ph ) % subState0 ( 1 : sizeDotState , me ) &
+ plasticState ( ph ) % dotState ( 1 : sizeDotState , me ) &
* crystallite_subdt ( co , ip , el )
2020-12-21 15:27:18 +05:30
2020-12-23 11:44:07 +05:30
broken = constitutive_deltaState ( crystallite_S ( 1 : 3 , 1 : 3 , co , ip , el ) , &
constitutive_mech_Fi ( ph ) % data ( 1 : 3 , 1 : 3 , me ) , co , ip , el , ph , me )
2020-12-21 15:27:18 +05:30
if ( broken ) return
2020-12-23 11:44:07 +05:30
broken = integrateStress ( co , ip , el )
crystallite_converged ( co , ip , el ) = . not . broken
2020-12-21 15:27:18 +05:30
end subroutine integrateStateEuler
!--------------------------------------------------------------------------------------------------
!> @brief integrate stress, state with 1st order Euler method with adaptive step size
!--------------------------------------------------------------------------------------------------
2020-12-23 11:44:07 +05:30
subroutine integrateStateAdaptiveEuler ( co , ip , el )
2020-12-21 15:27:18 +05:30
integer , intent ( in ) :: &
2020-12-23 11:44:07 +05:30
el , & !< element index in element loop
ip , & !< integration point index in ip loop
co !< grain index in grain loop
2020-12-21 15:27:18 +05:30
integer :: &
2020-12-23 11:44:07 +05:30
ph , &
me , &
2020-12-21 15:27:18 +05:30
sizeDotState
logical :: &
broken
real ( pReal ) , dimension ( constitutive_plasticity_maxSizeDotState ) :: residuum_plastic
2020-12-23 11:44:07 +05:30
ph = material_phaseAt ( co , el )
me = material_phaseMemberAt ( co , ip , el )
2020-12-21 15:27:18 +05:30
2020-12-23 11:44:07 +05:30
broken = mech_collectDotState ( crystallite_subdt ( co , ip , el ) , co , ip , el , ph , me )
2020-12-21 15:27:18 +05:30
if ( broken ) return
2020-12-23 11:44:07 +05:30
sizeDotState = plasticState ( ph ) % sizeDotState
2020-12-21 15:27:18 +05:30
2020-12-23 11:44:07 +05:30
residuum_plastic ( 1 : sizeDotState ) = - plasticState ( ph ) % dotstate ( 1 : sizeDotState , me ) * 0.5_pReal * crystallite_subdt ( co , ip , el )
plasticState ( ph ) % state ( 1 : sizeDotState , me ) = plasticState ( ph ) % subState0 ( 1 : sizeDotState , me ) &
+ plasticState ( ph ) % dotstate ( 1 : sizeDotState , me ) * crystallite_subdt ( co , ip , el )
2020-12-21 15:27:18 +05:30
2020-12-23 11:44:07 +05:30
broken = constitutive_deltaState ( crystallite_S ( 1 : 3 , 1 : 3 , co , ip , el ) , &
constitutive_mech_Fi ( ph ) % data ( 1 : 3 , 1 : 3 , me ) , co , ip , el , ph , me )
2020-12-21 15:27:18 +05:30
if ( broken ) return
2020-12-23 11:44:07 +05:30
broken = integrateStress ( co , ip , el )
2020-12-21 15:27:18 +05:30
if ( broken ) return
2020-12-23 11:44:07 +05:30
broken = mech_collectDotState ( crystallite_subdt ( co , ip , el ) , co , ip , el , ph , me )
2020-12-21 15:27:18 +05:30
if ( broken ) return
2020-12-23 11:44:07 +05:30
sizeDotState = plasticState ( ph ) % sizeDotState
crystallite_converged ( co , ip , el ) = converged ( residuum_plastic ( 1 : sizeDotState ) &
+ 0.5_pReal * plasticState ( ph ) % dotState ( : , me ) * crystallite_subdt ( co , ip , el ) , &
plasticState ( ph ) % state ( 1 : sizeDotState , me ) , &
plasticState ( ph ) % atol ( 1 : sizeDotState ) )
2020-12-21 15:27:18 +05:30
end subroutine integrateStateAdaptiveEuler
!---------------------------------------------------------------------------------------------------
!> @brief Integrate state (including stress integration) with the classic Runge Kutta method
!---------------------------------------------------------------------------------------------------
2020-12-23 11:44:07 +05:30
subroutine integrateStateRK4 ( co , ip , el )
2020-12-21 15:27:18 +05:30
2020-12-23 11:44:07 +05:30
integer , intent ( in ) :: co , ip , el
2020-12-21 15:27:18 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , parameter :: &
A = reshape ( [ &
0.5_pReal , 0.0_pReal , 0.0_pReal , &
0.0_pReal , 0.5_pReal , 0.0_pReal , &
0.0_pReal , 0.0_pReal , 1.0_pReal ] , &
shape ( A ) )
real ( pReal ) , dimension ( 3 ) , parameter :: &
C = [ 0.5_pReal , 0.5_pReal , 1.0_pReal ]
real ( pReal ) , dimension ( 4 ) , parameter :: &
B = [ 1.0_pReal / 6.0_pReal , 1.0_pReal / 3.0_pReal , 1.0_pReal / 3.0_pReal , 1.0_pReal / 6.0_pReal ]
2020-12-23 11:44:07 +05:30
call integrateStateRK ( co , ip , el , A , B , C )
2020-12-21 15:27:18 +05:30
end subroutine integrateStateRK4
!---------------------------------------------------------------------------------------------------
!> @brief Integrate state (including stress integration) with the Cash-Carp method
!---------------------------------------------------------------------------------------------------
2020-12-23 11:44:07 +05:30
subroutine integrateStateRKCK45 ( co , ip , el )
2020-12-21 15:27:18 +05:30
2020-12-23 11:44:07 +05:30
integer , intent ( in ) :: co , ip , el
2020-12-21 15:27:18 +05:30
real ( pReal ) , dimension ( 5 , 5 ) , parameter :: &
A = reshape ( [ &
1._pReal / 5._pReal , . 0_pReal , . 0_pReal , . 0_pReal , . 0_pReal , &
3._pReal / 4 0._pReal , 9._pReal / 4 0._pReal , . 0_pReal , . 0_pReal , . 0_pReal , &
3_pReal / 1 0._pReal , - 9._pReal / 1 0._pReal , 6._pReal / 5._pReal , . 0_pReal , . 0_pReal , &
- 1 1._pReal / 5 4._pReal , 5._pReal / 2._pReal , - 7 0.0_pReal / 2 7.0_pReal , 3 5.0_pReal / 2 7.0_pReal , . 0_pReal , &
163 1._pReal / 5529 6._pReal , 17 5._pReal / 51 2._pReal , 57 5._pReal / 1382 4._pReal , 4427 5._pReal / 11059 2._pReal , 25 3._pReal / 409 6._pReal ] , &
shape ( A ) )
real ( pReal ) , dimension ( 5 ) , parameter :: &
C = [ 0.2_pReal , 0.3_pReal , 0.6_pReal , 1.0_pReal , 0.875_pReal ]
real ( pReal ) , dimension ( 6 ) , parameter :: &
B = &
[ 3 7.0_pReal / 37 8.0_pReal , . 0_pReal , 25 0.0_pReal / 62 1.0_pReal , &
12 5.0_pReal / 59 4.0_pReal , . 0_pReal , 51 2.0_pReal / 177 1.0_pReal ] , &
DB = B - &
[ 282 5.0_pReal / 2764 8.0_pReal , . 0_pReal , 1857 5.0_pReal / 4838 4.0_pReal , &
1352 5.0_pReal / 5529 6.0_pReal , 27 7.0_pReal / 1433 6.0_pReal , 1._pReal / 4._pReal ]
2020-12-23 11:44:07 +05:30
call integrateStateRK ( co , ip , el , A , B , C , DB )
2020-12-21 15:27:18 +05:30
end subroutine integrateStateRKCK45
!--------------------------------------------------------------------------------------------------
!> @brief Integrate state (including stress integration) with an explicit Runge-Kutta method or an
!! embedded explicit Runge-Kutta method
!--------------------------------------------------------------------------------------------------
2020-12-23 11:44:07 +05:30
subroutine integrateStateRK ( co , ip , el , A , B , CC , DB )
2020-12-21 15:27:18 +05:30
real ( pReal ) , dimension ( : , : ) , intent ( in ) :: A
real ( pReal ) , dimension ( : ) , intent ( in ) :: B , CC
real ( pReal ) , dimension ( : ) , intent ( in ) , optional :: DB
integer , intent ( in ) :: &
2020-12-23 11:44:07 +05:30
el , & !< element index in element loop
ip , & !< integration point index in ip loop
co !< grain index in grain loop
2020-12-21 15:27:18 +05:30
integer :: &
stage , & ! stage index in integration stage loop
n , &
2020-12-23 11:44:07 +05:30
ph , &
me , &
2020-12-21 15:27:18 +05:30
sizeDotState
logical :: &
broken
real ( pReal ) , dimension ( constitutive_plasticity_maxSizeDotState , size ( B ) ) :: plastic_RKdotState
2020-12-23 11:44:07 +05:30
ph = material_phaseAt ( co , el )
me = material_phaseMemberAt ( co , ip , el )
2020-12-21 15:27:18 +05:30
2020-12-23 11:44:07 +05:30
broken = mech_collectDotState ( crystallite_subdt ( co , ip , el ) , co , ip , el , ph , me )
2020-12-21 15:27:18 +05:30
if ( broken ) return
do stage = 1 , size ( A , 1 )
2020-12-23 11:44:07 +05:30
sizeDotState = plasticState ( ph ) % sizeDotState
plastic_RKdotState ( 1 : sizeDotState , stage ) = plasticState ( ph ) % dotState ( : , me )
plasticState ( ph ) % dotState ( : , me ) = A ( 1 , stage ) * plastic_RKdotState ( 1 : sizeDotState , 1 )
2020-12-21 15:27:18 +05:30
do n = 2 , stage
2020-12-23 11:44:07 +05:30
sizeDotState = plasticState ( ph ) % sizeDotState
plasticState ( ph ) % dotState ( : , me ) = plasticState ( ph ) % dotState ( : , me ) &
2020-12-21 15:27:18 +05:30
+ A ( n , stage ) * plastic_RKdotState ( 1 : sizeDotState , n )
enddo
2020-12-23 11:44:07 +05:30
sizeDotState = plasticState ( ph ) % sizeDotState
plasticState ( ph ) % state ( 1 : sizeDotState , me ) = plasticState ( ph ) % subState0 ( 1 : sizeDotState , me ) &
+ plasticState ( ph ) % dotState ( 1 : sizeDotState , me ) &
* crystallite_subdt ( co , ip , el )
2020-12-21 15:27:18 +05:30
2020-12-23 11:44:07 +05:30
broken = integrateStress ( co , ip , el , CC ( stage ) )
2020-12-21 15:27:18 +05:30
if ( broken ) exit
2020-12-23 11:44:07 +05:30
broken = mech_collectDotState ( crystallite_subdt ( co , ip , el ) * CC ( stage ) , co , ip , el , ph , me )
2020-12-21 15:27:18 +05:30
if ( broken ) exit
enddo
if ( broken ) return
2020-12-23 11:44:07 +05:30
sizeDotState = plasticState ( ph ) % sizeDotState
2020-12-21 15:27:18 +05:30
2020-12-23 11:44:07 +05:30
plastic_RKdotState ( 1 : sizeDotState , size ( B ) ) = plasticState ( ph ) % dotState ( : , me )
plasticState ( ph ) % dotState ( : , me ) = matmul ( plastic_RKdotState ( 1 : sizeDotState , 1 : size ( B ) ) , B )
plasticState ( ph ) % state ( 1 : sizeDotState , me ) = plasticState ( ph ) % subState0 ( 1 : sizeDotState , me ) &
+ plasticState ( ph ) % dotState ( 1 : sizeDotState , me ) &
* crystallite_subdt ( co , ip , el )
2020-12-21 15:27:18 +05:30
if ( present ( DB ) ) &
broken = . not . converged ( matmul ( plastic_RKdotState ( 1 : sizeDotState , 1 : size ( DB ) ) , DB ) &
2020-12-23 11:44:07 +05:30
* crystallite_subdt ( co , ip , el ) , &
plasticState ( ph ) % state ( 1 : sizeDotState , me ) , &
plasticState ( ph ) % atol ( 1 : sizeDotState ) )
2020-12-21 15:27:18 +05:30
if ( broken ) return
2020-12-23 11:44:07 +05:30
broken = constitutive_deltaState ( crystallite_S ( 1 : 3 , 1 : 3 , co , ip , el ) , &
constitutive_mech_Fi ( ph ) % data ( 1 : 3 , 1 : 3 , me ) , co , ip , el , ph , me )
2020-12-21 15:27:18 +05:30
if ( broken ) return
2020-12-23 11:44:07 +05:30
broken = integrateStress ( co , ip , el )
crystallite_converged ( co , ip , el ) = . not . broken
2020-12-21 15:27:18 +05:30
end subroutine integrateStateRK
2020-12-22 13:24:32 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief writes crystallite results to HDF5 output file
!--------------------------------------------------------------------------------------------------
subroutine crystallite_results ( group , ph )
character ( len = * ) , intent ( in ) :: group
integer , intent ( in ) :: ph
integer :: ou
real ( pReal ) , allocatable , dimension ( : , : , : ) :: selected_tensors
real ( pReal ) , allocatable , dimension ( : , : ) :: selected_rotations
character ( len = : ) , allocatable :: structureLabel
call results_closeGroup ( results_addGroup ( group / / '/mechanics/' ) )
do ou = 1 , size ( output_constituent ( ph ) % label )
select case ( output_constituent ( ph ) % label ( ou ) )
case ( 'F' )
selected_tensors = select_tensors ( crystallite_partitionedF , ph )
call results_writeDataset ( group / / '/mechanics/' , selected_tensors , output_constituent ( ph ) % label ( ou ) , &
'deformation gradient' , '1' )
case ( 'F_e' )
selected_tensors = select_tensors ( crystallite_Fe , ph )
call results_writeDataset ( group / / '/mechanics/' , selected_tensors , output_constituent ( ph ) % label ( ou ) , &
'elastic deformation gradient' , '1' )
case ( 'F_p' )
2020-12-23 02:51:11 +05:30
call results_writeDataset ( group / / '/mechanics/' , constitutive_mech_Fp ( ph ) % data , output_constituent ( ph ) % label ( ou ) , &
2020-12-22 13:24:32 +05:30
'plastic deformation gradient' , '1' )
case ( 'F_i' )
call results_writeDataset ( group / / '/mechanics/' , constitutive_mech_Fi ( ph ) % data , output_constituent ( ph ) % label ( ou ) , &
'inelastic deformation gradient' , '1' )
case ( 'L_p' )
selected_tensors = select_tensors ( crystallite_Lp , ph )
call results_writeDataset ( group / / '/mechanics/' , selected_tensors , output_constituent ( ph ) % label ( ou ) , &
'plastic velocity gradient' , '1/s' )
case ( 'L_i' )
call results_writeDataset ( group / / '/mechanics/' , constitutive_mech_Li ( ph ) % data , output_constituent ( ph ) % label ( ou ) , &
'inelastic velocity gradient' , '1/s' )
case ( 'P' )
selected_tensors = select_tensors ( crystallite_P , ph )
call results_writeDataset ( group / / '/mechanics/' , selected_tensors , output_constituent ( ph ) % label ( ou ) , &
'First Piola-Kirchhoff stress' , 'Pa' )
case ( 'S' )
selected_tensors = select_tensors ( crystallite_S , ph )
call results_writeDataset ( group / / '/mechanics/' , selected_tensors , output_constituent ( ph ) % label ( ou ) , &
'Second Piola-Kirchhoff stress' , 'Pa' )
case ( 'O' )
select case ( lattice_structure ( ph ) )
case ( lattice_ISO_ID )
structureLabel = 'aP'
case ( lattice_FCC_ID )
structureLabel = 'cF'
case ( lattice_BCC_ID )
structureLabel = 'cI'
case ( lattice_BCT_ID )
structureLabel = 'tI'
case ( lattice_HEX_ID )
structureLabel = 'hP'
case ( lattice_ORT_ID )
structureLabel = 'oP'
end select
selected_rotations = select_rotations ( crystallite_orientation , ph )
call results_writeDataset ( group / / '/mechanics/' , selected_rotations , output_constituent ( ph ) % label ( ou ) , &
'crystal orientation as quaternion' , 'q_0 (q_1 q_2 q_3)' )
call results_addAttribute ( 'Lattice' , structureLabel , group / / '/mechanics/' / / output_constituent ( ph ) % label ( ou ) )
end select
enddo
contains
!------------------------------------------------------------------------------------------------
!> @brief select tensors for output
!------------------------------------------------------------------------------------------------
function select_tensors ( dataset , instance )
integer , intent ( in ) :: instance
real ( pReal ) , dimension ( : , : , : , : , : ) , intent ( in ) :: dataset
real ( pReal ) , allocatable , dimension ( : , : , : ) :: select_tensors
integer :: e , i , c , j
allocate ( select_tensors ( 3 , 3 , count ( material_phaseAt == instance ) * discretization_nIPs ) )
j = 0
do e = 1 , size ( material_phaseAt , 2 )
do i = 1 , discretization_nIPs
do c = 1 , size ( material_phaseAt , 1 ) !ToDo: this needs to be changed for varying Ngrains
if ( material_phaseAt ( c , e ) == instance ) then
j = j + 1
select_tensors ( 1 : 3 , 1 : 3 , j ) = dataset ( 1 : 3 , 1 : 3 , c , i , e )
endif
enddo
enddo
enddo
end function select_tensors
!--------------------------------------------------------------------------------------------------
!> @brief select rotations for output
!--------------------------------------------------------------------------------------------------
function select_rotations ( dataset , instance )
integer , intent ( in ) :: instance
type ( rotation ) , dimension ( : , : , : ) , intent ( in ) :: dataset
real ( pReal ) , allocatable , dimension ( : , : ) :: select_rotations
integer :: e , i , c , j
allocate ( select_rotations ( 4 , count ( material_phaseAt == instance ) * homogenization_maxNconstituents * discretization_nIPs ) )
j = 0
do e = 1 , size ( material_phaseAt , 2 )
do i = 1 , discretization_nIPs
do c = 1 , size ( material_phaseAt , 1 ) !ToDo: this needs to be changed for varying Ngrains
if ( material_phaseAt ( c , e ) == instance ) then
j = j + 1
select_rotations ( 1 : 4 , j ) = dataset ( c , i , e ) % asQuaternion ( )
endif
enddo
enddo
enddo
end function select_rotations
end subroutine crystallite_results
2020-12-22 14:33:19 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Backup data for homog cutback.
!--------------------------------------------------------------------------------------------------
module subroutine mech_initializeRestorationPoints ( ph , me )
2020-12-23 11:44:07 +05:30
integer , intent ( in ) :: ph , me
2020-12-22 14:33:19 +05:30
2020-12-22 23:43:30 +05:30
2020-12-22 14:33:19 +05:30
constitutive_mech_partionedFi0 ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = constitutive_mech_Fi0 ( ph ) % data ( 1 : 3 , 1 : 3 , me )
2020-12-23 02:51:11 +05:30
constitutive_mech_partionedFp0 ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = constitutive_mech_Fp0 ( ph ) % data ( 1 : 3 , 1 : 3 , me )
2020-12-22 14:33:19 +05:30
constitutive_mech_partionedLi0 ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = constitutive_mech_Li0 ( ph ) % data ( 1 : 3 , 1 : 3 , me )
plasticState ( ph ) % partitionedState0 ( : , me ) = plasticState ( ph ) % state0 ( : , me )
end subroutine mech_initializeRestorationPoints
2020-12-23 03:57:56 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Wind homog inc forward.
!--------------------------------------------------------------------------------------------------
module subroutine constitutive_mech_windForward ( ph , me )
integer , intent ( in ) :: ph , me
2020-12-23 11:44:07 +05:30
constitutive_mech_partionedFp0 ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = constitutive_mech_Fp ( ph ) % data ( 1 : 3 , 1 : 3 , me )
constitutive_mech_partionedFi0 ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = constitutive_mech_Fi ( ph ) % data ( 1 : 3 , 1 : 3 , me )
constitutive_mech_partionedLi0 ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = constitutive_mech_Li ( ph ) % data ( 1 : 3 , 1 : 3 , me )
2020-12-23 03:57:56 +05:30
2020-12-23 11:44:07 +05:30
plasticState ( ph ) % partitionedState0 ( : , me ) = plasticState ( ph ) % state ( : , me )
2020-12-23 03:57:56 +05:30
end subroutine constitutive_mech_windForward
2020-12-23 11:28:54 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Forward data after successful increment.
! ToDo: Any guessing for the current states possible?
!--------------------------------------------------------------------------------------------------
module subroutine constitutive_mech_forward ( )
integer :: ph
do ph = 1 , size ( plasticState )
plasticState ( ph ) % state0 = plasticState ( ph ) % state
constitutive_mech_Fi0 ( ph ) = constitutive_mech_Fi ( ph )
constitutive_mech_Fp0 ( ph ) = constitutive_mech_Fp ( ph )
constitutive_mech_Li0 ( ph ) = constitutive_mech_Li ( ph )
enddo
end subroutine constitutive_mech_forward
2020-11-03 05:19:17 +05:30
end submodule constitutive_mech
2020-07-09 04:31:08 +05:30