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
!----------------------------------------------------------------------------------------------------
2021-02-13 23:27:41 +05:30
submodule ( phase ) mechanical
2020-07-09 04:31:08 +05:30
2021-02-13 23:11:30 +05:30
2020-12-23 19:07:12 +05:30
enum , bind ( c ) ; enumerator :: &
ELASTICITY_UNDEFINED_ID , &
ELASTICITY_HOOKE_ID , &
STIFFNESS_DEGRADATION_UNDEFINED_ID , &
2021-01-08 04:20:06 +05:30
STIFFNESS_DEGRADATION_DAMAGE_ID , &
PLASTICITY_UNDEFINED_ID , &
PLASTICITY_NONE_ID , &
PLASTICITY_ISOTROPIC_ID , &
PLASTICITY_PHENOPOWERLAW_ID , &
PLASTICITY_KINEHARDENING_ID , &
PLASTICITY_DISLOTWIN_ID , &
PLASTICITY_DISLOTUNGSTEN_ID , &
2021-01-27 13:15:54 +05:30
PLASTICITY_NONLOCAL_ID , &
KINEMATICS_UNDEFINED_ID , &
KINEMATICS_CLEAVAGE_OPENING_ID , &
KINEMATICS_SLIPPLANE_OPENING_ID , &
KINEMATICS_THERMAL_EXPANSION_ID
2020-12-23 19:07:12 +05:30
end enum
2020-12-23 22:00:19 +05:30
2021-01-08 04:20:06 +05:30
integer ( kind ( ELASTICITY_UNDEFINED_ID ) ) , dimension ( : ) , allocatable :: &
2020-12-15 22:15:11 +05:30
phase_elasticity !< elasticity of each phase
2021-01-08 04:20:06 +05:30
integer ( kind ( STIFFNESS_DEGRADATION_UNDEFINED_ID ) ) , dimension ( : , : ) , allocatable :: &
2020-12-15 22:15:11 +05:30
phase_stiffnessDegradation !< active stiffness degradation mechanisms of each phase
2020-12-30 14:24:06 +05:30
type ( tTensorContainer ) , dimension ( : ) , allocatable :: &
! current value
2021-02-09 03:51:53 +05:30
phase_mechanical_Fe , &
phase_mechanical_Fi , &
phase_mechanical_Fp , &
phase_mechanical_F , &
phase_mechanical_Li , &
phase_mechanical_Lp , &
phase_mechanical_S , &
phase_mechanical_P , &
2020-12-30 14:24:06 +05:30
! converged value at end of last solver increment
2021-02-09 03:51:53 +05:30
phase_mechanical_Fi0 , &
phase_mechanical_Fp0 , &
phase_mechanical_F0 , &
phase_mechanical_Li0 , &
phase_mechanical_Lp0 , &
phase_mechanical_S0
2021-01-08 04:20:06 +05:30
integer ( kind ( PLASTICITY_undefined_ID ) ) , dimension ( : ) , allocatable :: &
phase_plasticity !< plasticity of each phase
2020-07-09 04:31:08 +05:30
interface
2021-01-27 05:02:44 +05:30
module subroutine eigendeformation_init ( phases )
class ( tNode ) , pointer :: phases
end subroutine eigendeformation_init
module subroutine plastic_init
end subroutine plastic_init
2021-02-14 11:36:14 +05:30
module subroutine plastic_isotropic_LiAndItsTangent ( Li , dLi_dMi , Mi , ph , me )
2021-01-27 11:40:53 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( out ) :: &
Li !< inleastic velocity gradient
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) , intent ( out ) :: &
dLi_dMi !< derivative of Li with respect to Mandel stress
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: &
Mi !< Mandel stress
integer , intent ( in ) :: &
2021-02-14 11:36:14 +05:30
ph , &
2021-01-27 11:40:53 +05:30
me
end subroutine plastic_isotropic_LiAndItsTangent
2020-07-09 04:31:08 +05:30
2021-01-26 16:47:00 +05:30
module function plastic_dotState ( subdt , co , ip , el , ph , me ) result ( broken )
2020-12-20 21:02:33 +05:30
integer , intent ( in ) :: &
2021-01-26 16:15:39 +05:30
co , & !< component-ID of integration point
ip , & !< integration point
el , & !< element
ph , &
2021-01-26 13:09:17 +05:30
me
2021-01-26 16:15:39 +05:30
real ( pReal ) , intent ( in ) :: &
subdt !< timestep
logical :: broken
2021-01-26 16:47:00 +05:30
end function plastic_dotState
2020-07-09 04:31:08 +05:30
2021-02-14 21:59:23 +05:30
module function plastic_deltaState ( ph , me ) result ( broken )
2020-07-09 04:31:08 +05:30
integer , intent ( in ) :: &
2021-01-26 16:15:39 +05:30
ph , &
2021-01-27 13:15:54 +05:30
me
2021-01-26 16:15:39 +05:30
logical :: &
broken
2021-01-26 16:47:00 +05:30
end function plastic_deltaState
2020-12-20 21:41:43 +05:30
2021-02-09 03:51:53 +05:30
module subroutine phase_LiAndItsTangents ( Li , dLi_dS , dLi_dFi , &
2021-02-13 11:45:57 +05:30
S , Fi , ph , me )
2021-01-27 11:40:53 +05:30
integer , intent ( in ) :: &
2021-02-13 11:45:57 +05:30
ph , me
2021-01-27 11:40:53 +05:30
real ( pReal ) , intent ( in ) , dimension ( 3 , 3 ) :: &
S !< 2nd Piola-Kirchhoff stress
real ( pReal ) , intent ( in ) , dimension ( 3 , 3 ) :: &
Fi !< intermediate deformation gradient
real ( pReal ) , intent ( out ) , dimension ( 3 , 3 ) :: &
Li !< intermediate velocity gradient
real ( pReal ) , intent ( out ) , dimension ( 3 , 3 , 3 , 3 ) :: &
dLi_dS , & !< derivative of Li with respect to S
dLi_dFi
2021-02-09 03:51:53 +05:30
end subroutine phase_LiAndItsTangents
2021-01-27 11:40:53 +05:30
2021-01-26 13:09:17 +05:30
2021-01-26 16:11:19 +05:30
module subroutine plastic_LpAndItsTangents ( Lp , dLp_dS , dLp_dFi , &
2021-02-14 21:59:23 +05:30
S , Fi , ph , me )
2021-01-26 13:09:17 +05:30
integer , intent ( in ) :: &
2021-02-14 21:59:23 +05:30
ph , me
2021-01-26 13:09:17 +05:30
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
2021-01-26 16:11:19 +05:30
end subroutine plastic_LpAndItsTangents
2021-01-26 13:09:17 +05:30
2021-01-27 04:36:41 +05:30
2021-02-14 11:36:14 +05:30
module subroutine plastic_isotropic_results ( ph , group )
integer , intent ( in ) :: ph
2020-07-12 18:52:40 +05:30
character ( len = * ) , intent ( in ) :: group
end subroutine plastic_isotropic_results
2021-02-14 05:20:42 +05:30
module subroutine plastic_phenopowerlaw_results ( ph , group )
integer , intent ( in ) :: ph
2020-07-12 18:52:40 +05:30
character ( len = * ) , intent ( in ) :: group
end subroutine plastic_phenopowerlaw_results
2021-02-14 05:20:42 +05:30
module subroutine plastic_kinehardening_results ( ph , group )
integer , intent ( in ) :: ph
2020-07-12 18:52:40 +05:30
character ( len = * ) , intent ( in ) :: group
end subroutine plastic_kinehardening_results
2021-02-14 11:36:14 +05:30
module subroutine plastic_dislotwin_results ( ph , group )
integer , intent ( in ) :: ph
2020-07-12 18:52:40 +05:30
character ( len = * ) , intent ( in ) :: group
end subroutine plastic_dislotwin_results
2021-02-14 11:36:14 +05:30
module subroutine plastic_dislotungsten_results ( ph , group )
integer , intent ( in ) :: ph
2020-07-12 18:52:40 +05:30
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
2021-02-14 19:59:10 +05:30
module subroutine plastic_nonlocal_results ( ph , group )
integer , intent ( in ) :: ph
2020-07-12 18:52:40 +05:30
character ( len = * ) , intent ( in ) :: group
end subroutine plastic_nonlocal_results
2021-01-19 14:44:34 +05:30
module function plastic_dislotwin_homogenizedC ( ph , me ) result ( homogenizedC )
real ( pReal ) , dimension ( 6 , 6 ) :: homogenizedC
integer , intent ( in ) :: ph , me
2020-12-23 22:00:19 +05:30
end function plastic_dislotwin_homogenizedC
2020-07-12 18:52:40 +05:30
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
2020-12-24 14:52:41 +05:30
procedure ( integrateStateFPI ) , pointer :: integrateState
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
!--------------------------------------------------------------------------------------------------
2021-02-09 03:51:53 +05:30
module subroutine mechanical_init ( phases )
2020-12-30 17:04:00 +05:30
class ( tNode ) , pointer :: &
phases
2020-07-09 04:31:08 +05:30
2020-11-03 02:09:33 +05:30
integer :: &
2020-12-29 23:32:22 +05:30
el , &
ip , &
co , &
2020-12-23 11:44:07 +05:30
ph , &
2020-12-29 23:32:22 +05:30
me , &
stiffDegradationCtr , &
Nconstituents
2020-11-03 02:09:33 +05:30
class ( tNode ) , pointer :: &
2020-12-22 23:54:00 +05:30
num_crystallite , &
2020-11-03 02:09:33 +05:30
phase , &
2020-11-03 03:16:46 +05:30
mech , &
2020-11-03 02:09:33 +05:30
elastic , &
2021-01-27 11:40:53 +05:30
stiffDegradation
2020-08-15 19:32:10 +05:30
2021-02-13 23:27:41 +05:30
print '(/,a)' , ' <<<+- phase:mechanical 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?
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
2021-02-09 03:51:53 +05:30
allocate ( phase_mechanical_Fe ( phases % length ) )
allocate ( phase_mechanical_Fi ( phases % length ) )
allocate ( phase_mechanical_Fi0 ( phases % length ) )
allocate ( phase_mechanical_Fp ( phases % length ) )
allocate ( phase_mechanical_Fp0 ( phases % length ) )
allocate ( phase_mechanical_F ( phases % length ) )
allocate ( phase_mechanical_F0 ( phases % length ) )
allocate ( phase_mechanical_Li ( phases % length ) )
allocate ( phase_mechanical_Li0 ( phases % length ) )
allocate ( phase_mechanical_Lp0 ( phases % length ) )
allocate ( phase_mechanical_Lp ( phases % length ) )
allocate ( phase_mechanical_S ( phases % length ) )
allocate ( phase_mechanical_P ( phases % length ) )
allocate ( phase_mechanical_S0 ( phases % length ) )
2020-12-29 23:32:22 +05:30
2020-12-23 11:44:07 +05:30
do ph = 1 , phases % length
2020-12-29 23:32:22 +05:30
Nconstituents = count ( material_phaseAt == ph ) * discretization_nIPs
2021-02-09 03:51:53 +05:30
allocate ( phase_mechanical_Fi ( ph ) % data ( 3 , 3 , Nconstituents ) )
allocate ( phase_mechanical_Fe ( ph ) % data ( 3 , 3 , Nconstituents ) )
allocate ( phase_mechanical_Fi0 ( ph ) % data ( 3 , 3 , Nconstituents ) )
allocate ( phase_mechanical_Fp ( ph ) % data ( 3 , 3 , Nconstituents ) )
allocate ( phase_mechanical_Fp0 ( ph ) % data ( 3 , 3 , Nconstituents ) )
allocate ( phase_mechanical_Li ( ph ) % data ( 3 , 3 , Nconstituents ) )
allocate ( phase_mechanical_Li0 ( ph ) % data ( 3 , 3 , Nconstituents ) )
allocate ( phase_mechanical_Lp0 ( ph ) % data ( 3 , 3 , Nconstituents ) )
allocate ( phase_mechanical_Lp ( ph ) % data ( 3 , 3 , Nconstituents ) )
allocate ( phase_mechanical_S ( ph ) % data ( 3 , 3 , Nconstituents ) , source = 0.0_pReal )
allocate ( phase_mechanical_P ( ph ) % data ( 3 , 3 , Nconstituents ) , source = 0.0_pReal )
allocate ( phase_mechanical_S0 ( ph ) % data ( 3 , 3 , Nconstituents ) , source = 0.0_pReal )
allocate ( phase_mechanical_F ( ph ) % data ( 3 , 3 , Nconstituents ) )
allocate ( phase_mechanical_F0 ( ph ) % data ( 3 , 3 , Nconstituents ) )
2020-12-30 04:44:48 +05:30
2020-12-23 11:44:07 +05:30
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-12-29 23:32:22 +05:30
!$OMP PARALLEL DO PRIVATE(ph,me)
do el = 1 , size ( material_phaseMemberAt , 3 ) ; do ip = 1 , size ( material_phaseMemberAt , 2 )
do co = 1 , homogenization_Nconstituents ( material_homogenizationAt ( el ) )
ph = material_phaseAt ( co , el )
me = material_phaseMemberAt ( co , ip , el )
2020-12-30 04:44:48 +05:30
2021-02-09 03:51:53 +05:30
phase_mechanical_Fp0 ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = material_orientation0 ( co , ip , el ) % asMatrix ( ) ! Fp reflects initial orientation (see 10.1016/j.actamat.2006.01.005)
phase_mechanical_Fp0 ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = phase_mechanical_Fp0 ( ph ) % data ( 1 : 3 , 1 : 3 , me ) &
/ math_det33 ( phase_mechanical_Fp0 ( ph ) % data ( 1 : 3 , 1 : 3 , me ) ) ** ( 1.0_pReal / 3.0_pReal )
phase_mechanical_Fi0 ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = math_I3
phase_mechanical_F0 ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = math_I3
2020-12-30 04:44:48 +05:30
2021-02-09 03:51:53 +05:30
phase_mechanical_Fe ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = math_inv33 ( matmul ( phase_mechanical_Fi0 ( ph ) % data ( 1 : 3 , 1 : 3 , me ) , &
phase_mechanical_Fp0 ( ph ) % data ( 1 : 3 , 1 : 3 , me ) ) ) ! assuming that euler angles are given in internal strain free configuration
phase_mechanical_Fp ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = phase_mechanical_Fp0 ( ph ) % data ( 1 : 3 , 1 : 3 , me )
phase_mechanical_Fi ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = phase_mechanical_Fi0 ( ph ) % data ( 1 : 3 , 1 : 3 , me )
phase_mechanical_F ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = phase_mechanical_F0 ( ph ) % data ( 1 : 3 , 1 : 3 , me )
2020-12-29 23:32:22 +05:30
enddo
enddo ; enddo
!$OMP END PARALLEL DO
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_localPlasticity ( phases % length ) , source = . true . )
2021-01-27 05:02:44 +05:30
call plastic_init ( )
2020-08-15 19:32:10 +05:30
2020-12-23 11:44:07 +05:30
do ph = 1 , phases % length
phase_elasticityInstance ( ph ) = count ( phase_elasticity ( 1 : ph ) == phase_elasticity ( 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
2021-01-26 01:01:12 +05:30
2021-01-27 05:02:44 +05:30
call eigendeformation_init ( phases )
2021-01-26 01:01:12 +05:30
2021-02-09 03:51:53 +05:30
end subroutine mechanical_init
2020-07-09 04:31:08 +05:30
2020-07-10 18:29:07 +05:30
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
!--------------------------------------------------------------------------------------------------
2021-02-09 03:51:53 +05:30
subroutine phase_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
2021-01-19 14:44:34 +05:30
2020-11-03 02:09:33 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: E
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) :: C
integer :: &
ho , & !< homogenization
2021-01-19 14:44:34 +05:30
d , & !< counter in degradation loop
i , j , ph , me
2020-11-03 02:09:33 +05:30
ho = material_homogenizationAt ( el )
2021-02-09 03:51:53 +05:30
C = math_66toSym3333 ( phase_homogenizedC ( material_phaseAt ( co , el ) , material_phaseMemberAt ( 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
2021-02-12 14:01:02 +05:30
C = C * phase_damage_get_phi ( co , 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
2021-02-09 03:51:53 +05:30
end subroutine phase_hooke_SandItsTangents
2020-11-03 02:09:33 +05:30
2021-02-09 03:51:53 +05:30
module subroutine mechanical_results ( group , ph )
2020-12-21 16:44:09 +05:30
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 )
2021-02-14 11:36:14 +05:30
call plastic_isotropic_results ( ph , group / / 'plastic/' )
2020-12-21 16:44:09 +05:30
case ( PLASTICITY_PHENOPOWERLAW_ID )
2021-02-14 05:20:42 +05:30
call plastic_phenopowerlaw_results ( ph , group / / 'plastic/' )
2020-12-21 16:44:09 +05:30
case ( PLASTICITY_KINEHARDENING_ID )
2021-02-14 05:20:42 +05:30
call plastic_kinehardening_results ( ph , group / / 'plastic/' )
2020-12-21 16:44:09 +05:30
case ( PLASTICITY_DISLOTWIN_ID )
2021-02-14 11:36:14 +05:30
call plastic_dislotwin_results ( ph , group / / 'plastic/' )
2020-12-21 16:44:09 +05:30
case ( PLASTICITY_DISLOTUNGSTEN_ID )
2021-02-14 11:36:14 +05:30
call plastic_dislotungsten_results ( ph , group / / 'plastic/' )
2020-12-21 16:44:09 +05:30
case ( PLASTICITY_NONLOCAL_ID )
2021-02-14 19:59:10 +05:30
call plastic_nonlocal_results ( 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 )
2021-02-09 03:51:53 +05:30
end subroutine mechanical_results
2020-12-21 16:44:09 +05:30
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-29 04:43:49 +05:30
function integrateStress ( F , subFp0 , subFi0 , Delta_t , co , ip , el ) result ( broken )
2020-12-21 15:27:18 +05:30
2020-12-29 04:43:49 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: F , subFp0 , subFi0
2020-12-24 15:06:48 +05:30
real ( pReal ) , intent ( in ) :: Delta_t
2020-12-21 15:27:18 +05:30
integer , intent ( in ) :: el , & ! element index
ip , & ! integration point index
2020-12-24 15:06:48 +05:30
co ! grain index
2020-12-21 15:27:18 +05:30
2020-12-24 15:06:48 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: Fp_new , & ! plastic deformation gradient at end of timestep
2020-12-21 15:27:18 +05:30
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 , &
atol_Lp , &
atol_Li , &
devNull
integer NiterationStressLp , & ! number of stress integrations
NiterationStressLi , & ! number of inner stress integrations
ierr , & ! error indicator for LAPACK
o , &
p , &
2020-12-24 13:23:02 +05:30
ph , &
me , &
2020-12-21 15:27:18 +05:30
jacoCounterLp , &
jacoCounterLi ! counters to check for Jacobian update
logical :: error , broken
2020-12-29 22:57:24 +05:30
broken = . true .
2020-12-21 15:27:18 +05:30
2020-12-24 13:23:02 +05:30
ph = material_phaseAt ( co , el )
me = material_phaseMemberAt ( co , ip , el )
2020-12-30 04:44:48 +05:30
2021-01-26 16:47:00 +05:30
call plastic_dependentState ( co , ip , el )
2020-12-21 15:27:18 +05:30
2021-02-09 03:51:53 +05:30
Lpguess = phase_mechanical_Lp ( ph ) % data ( 1 : 3 , 1 : 3 , me ) ! take as first guess
Liguess = phase_mechanical_Li ( ph ) % data ( 1 : 3 , 1 : 3 , me ) ! take as first guess
2020-12-21 15:27:18 +05:30
2020-12-29 04:43:49 +05:30
call math_invert33 ( invFp_current , devNull , error , subFp0 )
2020-12-21 15:27:18 +05:30
if ( error ) return ! error
2020-12-29 04:43:49 +05:30
call math_invert33 ( invFi_current , devNull , error , subFi0 )
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
2020-12-24 15:06:48 +05:30
invFi_new = matmul ( invFi_current , math_I3 - Delta_t * Liguess )
2020-12-21 15:27:18 +05:30
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
2020-12-24 15:06:48 +05:30
B = math_I3 - Delta_t * Lpguess
2020-12-21 15:27:18 +05:30
Fe = matmul ( matmul ( A , B ) , invFi_new )
2021-02-09 03:51:53 +05:30
call phase_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
2021-01-26 16:11:19 +05:30
call plastic_LpAndItsTangents ( Lp_constitutive , dLp_dS , dLp_dFi , &
2021-02-14 21:59:23 +05:30
S , Fi_new , ph , me )
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
2020-12-24 15:06:48 +05:30
dFe_dLp ( o , 1 : 3 , p , 1 : 3 ) = - Delta_t * A ( o , p ) * transpose ( invFi_new ) ! dFe_dLp(i,j,k,l) = -Delta_t * A(i,k) invFi(l,j)
2020-12-21 15:27:18 +05:30
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
2021-02-09 03:51:53 +05:30
call phase_LiAndItsTangents ( Li_constitutive , dLi_dS , dLi_dFi , &
2021-02-13 11:45:57 +05:30
S , Fi_new , ph , me )
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
2020-12-24 15:06:48 +05:30
dFe_dLi ( 1 : 3 , o , 1 : 3 , p ) = - Delta_t * math_I3 ( o , p ) * temp_33 ! dFe_dLp(i,j,k,l) = -Delta_t * A(i,k) invFi(l,j)
dFi_dLi ( 1 : 3 , o , 1 : 3 , p ) = - Delta_t * math_I3 ( o , p ) * invFi_current
2020-12-21 15:27:18 +05:30
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
2021-02-09 03:51:53 +05:30
phase_mechanical_P ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = matmul ( matmul ( F , invFp_new ) , matmul ( S , transpose ( invFp_new ) ) )
phase_mechanical_S ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = S
phase_mechanical_Lp ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = Lpguess
phase_mechanical_Li ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = Liguess
phase_mechanical_Fp ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = Fp_new / math_det33 ( Fp_new ) ** ( 1.0_pReal / 3.0_pReal ) ! regularize
phase_mechanical_Fi ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = Fi_new
phase_mechanical_Fe ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = 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-29 04:43:49 +05:30
function integrateStateFPI ( F_0 , F , subFp0 , subFi0 , subState0 , Delta_t , co , ip , el ) result ( broken )
2020-12-21 15:27:18 +05:30
2020-12-29 04:43:49 +05:30
real ( pReal ) , intent ( in ) , dimension ( 3 , 3 ) :: F_0 , F , subFp0 , subFi0
real ( pReal ) , intent ( in ) , dimension ( : ) :: subState0
2020-12-24 15:50:34 +05:30
real ( pReal ) , intent ( in ) :: Delta_t
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-29 02:11:48 +05:30
logical :: &
broken
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-29 03:03:04 +05:30
sizeDotState
2020-12-21 15:27:18 +05:30
real ( pReal ) :: &
zeta
2021-02-09 03:51:53 +05:30
real ( pReal ) , dimension ( phase_plasticity_maxSizeDotState ) :: &
2020-12-21 15:27:18 +05:30
r ! state residuum
2021-02-09 03:51:53 +05:30
real ( pReal ) , dimension ( phase_plasticity_maxSizeDotState , 2 ) :: &
2020-12-29 03:03:04 +05:30
dotState
2020-12-29 02:11:48 +05:30
2020-12-21 15:27:18 +05:30
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
2021-01-26 16:47:00 +05:30
broken = plastic_dotState ( Delta_t , co , ip , el , ph , me )
2020-12-21 15:27:18 +05:30
if ( broken ) return
2020-12-29 03:03:04 +05:30
sizeDotState = plasticState ( ph ) % sizeDotState
2020-12-29 04:43:49 +05:30
plasticState ( ph ) % state ( 1 : sizeDotState , me ) = subState0 &
2020-12-29 03:03:04 +05:30
+ plasticState ( ph ) % dotState ( 1 : sizeDotState , me ) * Delta_t
dotState ( 1 : sizeDotState , 2 ) = 0.0_pReal
2020-12-21 15:27:18 +05:30
iteration : do NiterationState = 1 , num % nState
2020-12-29 03:03:04 +05:30
if ( nIterationState > 1 ) dotState ( 1 : sizeDotState , 2 ) = dotState ( 1 : sizeDotState , 1 )
dotState ( 1 : sizeDotState , 1 ) = plasticState ( ph ) % dotState ( : , me )
2020-12-21 15:27:18 +05:30
2020-12-29 04:43:49 +05:30
broken = integrateStress ( F , subFp0 , subFi0 , Delta_t , co , ip , el )
2020-12-21 15:27:18 +05:30
if ( broken ) exit iteration
2021-01-26 16:47:00 +05:30
broken = plastic_dotState ( Delta_t , co , ip , el , ph , me )
2020-12-21 15:27:18 +05:30
if ( broken ) exit iteration
2020-12-29 03:03:04 +05:30
zeta = damper ( plasticState ( ph ) % dotState ( : , me ) , dotState ( 1 : sizeDotState , 1 ) , &
dotState ( 1 : sizeDotState , 2 ) )
2020-12-23 11:44:07 +05:30
plasticState ( ph ) % dotState ( : , me ) = plasticState ( ph ) % dotState ( : , me ) * zeta &
2020-12-29 03:03:04 +05:30
+ dotState ( 1 : sizeDotState , 1 ) * ( 1.0_pReal - zeta )
r ( 1 : sizeDotState ) = plasticState ( ph ) % state ( 1 : sizeDotState , me ) &
2020-12-29 04:43:49 +05:30
- subState0 &
2020-12-29 03:03:04 +05:30
- plasticState ( ph ) % dotState ( 1 : sizeDotState , me ) * Delta_t
plasticState ( ph ) % state ( 1 : sizeDotState , me ) = plasticState ( ph ) % state ( 1 : sizeDotState , me ) &
- r ( 1 : sizeDotState )
if ( converged ( r ( 1 : sizeDotState ) , plasticState ( ph ) % state ( 1 : sizeDotState , me ) , plasticState ( ph ) % atol ( 1 : sizeDotState ) ) ) then
2021-02-14 21:59:23 +05:30
broken = plastic_deltaState ( ph , me )
2020-12-21 15:27:18 +05:30
exit iteration
endif
enddo iteration
2021-01-06 23:11:25 +05:30
2020-12-21 15:27:18 +05:30
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
2020-12-27 20:43:31 +05:30
end function integrateStateFPI
2020-12-21 15:27:18 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief integrate state with 1st order explicit Euler method
!--------------------------------------------------------------------------------------------------
2020-12-29 04:43:49 +05:30
function integrateStateEuler ( F_0 , F , subFp0 , subFi0 , subState0 , Delta_t , co , ip , el ) result ( broken )
2020-12-21 15:27:18 +05:30
2020-12-29 04:43:49 +05:30
real ( pReal ) , intent ( in ) , dimension ( 3 , 3 ) :: F_0 , F , subFp0 , subFi0
real ( pReal ) , intent ( in ) , dimension ( : ) :: subState0
2020-12-24 15:50:34 +05:30
real ( pReal ) , intent ( in ) :: Delta_t
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-29 02:11:48 +05:30
logical :: &
broken
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
2020-12-29 02:11:48 +05:30
2020-12-21 15:27:18 +05:30
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
2021-01-26 16:47:00 +05:30
broken = plastic_dotState ( Delta_t , 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-29 04:43:49 +05:30
plasticState ( ph ) % state ( 1 : sizeDotState , me ) = subState0 &
+ plasticState ( ph ) % dotState ( 1 : sizeDotState , me ) * Delta_t
2020-12-21 15:27:18 +05:30
2021-02-14 21:59:23 +05:30
broken = plastic_deltaState ( ph , me )
2020-12-21 15:27:18 +05:30
if ( broken ) return
2020-12-29 04:43:49 +05:30
broken = integrateStress ( F , subFp0 , subFi0 , Delta_t , co , ip , el )
2020-12-21 15:27:18 +05:30
2020-12-27 20:43:31 +05:30
end function integrateStateEuler
2020-12-21 15:27:18 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief integrate stress, state with 1st order Euler method with adaptive step size
!--------------------------------------------------------------------------------------------------
2020-12-29 04:43:49 +05:30
function integrateStateAdaptiveEuler ( F_0 , F , subFp0 , subFi0 , subState0 , Delta_t , co , ip , el ) result ( broken )
2020-12-21 15:27:18 +05:30
2020-12-29 04:43:49 +05:30
real ( pReal ) , intent ( in ) , dimension ( 3 , 3 ) :: F_0 , F , subFp0 , subFi0
real ( pReal ) , intent ( in ) , dimension ( : ) :: subState0
2020-12-24 15:50:34 +05:30
real ( pReal ) , intent ( in ) :: Delta_t
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-29 02:11:48 +05:30
logical :: &
broken
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
2021-02-09 03:51:53 +05:30
real ( pReal ) , dimension ( phase_plasticity_maxSizeDotState ) :: residuum_plastic
2020-12-21 15:27:18 +05:30
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
2021-01-26 16:47:00 +05:30
broken = plastic_dotState ( Delta_t , 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-24 15:50:34 +05:30
residuum_plastic ( 1 : sizeDotState ) = - plasticState ( ph ) % dotstate ( 1 : sizeDotState , me ) * 0.5_pReal * Delta_t
2020-12-29 04:43:49 +05:30
plasticState ( ph ) % state ( 1 : sizeDotState , me ) = subState0 &
2020-12-29 02:11:48 +05:30
+ plasticState ( ph ) % dotstate ( 1 : sizeDotState , me ) * Delta_t
2020-12-21 15:27:18 +05:30
2021-02-14 21:59:23 +05:30
broken = plastic_deltaState ( ph , me )
2020-12-21 15:27:18 +05:30
if ( broken ) return
2020-12-29 04:43:49 +05:30
broken = integrateStress ( F , subFp0 , subFi0 , Delta_t , co , ip , el )
2020-12-21 15:27:18 +05:30
if ( broken ) return
2021-01-26 16:47:00 +05:30
broken = plastic_dotState ( Delta_t , co , ip , el , ph , me )
2020-12-21 15:27:18 +05:30
if ( broken ) return
2020-12-27 20:43:31 +05:30
broken = . not . converged ( residuum_plastic ( 1 : sizeDotState ) + 0.5_pReal * plasticState ( ph ) % dotState ( : , me ) * Delta_t , &
plasticState ( ph ) % state ( 1 : sizeDotState , me ) , &
plasticState ( ph ) % atol ( 1 : sizeDotState ) )
2020-12-21 15:27:18 +05:30
2020-12-27 20:43:31 +05:30
end function integrateStateAdaptiveEuler
2020-12-21 15:27:18 +05:30
!---------------------------------------------------------------------------------------------------
!> @brief Integrate state (including stress integration) with the classic Runge Kutta method
!---------------------------------------------------------------------------------------------------
2020-12-29 04:43:49 +05:30
function integrateStateRK4 ( F_0 , F , subFp0 , subFi0 , subState0 , Delta_t , co , ip , el ) result ( broken )
2020-12-21 15:27:18 +05:30
2020-12-29 04:43:49 +05:30
real ( pReal ) , intent ( in ) , dimension ( 3 , 3 ) :: F_0 , F , subFp0 , subFi0
real ( pReal ) , intent ( in ) , dimension ( : ) :: subState0
2020-12-24 15:50:34 +05:30
real ( pReal ) , intent ( in ) :: Delta_t
2020-12-23 11:44:07 +05:30
integer , intent ( in ) :: co , ip , el
2020-12-27 20:43:31 +05:30
logical :: broken
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-29 02:11:48 +05:30
2020-12-29 04:43:49 +05:30
broken = integrateStateRK ( F_0 , F , subFp0 , subFi0 , subState0 , Delta_t , co , ip , el , A , B , C )
2020-12-21 15:27:18 +05:30
2020-12-27 20:43:31 +05:30
end function integrateStateRK4
2020-12-21 15:27:18 +05:30
!---------------------------------------------------------------------------------------------------
!> @brief Integrate state (including stress integration) with the Cash-Carp method
!---------------------------------------------------------------------------------------------------
2020-12-29 04:43:49 +05:30
function integrateStateRKCK45 ( F_0 , F , subFp0 , subFi0 , subState0 , Delta_t , co , ip , el ) result ( broken )
2020-12-21 15:27:18 +05:30
2020-12-29 04:43:49 +05:30
real ( pReal ) , intent ( in ) , dimension ( 3 , 3 ) :: F_0 , F , subFp0 , subFi0
real ( pReal ) , intent ( in ) , dimension ( : ) :: subState0
2020-12-24 15:50:34 +05:30
real ( pReal ) , intent ( in ) :: Delta_t
2020-12-23 11:44:07 +05:30
integer , intent ( in ) :: co , ip , el
2020-12-27 20:43:31 +05:30
logical :: broken
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-29 02:11:48 +05:30
2020-12-29 04:43:49 +05:30
broken = integrateStateRK ( F_0 , F , subFp0 , subFi0 , subState0 , Delta_t , co , ip , el , A , B , C , DB )
2020-12-21 15:27:18 +05:30
2020-12-27 20:43:31 +05:30
end function integrateStateRKCK45
2020-12-21 15:27:18 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Integrate state (including stress integration) with an explicit Runge-Kutta method or an
!! embedded explicit Runge-Kutta method
!--------------------------------------------------------------------------------------------------
2020-12-29 04:43:49 +05:30
function integrateStateRK ( F_0 , F , subFp0 , subFi0 , subState0 , Delta_t , co , ip , el , A , B , C , DB ) result ( broken )
2020-12-21 15:27:18 +05:30
2020-12-29 04:43:49 +05:30
real ( pReal ) , intent ( in ) , dimension ( 3 , 3 ) :: F_0 , F , subFp0 , subFi0
real ( pReal ) , intent ( in ) , dimension ( : ) :: subState0
2020-12-24 15:50:34 +05:30
real ( pReal ) , intent ( in ) :: Delta_t
2020-12-21 15:27:18 +05:30
real ( pReal ) , dimension ( : , : ) , intent ( in ) :: A
2020-12-24 13:23:02 +05:30
real ( pReal ) , dimension ( : ) , intent ( in ) :: B , C
2020-12-21 15:27:18 +05:30
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-27 20:43:31 +05:30
logical :: broken
2020-12-23 12:42:56 +05:30
2020-12-27 20:43:31 +05:30
integer :: &
2020-12-21 15:27:18 +05:30
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
2021-02-09 03:51:53 +05:30
real ( pReal ) , dimension ( phase_plasticity_maxSizeDotState , size ( B ) ) :: plastic_RKdotState
2020-12-23 12:42:56 +05:30
2020-12-21 15:27:18 +05:30
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
2021-01-26 16:47:00 +05:30
broken = plastic_dotState ( Delta_t , co , ip , el , ph , me )
2020-12-21 15:27:18 +05:30
if ( broken ) return
2020-12-29 02:11:48 +05:30
sizeDotState = plasticState ( ph ) % sizeDotState
2020-12-24 15:50:34 +05:30
do stage = 1 , size ( A , 1 )
2020-12-29 02:11:48 +05:30
2020-12-23 11:44:07 +05:30
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
plasticState ( ph ) % dotState ( : , me ) = plasticState ( ph ) % dotState ( : , me ) &
2020-12-29 02:11:48 +05:30
+ A ( n , stage ) * plastic_RKdotState ( 1 : sizeDotState , n )
2020-12-21 15:27:18 +05:30
enddo
2020-12-29 04:43:49 +05:30
plasticState ( ph ) % state ( 1 : sizeDotState , me ) = subState0 &
2020-12-24 15:50:34 +05:30
+ plasticState ( ph ) % dotState ( 1 : sizeDotState , me ) * Delta_t
2020-12-21 15:27:18 +05:30
2020-12-29 04:43:49 +05:30
broken = integrateStress ( F_0 + ( F - F_0 ) * Delta_t * C ( stage ) , subFp0 , subFi0 , Delta_t * C ( stage ) , co , ip , el )
2020-12-21 15:27:18 +05:30
if ( broken ) exit
2021-01-26 16:47:00 +05:30
broken = plastic_dotState ( Delta_t * C ( 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
plastic_RKdotState ( 1 : sizeDotState , size ( B ) ) = plasticState ( ph ) % dotState ( : , me )
plasticState ( ph ) % dotState ( : , me ) = matmul ( plastic_RKdotState ( 1 : sizeDotState , 1 : size ( B ) ) , B )
2020-12-29 04:43:49 +05:30
plasticState ( ph ) % state ( 1 : sizeDotState , me ) = subState0 &
2020-12-24 15:50:34 +05:30
+ plasticState ( ph ) % dotState ( 1 : sizeDotState , me ) * Delta_t
2020-12-21 15:27:18 +05:30
if ( present ( DB ) ) &
2020-12-24 15:50:34 +05:30
broken = . not . converged ( matmul ( plastic_RKdotState ( 1 : sizeDotState , 1 : size ( DB ) ) , DB ) * Delta_t , &
plasticState ( ph ) % state ( 1 : sizeDotState , me ) , &
plasticState ( ph ) % atol ( 1 : sizeDotState ) )
2020-12-21 15:27:18 +05:30
if ( broken ) return
2021-02-14 21:59:23 +05:30
broken = plastic_deltaState ( ph , me )
2020-12-21 15:27:18 +05:30
if ( broken ) return
2020-12-29 04:43:49 +05:30
broken = integrateStress ( F , subFp0 , subFi0 , Delta_t , co , ip , el )
2020-12-21 15:27:18 +05:30
2020-12-27 20:43:31 +05:30
end function integrateStateRK
2020-12-21 15:27:18 +05:30
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_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' )
2021-02-09 03:51:53 +05:30
call results_writeDataset ( group / / '/mechanics/' , phase_mechanical_F ( ph ) % data , 'F' , &
2020-12-22 13:24:32 +05:30
'deformation gradient' , '1' )
case ( 'F_e' )
2021-02-09 03:51:53 +05:30
call results_writeDataset ( group / / '/mechanics/' , phase_mechanical_Fe ( ph ) % data , 'F_e' , &
2020-12-22 13:24:32 +05:30
'elastic deformation gradient' , '1' )
case ( 'F_p' )
2021-02-09 03:51:53 +05:30
call results_writeDataset ( group / / '/mechanics/' , phase_mechanical_Fp ( ph ) % data , 'F_p' , &
2020-12-22 13:24:32 +05:30
'plastic deformation gradient' , '1' )
case ( 'F_i' )
2021-02-09 03:51:53 +05:30
call results_writeDataset ( group / / '/mechanics/' , phase_mechanical_Fi ( ph ) % data , 'F_i' , &
2020-12-22 13:24:32 +05:30
'inelastic deformation gradient' , '1' )
case ( 'L_p' )
2021-02-09 03:51:53 +05:30
call results_writeDataset ( group / / '/mechanics/' , phase_mechanical_Lp ( ph ) % data , 'L_p' , &
2020-12-22 13:24:32 +05:30
'plastic velocity gradient' , '1/s' )
case ( 'L_i' )
2021-02-09 03:51:53 +05:30
call results_writeDataset ( group / / '/mechanics/' , phase_mechanical_Li ( ph ) % data , 'L_i' , &
2020-12-22 13:24:32 +05:30
'inelastic velocity gradient' , '1/s' )
case ( 'P' )
2021-02-09 03:51:53 +05:30
call results_writeDataset ( group / / '/mechanics/' , phase_mechanical_P ( ph ) % data , 'P' , &
2020-12-22 13:24:32 +05:30
'First Piola-Kirchhoff stress' , 'Pa' )
case ( 'S' )
2021-02-09 03:51:53 +05:30
call results_writeDataset ( group / / '/mechanics/' , phase_mechanical_S ( ph ) % data , 'S' , &
2020-12-22 13:24:32 +05:30
'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 rotations for output
!--------------------------------------------------------------------------------------------------
2020-12-24 13:23:02 +05:30
function select_rotations ( dataset , ph )
2020-12-22 13:24:32 +05:30
2020-12-24 13:23:02 +05:30
integer , intent ( in ) :: ph
2020-12-22 13:24:32 +05:30
type ( rotation ) , dimension ( : , : , : ) , intent ( in ) :: dataset
real ( pReal ) , allocatable , dimension ( : , : ) :: select_rotations
2020-12-24 13:23:02 +05:30
integer :: el , ip , co , j
2020-12-22 13:24:32 +05:30
2020-12-24 13:23:02 +05:30
allocate ( select_rotations ( 4 , count ( material_phaseAt == ph ) * homogenization_maxNconstituents * discretization_nIPs ) )
2020-12-22 13:24:32 +05:30
j = 0
2020-12-24 13:23:02 +05:30
do el = 1 , size ( material_phaseAt , 2 )
do ip = 1 , discretization_nIPs
do co = 1 , size ( material_phaseAt , 1 ) !ToDo: this needs to be changed for varying Ngrains
if ( material_phaseAt ( co , el ) == ph ) then
2020-12-22 13:24:32 +05:30
j = j + 1
2020-12-24 13:23:02 +05:30
select_rotations ( 1 : 4 , j ) = dataset ( co , ip , el ) % asQuaternion ( )
2020-12-22 13:24:32 +05:30
endif
enddo
enddo
enddo
end function select_rotations
end subroutine crystallite_results
2020-12-23 03:57:56 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Wind homog inc forward.
!--------------------------------------------------------------------------------------------------
2021-02-09 03:51:53 +05:30
module subroutine mechanical_windForward ( ph , me )
2020-12-23 03:57:56 +05:30
integer , intent ( in ) :: ph , me
2021-02-09 03:51:53 +05:30
phase_mechanical_Fp0 ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = phase_mechanical_Fp ( ph ) % data ( 1 : 3 , 1 : 3 , me )
phase_mechanical_Fi0 ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = phase_mechanical_Fi ( ph ) % data ( 1 : 3 , 1 : 3 , me )
phase_mechanical_F0 ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = phase_mechanical_F ( ph ) % data ( 1 : 3 , 1 : 3 , me )
phase_mechanical_Li0 ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = phase_mechanical_Li ( ph ) % data ( 1 : 3 , 1 : 3 , me )
phase_mechanical_Lp0 ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = phase_mechanical_Lp ( ph ) % data ( 1 : 3 , 1 : 3 , me )
phase_mechanical_S0 ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = phase_mechanical_S ( ph ) % data ( 1 : 3 , 1 : 3 , me )
2020-12-23 03:57:56 +05:30
2021-01-25 18:33:49 +05:30
plasticState ( ph ) % State0 ( : , me ) = plasticState ( ph ) % state ( : , me )
2020-12-23 03:57:56 +05:30
2021-02-09 03:51:53 +05:30
end subroutine mechanical_windForward
2020-12-23 03:57:56 +05:30
2020-12-23 11:28:54 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Forward data after successful increment.
! ToDo: Any guessing for the current states possible?
!--------------------------------------------------------------------------------------------------
2021-02-09 03:51:53 +05:30
module subroutine mechanical_forward ( )
2020-12-23 11:28:54 +05:30
integer :: ph
do ph = 1 , size ( plasticState )
2021-02-09 03:51:53 +05:30
phase_mechanical_Fi0 ( ph ) = phase_mechanical_Fi ( ph )
phase_mechanical_Fp0 ( ph ) = phase_mechanical_Fp ( ph )
phase_mechanical_F0 ( ph ) = phase_mechanical_F ( ph )
phase_mechanical_Li0 ( ph ) = phase_mechanical_Li ( ph )
phase_mechanical_Lp0 ( ph ) = phase_mechanical_Lp ( ph )
phase_mechanical_S0 ( ph ) = phase_mechanical_S ( ph )
2020-12-29 13:56:24 +05:30
plasticState ( ph ) % state0 = plasticState ( ph ) % state
2020-12-23 11:28:54 +05:30
enddo
2021-02-09 03:51:53 +05:30
end subroutine mechanical_forward
2020-12-23 11:28:54 +05:30
2020-12-23 22:00:19 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief returns the homogenize elasticity matrix
!> ToDo: homogenizedC66 would be more consistent
!--------------------------------------------------------------------------------------------------
2021-02-09 03:51:53 +05:30
module function phase_homogenizedC ( ph , me ) result ( C )
2020-12-23 22:00:19 +05:30
real ( pReal ) , dimension ( 6 , 6 ) :: C
2021-01-19 14:44:34 +05:30
integer , intent ( in ) :: ph , me
2020-12-23 22:00:19 +05:30
2021-02-09 03:51:53 +05:30
plasticType : select case ( phase_plasticity ( ph ) )
case ( PLASTICITY_DISLOTWIN_ID ) plasticType
2021-01-19 14:44:34 +05:30
C = plastic_dislotwin_homogenizedC ( ph , me )
2021-02-09 03:51:53 +05:30
case default plasticType
2021-01-19 14:44:34 +05:30
C = lattice_C66 ( 1 : 6 , 1 : 6 , ph )
2021-02-09 03:51:53 +05:30
end select plasticType
2020-12-23 22:00:19 +05:30
2021-02-09 03:51:53 +05:30
end function phase_homogenizedC
2020-12-23 22:00:19 +05:30
2020-12-24 14:52:41 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief calculate stress (P)
!--------------------------------------------------------------------------------------------------
2020-12-27 20:43:31 +05:30
module function crystallite_stress ( dt , co , ip , el ) result ( converged_ )
2020-12-24 14:52:41 +05:30
real ( pReal ) , intent ( in ) :: dt
integer , intent ( in ) :: &
co , &
ip , &
el
2020-12-27 20:43:31 +05:30
logical :: converged_
2020-12-24 14:52:41 +05:30
real ( pReal ) :: &
formerSubStep
integer :: &
2021-02-13 13:02:43 +05:30
ph , me , sizeDotState
2020-12-24 14:52:41 +05:30
logical :: todo
real ( pReal ) :: subFrac , subStep
real ( pReal ) , dimension ( 3 , 3 ) :: &
2020-12-29 04:43:49 +05:30
subFp0 , &
subFi0 , &
subLp0 , &
subLi0 , &
2020-12-28 03:26:21 +05:30
subF0 , &
subF
2020-12-29 04:43:49 +05:30
real ( pReal ) , dimension ( : ) , allocatable :: subState0
2020-12-24 14:52:41 +05:30
ph = material_phaseAt ( co , el )
me = material_phaseMemberAt ( co , ip , el )
2020-12-29 04:43:49 +05:30
sizeDotState = plasticState ( ph ) % sizeDotState
2021-02-09 03:51:53 +05:30
subLi0 = phase_mechanical_Li0 ( ph ) % data ( 1 : 3 , 1 : 3 , me )
subLp0 = phase_mechanical_Lp0 ( ph ) % data ( 1 : 3 , 1 : 3 , me )
2021-01-25 18:33:49 +05:30
subState0 = plasticState ( ph ) % State0 ( : , me )
2020-12-29 04:43:49 +05:30
2021-02-13 23:11:30 +05:30
if ( damageState ( ph ) % sizeState > 0 ) &
2021-02-13 13:02:43 +05:30
damageState ( ph ) % subState0 ( : , me ) = damageState ( ph ) % state0 ( : , me )
2020-12-24 14:52:41 +05:30
2021-02-09 03:51:53 +05:30
subFp0 = phase_mechanical_Fp0 ( ph ) % data ( 1 : 3 , 1 : 3 , me )
subFi0 = phase_mechanical_Fi0 ( ph ) % data ( 1 : 3 , 1 : 3 , me )
subF0 = phase_mechanical_F0 ( ph ) % data ( 1 : 3 , 1 : 3 , me )
2020-12-24 14:52:41 +05:30
subFrac = 0.0_pReal
subStep = 1.0_pReal / num % subStepSizeCryst
todo = . true .
2020-12-27 20:43:31 +05:30
converged_ = . false . ! pretend failed step of 1/subStepSizeCryst
2020-12-24 14:52:41 +05:30
todo = . true .
cutbackLooping : do while ( todo )
2020-12-27 20:43:31 +05:30
if ( converged_ ) then
2020-12-24 14:52:41 +05:30
formerSubStep = subStep
subFrac = subFrac + subStep
subStep = min ( 1.0_pReal - subFrac , num % stepIncreaseCryst * subStep )
todo = subStep > 0.0_pReal ! still time left to integrate on?
if ( todo ) then
2020-12-28 03:26:21 +05:30
subF0 = subF
2021-02-09 03:51:53 +05:30
subLp0 = phase_mechanical_Lp ( ph ) % data ( 1 : 3 , 1 : 3 , me )
subLi0 = phase_mechanical_Li ( ph ) % data ( 1 : 3 , 1 : 3 , me )
subFp0 = phase_mechanical_Fp ( ph ) % data ( 1 : 3 , 1 : 3 , me )
subFi0 = phase_mechanical_Fi ( ph ) % data ( 1 : 3 , 1 : 3 , me )
2020-12-29 04:43:49 +05:30
subState0 = plasticState ( ph ) % state ( : , me )
2021-02-13 23:11:30 +05:30
if ( damageState ( ph ) % sizeState > 0 ) &
2021-02-13 13:02:43 +05:30
damageState ( ph ) % subState0 ( : , me ) = damageState ( ph ) % state ( : , me )
2020-12-24 14:52:41 +05:30
endif
!--------------------------------------------------------------------------------------------------
! cut back (reduced time and restore)
else
subStep = num % subStepSizeCryst * subStep
2021-02-09 03:51:53 +05:30
phase_mechanical_Fp ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = subFp0
phase_mechanical_Fi ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = subFi0
phase_mechanical_S ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = phase_mechanical_S0 ( ph ) % data ( 1 : 3 , 1 : 3 , me ) ! why no subS0 ? is S0 of any use?
2020-12-29 04:43:49 +05:30
if ( subStep < 1.0_pReal ) then ! actual (not initial) cutback
2021-02-09 03:51:53 +05:30
phase_mechanical_Lp ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = subLp0
phase_mechanical_Li ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = subLi0
2020-12-24 14:52:41 +05:30
endif
2020-12-29 04:43:49 +05:30
plasticState ( ph ) % state ( : , me ) = subState0
2021-02-13 23:11:30 +05:30
if ( damageState ( ph ) % sizeState > 0 ) &
2021-02-13 13:02:43 +05:30
damageState ( ph ) % state ( : , me ) = damageState ( ph ) % subState0 ( : , me )
2020-12-24 14:52:41 +05:30
todo = subStep > num % subStepMinCryst ! still on track or already done (beyond repair)
endif
!--------------------------------------------------------------------------------------------------
! prepare for integration
if ( todo ) then
2020-12-28 03:26:21 +05:30
subF = subF0 &
2021-02-09 03:51:53 +05:30
+ subStep * ( phase_mechanical_F ( ph ) % data ( 1 : 3 , 1 : 3 , me ) - phase_mechanical_F0 ( ph ) % data ( 1 : 3 , 1 : 3 , me ) )
phase_mechanical_Fe ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = matmul ( subF , math_inv33 ( matmul ( phase_mechanical_Fi ( ph ) % data ( 1 : 3 , 1 : 3 , me ) , &
phase_mechanical_Fp ( ph ) % data ( 1 : 3 , 1 : 3 , me ) ) ) )
2020-12-29 04:43:49 +05:30
converged_ = . not . integrateState ( subF0 , subF , subFp0 , subFi0 , subState0 ( 1 : sizeDotState ) , subStep * dt , co , ip , el )
2021-01-08 04:02:54 +05:30
converged_ = converged_ . and . . not . integrateDamageState ( subStep * dt , co , ip , el )
2020-12-24 14:52:41 +05:30
endif
enddo cutbackLooping
end function crystallite_stress
2020-12-28 02:15:31 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Restore data after homog cutback.
!--------------------------------------------------------------------------------------------------
2021-02-09 03:51:53 +05:30
module subroutine mechanical_restore ( ce , includeL )
2020-12-28 02:15:31 +05:30
2021-01-24 14:09:40 +05:30
integer , intent ( in ) :: ce
2020-12-28 02:15:31 +05:30
logical , intent ( in ) :: &
includeL !< protect agains fake cutback
2020-12-29 05:14:42 +05:30
2020-12-28 02:15:31 +05:30
integer :: &
2020-12-29 05:14:42 +05:30
co , ph , me
2020-12-28 02:15:31 +05:30
2021-01-24 14:09:40 +05:30
do co = 1 , homogenization_Nconstituents ( material_homogenizationAt2 ( ce ) )
ph = material_phaseAt2 ( co , ce )
me = material_phaseMemberAt2 ( co , ce )
2020-12-28 02:15:31 +05:30
if ( includeL ) then
2021-02-09 03:51:53 +05:30
phase_mechanical_Lp ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = phase_mechanical_Lp0 ( ph ) % data ( 1 : 3 , 1 : 3 , me )
phase_mechanical_Li ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = phase_mechanical_Li0 ( ph ) % data ( 1 : 3 , 1 : 3 , me )
2020-12-28 02:15:31 +05:30
endif ! maybe protecting everything from overwriting makes more sense
2021-02-09 03:51:53 +05:30
phase_mechanical_Fp ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = phase_mechanical_Fp0 ( ph ) % data ( 1 : 3 , 1 : 3 , me )
phase_mechanical_Fi ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = phase_mechanical_Fi0 ( ph ) % data ( 1 : 3 , 1 : 3 , me )
phase_mechanical_S ( ph ) % data ( 1 : 3 , 1 : 3 , me ) = phase_mechanical_S0 ( ph ) % data ( 1 : 3 , 1 : 3 , me )
2020-12-28 02:15:31 +05:30
2021-01-25 18:33:49 +05:30
plasticState ( ph ) % state ( : , me ) = plasticState ( ph ) % State0 ( : , me )
2020-12-28 02:15:31 +05:30
enddo
2021-02-09 03:51:53 +05:30
end subroutine mechanical_restore
2020-12-28 02:15:31 +05:30
2020-12-30 02:01:22 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief Calculate tangent (dPdF).
!--------------------------------------------------------------------------------------------------
2021-02-09 03:51:53 +05:30
module function phase_mechanical_dPdF ( dt , co , ip , el ) result ( dPdF )
2020-12-30 02:01:22 +05:30
real ( pReal ) , intent ( in ) :: dt
integer , intent ( in ) :: &
co , & !< counter in constituent loop
ip , & !< counter in integration point loop
el !< counter in element loop
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) :: dPdF
integer :: &
o , &
p , ph , me
real ( pReal ) , dimension ( 3 , 3 ) :: devNull , &
invSubFp0 , invSubFi0 , invFp , invFi , &
temp_33_1 , temp_33_2 , temp_33_3
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) :: dSdFe , &
dSdF , &
dSdFi , &
dLidS , & ! tangent in lattice configuration
dLidFi , &
dLpdS , &
dLpdFi , &
dFidS , &
dFpinvdF , &
rhs_3333 , &
lhs_3333 , &
temp_3333
real ( pReal ) , dimension ( 9 , 9 ) :: temp_99
logical :: error
ph = material_phaseAt ( co , el )
me = material_phaseMemberAt ( co , ip , el )
2021-02-09 03:51:53 +05:30
call phase_hooke_SandItsTangents ( devNull , dSdFe , dSdFi , &
phase_mechanical_Fe ( ph ) % data ( 1 : 3 , 1 : 3 , me ) , &
phase_mechanical_Fi ( ph ) % data ( 1 : 3 , 1 : 3 , me ) , co , ip , el )
call phase_LiAndItsTangents ( devNull , dLidS , dLidFi , &
phase_mechanical_S ( ph ) % data ( 1 : 3 , 1 : 3 , me ) , &
phase_mechanical_Fi ( ph ) % data ( 1 : 3 , 1 : 3 , me ) , &
2021-02-13 11:45:57 +05:30
ph , me )
2020-12-30 02:01:22 +05:30
2021-02-09 03:51:53 +05:30
invFp = math_inv33 ( phase_mechanical_Fp ( ph ) % data ( 1 : 3 , 1 : 3 , me ) )
invFi = math_inv33 ( phase_mechanical_Fi ( ph ) % data ( 1 : 3 , 1 : 3 , me ) )
invSubFp0 = math_inv33 ( phase_mechanical_Fp0 ( ph ) % data ( 1 : 3 , 1 : 3 , me ) )
invSubFi0 = math_inv33 ( phase_mechanical_Fi0 ( ph ) % data ( 1 : 3 , 1 : 3 , me ) )
2020-12-30 02:01:22 +05:30
if ( sum ( abs ( dLidS ) ) < tol_math_check ) then
dFidS = 0.0_pReal
else
lhs_3333 = 0.0_pReal ; rhs_3333 = 0.0_pReal
do o = 1 , 3 ; do p = 1 , 3
lhs_3333 ( 1 : 3 , 1 : 3 , o , p ) = lhs_3333 ( 1 : 3 , 1 : 3 , o , p ) &
+ matmul ( invSubFi0 , dLidFi ( 1 : 3 , 1 : 3 , o , p ) ) * dt
lhs_3333 ( 1 : 3 , o , 1 : 3 , p ) = lhs_3333 ( 1 : 3 , o , 1 : 3 , p ) &
+ invFi * invFi ( p , o )
rhs_3333 ( 1 : 3 , 1 : 3 , o , p ) = rhs_3333 ( 1 : 3 , 1 : 3 , o , p ) &
- matmul ( invSubFi0 , dLidS ( 1 : 3 , 1 : 3 , o , p ) ) * dt
enddo ; enddo
call math_invert ( temp_99 , error , math_3333to99 ( lhs_3333 ) )
if ( error ) then
call IO_warning ( warning_ID = 600 , el = el , ip = ip , g = co , &
ext_msg = 'inversion error in analytic tangent calculation' )
dFidS = 0.0_pReal
else
dFidS = math_mul3333xx3333 ( math_99to3333 ( temp_99 ) , rhs_3333 )
endif
dLidS = math_mul3333xx3333 ( dLidFi , dFidS ) + dLidS
endif
2021-01-26 16:11:19 +05:30
call plastic_LpAndItsTangents ( devNull , dLpdS , dLpdFi , &
2021-02-09 03:51:53 +05:30
phase_mechanical_S ( ph ) % data ( 1 : 3 , 1 : 3 , me ) , &
2021-02-14 21:59:23 +05:30
phase_mechanical_Fi ( ph ) % data ( 1 : 3 , 1 : 3 , me ) , ph , me )
2020-12-30 02:01:22 +05:30
dLpdS = math_mul3333xx3333 ( dLpdFi , dFidS ) + dLpdS
!--------------------------------------------------------------------------------------------------
! calculate dSdF
temp_33_1 = transpose ( matmul ( invFp , invFi ) )
2021-02-09 03:51:53 +05:30
temp_33_2 = matmul ( phase_mechanical_F ( ph ) % data ( 1 : 3 , 1 : 3 , me ) , invSubFp0 )
temp_33_3 = matmul ( matmul ( phase_mechanical_F ( ph ) % data ( 1 : 3 , 1 : 3 , me ) , invFp ) , invSubFi0 )
2020-12-30 02:01:22 +05:30
do o = 1 , 3 ; do p = 1 , 3
rhs_3333 ( p , o , 1 : 3 , 1 : 3 ) = matmul ( dSdFe ( p , o , 1 : 3 , 1 : 3 ) , temp_33_1 )
temp_3333 ( 1 : 3 , 1 : 3 , p , o ) = matmul ( matmul ( temp_33_2 , dLpdS ( 1 : 3 , 1 : 3 , p , o ) ) , invFi ) &
+ matmul ( temp_33_3 , dLidS ( 1 : 3 , 1 : 3 , p , o ) )
enddo ; enddo
lhs_3333 = math_mul3333xx3333 ( dSdFe , temp_3333 ) * dt &
+ math_mul3333xx3333 ( dSdFi , dFidS )
call math_invert ( temp_99 , error , math_eye ( 9 ) + math_3333to99 ( lhs_3333 ) )
if ( error ) then
call IO_warning ( warning_ID = 600 , el = el , ip = ip , g = co , &
ext_msg = 'inversion error in analytic tangent calculation' )
dSdF = rhs_3333
else
dSdF = math_mul3333xx3333 ( math_99to3333 ( temp_99 ) , rhs_3333 )
endif
!--------------------------------------------------------------------------------------------------
! calculate dFpinvdF
temp_3333 = math_mul3333xx3333 ( dLpdS , dSdF )
do o = 1 , 3 ; do p = 1 , 3
dFpinvdF ( 1 : 3 , 1 : 3 , p , o ) = - matmul ( invSubFp0 , matmul ( temp_3333 ( 1 : 3 , 1 : 3 , p , o ) , invFi ) ) * dt
enddo ; enddo
!--------------------------------------------------------------------------------------------------
! assemble dPdF
2021-02-09 03:51:53 +05:30
temp_33_1 = matmul ( phase_mechanical_S ( ph ) % data ( 1 : 3 , 1 : 3 , me ) , transpose ( invFp ) )
temp_33_2 = matmul ( phase_mechanical_F ( ph ) % data ( 1 : 3 , 1 : 3 , me ) , invFp )
temp_33_3 = matmul ( temp_33_2 , phase_mechanical_S ( ph ) % data ( 1 : 3 , 1 : 3 , me ) )
2020-12-30 02:01:22 +05:30
dPdF = 0.0_pReal
do p = 1 , 3
dPdF ( p , 1 : 3 , p , 1 : 3 ) = transpose ( matmul ( invFp , temp_33_1 ) )
enddo
do o = 1 , 3 ; do p = 1 , 3
dPdF ( 1 : 3 , 1 : 3 , p , o ) = dPdF ( 1 : 3 , 1 : 3 , p , o ) &
2021-02-09 03:51:53 +05:30
+ matmul ( matmul ( phase_mechanical_F ( ph ) % data ( 1 : 3 , 1 : 3 , me ) , dFpinvdF ( 1 : 3 , 1 : 3 , p , o ) ) , temp_33_1 ) &
2020-12-30 02:01:22 +05:30
+ matmul ( matmul ( temp_33_2 , dSdF ( 1 : 3 , 1 : 3 , p , o ) ) , transpose ( invFp ) ) &
+ matmul ( temp_33_3 , transpose ( dFpinvdF ( 1 : 3 , 1 : 3 , p , o ) ) )
enddo ; enddo
2021-02-09 03:51:53 +05:30
end function phase_mechanical_dPdF
2020-12-30 02:01:22 +05:30
2020-12-30 04:44:48 +05:30
2021-02-09 03:51:53 +05:30
module subroutine mechanical_restartWrite ( groupHandle , ph )
2020-12-30 04:44:48 +05:30
integer ( HID_T ) , intent ( in ) :: groupHandle
integer , intent ( in ) :: ph
call HDF5_write ( groupHandle , plasticState ( ph ) % state , 'omega' )
2021-02-09 03:51:53 +05:30
call HDF5_write ( groupHandle , phase_mechanical_Fi ( ph ) % data , 'F_i' )
call HDF5_write ( groupHandle , phase_mechanical_Li ( ph ) % data , 'L_i' )
call HDF5_write ( groupHandle , phase_mechanical_Lp ( ph ) % data , 'L_p' )
call HDF5_write ( groupHandle , phase_mechanical_Fp ( ph ) % data , 'F_p' )
call HDF5_write ( groupHandle , phase_mechanical_S ( ph ) % data , 'S' )
call HDF5_write ( groupHandle , phase_mechanical_F ( ph ) % data , 'F' )
2020-12-30 04:44:48 +05:30
2021-02-09 03:51:53 +05:30
end subroutine mechanical_restartWrite
2020-12-30 04:44:48 +05:30
2021-02-09 03:51:53 +05:30
module subroutine mechanical_restartRead ( groupHandle , ph )
2020-12-30 04:44:48 +05:30
integer ( HID_T ) , intent ( in ) :: groupHandle
integer , intent ( in ) :: ph
call HDF5_read ( groupHandle , plasticState ( ph ) % state0 , 'omega' )
2021-02-09 03:51:53 +05:30
call HDF5_read ( groupHandle , phase_mechanical_Fi0 ( ph ) % data , 'F_i' )
call HDF5_read ( groupHandle , phase_mechanical_Li0 ( ph ) % data , 'L_i' )
call HDF5_read ( groupHandle , phase_mechanical_Lp0 ( ph ) % data , 'L_p' )
call HDF5_read ( groupHandle , phase_mechanical_Fp0 ( ph ) % data , 'F_p' )
call HDF5_read ( groupHandle , phase_mechanical_S0 ( ph ) % data , 'S' )
call HDF5_read ( groupHandle , phase_mechanical_F0 ( ph ) % data , 'F' )
2020-12-30 04:44:48 +05:30
2021-02-09 03:51:53 +05:30
end subroutine mechanical_restartRead
2020-12-30 04:44:48 +05:30
2020-12-30 14:24:06 +05:30
2020-12-30 18:27:37 +05:30
!----------------------------------------------------------------------------------------------
!< @brief Get first Piola-Kichhoff stress (for use by non-mech physics)
!----------------------------------------------------------------------------------------------
2021-02-09 03:51:53 +05:30
module function mechanical_S ( ph , me ) result ( S )
2020-12-30 14:24:06 +05:30
2020-12-30 17:15:08 +05:30
integer , intent ( in ) :: ph , me
2020-12-30 14:24:06 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: S
2021-02-09 03:51:53 +05:30
S = phase_mechanical_S ( ph ) % data ( 1 : 3 , 1 : 3 , me )
2020-12-30 14:24:06 +05:30
2021-02-09 03:51:53 +05:30
end function mechanical_S
2020-12-30 14:24:06 +05:30
2020-12-30 18:27:37 +05:30
!----------------------------------------------------------------------------------------------
!< @brief Get plastic velocity gradient (for use by non-mech physics)
!----------------------------------------------------------------------------------------------
2021-02-09 03:51:53 +05:30
module function mechanical_L_p ( ph , me ) result ( L_p )
2020-12-30 14:24:06 +05:30
2020-12-30 17:15:08 +05:30
integer , intent ( in ) :: ph , me
real ( pReal ) , dimension ( 3 , 3 ) :: L_p
2020-12-30 14:24:06 +05:30
2021-02-09 03:51:53 +05:30
L_p = phase_mechanical_Lp ( ph ) % data ( 1 : 3 , 1 : 3 , me )
2020-12-30 14:24:06 +05:30
2021-02-09 03:51:53 +05:30
end function mechanical_L_p
2020-12-30 14:24:06 +05:30
2020-12-30 18:27:37 +05:30
!----------------------------------------------------------------------------------------------
!< @brief Get deformation gradient (for use by homogenization)
!----------------------------------------------------------------------------------------------
2021-02-09 03:51:53 +05:30
module function phase_mechanical_getF ( co , ip , el ) result ( F )
2020-12-30 14:24:06 +05:30
integer , intent ( in ) :: co , ip , el
real ( pReal ) , dimension ( 3 , 3 ) :: F
2021-02-09 03:51:53 +05:30
F = phase_mechanical_F ( material_phaseAt ( co , el ) ) % data ( 1 : 3 , 1 : 3 , material_phaseMemberAt ( co , ip , el ) )
2020-12-30 14:24:06 +05:30
2021-02-09 03:51:53 +05:30
end function phase_mechanical_getF
2020-12-30 14:24:06 +05:30
2020-12-30 18:27:37 +05:30
!----------------------------------------------------------------------------------------------
!< @brief Get elastic deformation gradient (for use by non-mech physics)
!----------------------------------------------------------------------------------------------
2021-02-09 03:51:53 +05:30
module function mechanical_F_e ( ph , me ) result ( F_e )
2020-12-30 14:24:06 +05:30
2020-12-30 17:15:08 +05:30
integer , intent ( in ) :: ph , me
2020-12-30 14:24:06 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: F_e
2021-02-09 03:51:53 +05:30
F_e = phase_mechanical_Fe ( ph ) % data ( 1 : 3 , 1 : 3 , me )
2020-12-30 14:24:06 +05:30
2021-02-09 03:51:53 +05:30
end function mechanical_F_e
2020-12-30 14:24:06 +05:30
2020-12-30 15:33:13 +05:30
2020-12-30 18:27:37 +05:30
!----------------------------------------------------------------------------------------------
!< @brief Get second Piola-Kichhoff stress (for use by homogenization)
!----------------------------------------------------------------------------------------------
2021-02-09 03:51:53 +05:30
module function phase_mechanical_getP ( co , ip , el ) result ( P )
2020-12-30 15:33:13 +05:30
integer , intent ( in ) :: co , ip , el
real ( pReal ) , dimension ( 3 , 3 ) :: P
2021-02-09 03:51:53 +05:30
P = phase_mechanical_P ( material_phaseAt ( co , el ) ) % data ( 1 : 3 , 1 : 3 , material_phaseMemberAt ( co , ip , el ) )
2020-12-30 15:33:13 +05:30
2021-02-09 03:51:53 +05:30
end function phase_mechanical_getP
2020-12-30 15:33:13 +05:30
2020-12-30 14:24:06 +05:30
! setter for homogenization
2021-02-09 03:51:53 +05:30
module subroutine phase_mechanical_setF ( F , co , ip , el )
2020-12-30 14:24:06 +05:30
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: F
integer , intent ( in ) :: co , ip , el
2021-02-09 03:51:53 +05:30
phase_mechanical_F ( material_phaseAt ( co , el ) ) % data ( 1 : 3 , 1 : 3 , material_phaseMemberAt ( co , ip , el ) ) = F
2020-12-30 14:24:06 +05:30
2021-02-09 03:51:53 +05:30
end subroutine phase_mechanical_setF
2020-12-30 14:24:06 +05:30
2021-01-08 04:02:54 +05:30
2021-02-13 23:27:41 +05:30
end submodule mechanical
