2013-03-06 20:11:15 +05:30
!--------------------------------------------------------------------------------------------------
! $Id$
!--------------------------------------------------------------------------------------------------
!> @author Franz Roters, Max-Planck-Institut für Eisenforschung GmbH
!> @author Philip Eisenlohr, Max-Planck-Institut für Eisenforschung GmbH
!> @brief CPFEM engine
!--------------------------------------------------------------------------------------------------
2013-03-01 17:18:29 +05:30
module CPFEM
2013-03-06 20:11:15 +05:30
use prec , only : &
pReal , &
pInt
2009-06-15 18:41:21 +05:30
2013-03-01 17:18:29 +05:30
implicit none
2013-10-19 00:27:28 +05:30
private
2014-07-24 17:49:15 +05:30
#if defined(Marc4DAMASK) || defined(Abaqus)
2013-10-23 15:59:38 +05:30
real ( pReal ) , parameter , private :: &
CPFEM_odd_stress = 1e15_pReal , & !< return value for stress in case of ping pong dummy cycle
CPFEM_odd_jacobian = 1e50_pReal !< return value for jacobian in case of ping pong dummy cycle
real ( pReal ) , dimension ( : , : , : ) , allocatable , private :: &
CPFEM_cs !< Cauchy stress
real ( pReal ) , dimension ( : , : , : , : ) , allocatable , private :: &
CPFEM_dcsdE !< Cauchy stress tangent
real ( pReal ) , dimension ( : , : , : , : ) , allocatable , private :: &
CPFEM_dcsdE_knownGood !< known good tangent
2014-07-24 17:49:15 +05:30
#endif
2013-10-23 15:59:38 +05:30
logical , public , protected :: &
CPFEM_init_done = . false . , & !< remember whether init has been done already
CPFEM_calc_done = . false . !< remember whether first ip has already calced the results
integer ( pInt ) , parameter , public :: &
2013-03-01 17:18:29 +05:30
CPFEM_CALCRESULTS = 2_pInt ** 0_pInt , &
CPFEM_AGERESULTS = 2_pInt ** 1_pInt , &
CPFEM_BACKUPJACOBIAN = 2_pInt ** 2_pInt , &
CPFEM_RESTOREJACOBIAN = 2_pInt ** 3_pInt , &
2013-08-02 16:50:11 +05:30
CPFEM_COLLECT = 2_pInt ** 4_pInt
2013-10-23 15:59:38 +05:30
public :: &
CPFEM_general , &
2013-03-28 13:10:30 +05:30
CPFEM_initAll
2013-03-01 17:18:29 +05:30
contains
2009-03-04 19:31:36 +05:30
2010-11-03 20:09:18 +05:30
2013-03-06 20:11:15 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief call (thread safe) all module initializations
!--------------------------------------------------------------------------------------------------
2013-10-16 18:34:59 +05:30
subroutine CPFEM_initAll ( temperature , el , ip )
2013-03-06 20:11:15 +05:30
use prec , only : &
prec_init
use numerics , only : &
numerics_init
use debug , only : &
debug_init
use FEsolving , only : &
FE_init
use math , only : &
math_init
use mesh , only : &
mesh_init
use lattice , only : &
lattice_init
use material , only : &
material_init
use constitutive , only : &
constitutive_init
use crystallite , only : &
crystallite_init
use homogenization , only : &
homogenization_init
use IO , only : &
IO_init
2013-03-01 17:18:29 +05:30
use DAMASK_interface
2014-06-06 06:08:29 +05:30
#ifdef FEM
use FEZoo , only : &
FEZoo_init
#endif
2014-06-25 04:49:21 +05:30
use constitutive_thermal , only : &
constitutive_thermal_init
use constitutive_damage , only : &
constitutive_damage_init
2013-03-01 17:18:29 +05:30
implicit none
2013-10-16 18:34:59 +05:30
integer ( pInt ) , intent ( in ) :: el , & ! FE el number
ip ! FE integration point number
real ( pReal ) , intent ( in ) :: temperature ! temperature
2013-03-01 17:18:29 +05:30
2014-05-14 19:27:25 +05:30
!$OMP CRITICAL (init)
if ( . not . CPFEM_init_done ) then
2014-07-24 17:49:15 +05:30
#if defined(Spectral) || defined(FEM)
call DAMASK_interface_init ! Spectral and FEM interface to commandline
2012-07-30 21:21:48 +05:30
#endif
2013-03-01 17:18:29 +05:30
call prec_init
call IO_init
2014-06-06 06:08:29 +05:30
#ifdef FEM
2014-07-24 17:49:15 +05:30
call FEZoo_init
2014-06-06 06:08:29 +05:30
#endif
2013-03-01 17:18:29 +05:30
call numerics_init
call debug_init
call math_init
call FE_init
2014-05-14 19:27:25 +05:30
call mesh_init ( ip , el ) ! pass on coordinates to alter calcMode of first ip
2013-03-01 17:18:29 +05:30
call lattice_init
call material_init
call constitutive_init
2014-06-25 04:49:21 +05:30
call constitutive_thermal_init
call constitutive_damage_init
2013-10-16 18:34:59 +05:30
call crystallite_init ( temperature ) ! (have to) use temperature of first ip for whole model
call homogenization_init
2013-03-01 17:18:29 +05:30
call CPFEM_init
2014-07-24 17:49:15 +05:30
#if defined(Marc4DAMASK) || defined(Abaqus)
call DAMASK_interface_init ! Spectral solver and FEM init is already done
2012-06-15 21:40:21 +05:30
#endif
2013-03-01 17:18:29 +05:30
CPFEM_init_done = . true .
endif
2014-05-14 19:27:25 +05:30
!$OMP END CRITICAL (init)
2010-11-03 20:09:18 +05:30
2012-03-09 01:55:28 +05:30
end subroutine CPFEM_initAll
2010-11-03 20:09:18 +05:30
2013-03-06 20:11:15 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief allocate the arrays defined in module CPFEM and initialize them
!--------------------------------------------------------------------------------------------------
2012-03-09 01:55:28 +05:30
subroutine CPFEM_init
2013-03-01 17:18:29 +05:30
use , intrinsic :: iso_fortran_env ! to get compiler_version and compiler_options (at least for gfortran 4.6 at the moment)
use prec , only : &
pInt
use IO , only : &
2013-09-18 19:37:55 +05:30
IO_read_realFile , &
IO_read_intFile , &
2013-03-01 17:18:29 +05:30
IO_timeStamp , &
IO_error
use numerics , only : &
2013-08-02 16:50:11 +05:30
DAMASK_NumThreadsInt
2013-03-01 17:18:29 +05:30
use debug , only : &
debug_level , &
debug_CPFEM , &
debug_levelBasic , &
debug_levelExtensive
use FEsolving , only : &
symmetricSolver , &
restartRead , &
modelName
use mesh , only : &
mesh_NcpElems , &
mesh_maxNips
use material , only : &
homogenization_maxNgrains , &
2014-07-24 14:08:52 +05:30
material_phase , &
2014-08-21 23:18:20 +05:30
#ifdef NEWSTATE
homogState , &
mappingHomogenization , &
#endif
2014-07-24 14:08:52 +05:30
phase_plasticity , &
plasticState
2013-03-01 17:18:29 +05:30
use crystallite , only : &
crystallite_F0 , &
crystallite_Fp0 , &
crystallite_Lp0 , &
crystallite_dPdF0 , &
crystallite_Tstar0_v , &
crystallite_localPlasticity
use homogenization , only : &
2014-08-21 23:18:20 +05:30
#ifndef NEWSTATE
homogenization_state0 , &
#endif
homogenization_sizeState
2013-03-01 17:18:29 +05:30
implicit none
2014-07-24 14:08:52 +05:30
integer ( pInt ) :: i , j , k , l , m , ph
2013-03-01 17:18:29 +05:30
2013-03-27 17:58:55 +05:30
write ( 6 , '(/,a)' ) ' <<<+- CPFEM init -+>>>'
write ( 6 , '(a)' ) ' $Id$'
write ( 6 , '(a15,a)' ) ' Current time: ' , IO_timeStamp ( )
#include "compilation_info.f90"
2014-07-24 17:49:15 +05:30
#if defined(Marc4DAMASK) || defined(Abaqus)
! initialize stress and jacobian to zero
2013-03-01 17:18:29 +05:30
allocate ( CPFEM_cs ( 6 , mesh_maxNips , mesh_NcpElems ) ) ; CPFEM_cs = 0.0_pReal
allocate ( CPFEM_dcsdE ( 6 , 6 , mesh_maxNips , mesh_NcpElems ) ) ; CPFEM_dcsdE = 0.0_pReal
allocate ( CPFEM_dcsdE_knownGood ( 6 , 6 , mesh_maxNips , mesh_NcpElems ) ) ; CPFEM_dcsdE_knownGood = 0.0_pReal
2014-07-24 17:49:15 +05:30
#endif
2013-03-01 17:18:29 +05:30
! *** restore the last converged values of each essential variable from the binary file
if ( restartRead ) then
if ( iand ( debug_level ( debug_CPFEM ) , debug_levelExtensive ) / = 0_pInt ) then
2014-05-16 19:31:27 +05:30
write ( 6 , '(a)' ) '<< CPFEM >> restored state variables of last converged step from binary files'
flush ( 6 )
2013-03-01 17:18:29 +05:30
endif
2013-09-18 19:37:55 +05:30
call IO_read_intFile ( 777 , 'recordedPhase' , modelName , size ( material_phase ) )
2013-03-01 17:18:29 +05:30
read ( 777 , rec = 1 ) material_phase
close ( 777 )
2013-09-18 19:37:55 +05:30
call IO_read_realFile ( 777 , 'convergedF' , modelName , size ( crystallite_F0 ) )
2013-03-01 17:18:29 +05:30
read ( 777 , rec = 1 ) crystallite_F0
close ( 777 )
2013-09-18 19:37:55 +05:30
call IO_read_realFile ( 777 , 'convergedFp' , modelName , size ( crystallite_Fp0 ) )
2013-03-01 17:18:29 +05:30
read ( 777 , rec = 1 ) crystallite_Fp0
close ( 777 )
2013-09-18 19:37:55 +05:30
call IO_read_realFile ( 777 , 'convergedLp' , modelName , size ( crystallite_Lp0 ) )
2013-03-01 17:18:29 +05:30
read ( 777 , rec = 1 ) crystallite_Lp0
close ( 777 )
2013-09-18 19:37:55 +05:30
call IO_read_realFile ( 777 , 'convergeddPdF' , modelName , size ( crystallite_dPdF0 ) )
2013-03-01 17:18:29 +05:30
read ( 777 , rec = 1 ) crystallite_dPdF0
close ( 777 )
2013-09-18 19:37:55 +05:30
call IO_read_realFile ( 777 , 'convergedTstar' , modelName , size ( crystallite_Tstar0_v ) )
2013-03-01 17:18:29 +05:30
read ( 777 , rec = 1 ) crystallite_Tstar0_v
close ( 777 )
2014-07-02 17:57:39 +05:30
2013-09-18 19:37:55 +05:30
call IO_read_realFile ( 777 , 'convergedStateConst' , modelName )
2013-03-01 17:18:29 +05:30
m = 0_pInt
2014-07-24 14:08:52 +05:30
readInstances : do ph = 1_pInt , size ( phase_plasticity )
do k = 1_pInt , plasticState ( ph ) % sizeState
2014-08-04 23:20:01 +05:30
do l = 1 , size ( plasticState ( ph ) % state0 ( 1 , : ) )
2014-07-24 14:08:52 +05:30
m = m + 1_pInt
read ( 777 , rec = m ) plasticState ( ph ) % state0 ( k , l )
enddo ; enddo
enddo readInstances
2013-03-01 17:18:29 +05:30
close ( 777 )
2013-09-18 19:37:55 +05:30
call IO_read_realFile ( 777 , 'convergedStateHomog' , modelName )
2013-03-01 17:18:29 +05:30
m = 0_pInt
do k = 1 , mesh_NcpElems ; do j = 1 , mesh_maxNips
do l = 1 , homogenization_sizeState ( j , k )
m = m + 1_pInt
2014-08-21 23:18:20 +05:30
#ifdef NEWSTATE
read ( 777 , rec = m ) homogState ( mappingHomogenization ( 2 , j , k ) ) % state0 ( l , mappingHomogenization ( 1 , j , k ) )
#else
2013-03-01 17:18:29 +05:30
read ( 777 , rec = m ) homogenization_state0 ( j , k ) % p ( l )
2014-08-21 23:18:20 +05:30
#endif
2013-03-01 17:18:29 +05:30
enddo
2014-08-04 23:20:01 +05:30
2013-03-01 17:18:29 +05:30
enddo ; enddo
close ( 777 )
2014-07-24 17:49:15 +05:30
#if defined(Marc4DAMASK) || defined(Abaqus)
2013-09-18 19:37:55 +05:30
call IO_read_realFile ( 777 , 'convergeddcsdE' , modelName , size ( CPFEM_dcsdE ) )
2013-03-01 17:18:29 +05:30
read ( 777 , rec = 1 ) CPFEM_dcsdE
close ( 777 )
2014-07-24 17:49:15 +05:30
#endif
2013-03-01 17:18:29 +05:30
restartRead = . false .
endif
2014-07-24 17:49:15 +05:30
#if defined(Marc4DAMASK) || defined(Abaqus)
2013-03-01 17:18:29 +05:30
if ( iand ( debug_level ( debug_CPFEM ) , debug_levelBasic ) / = 0 ) then
2013-10-16 18:34:59 +05:30
write ( 6 , '(a32,1x,6(i8,1x))' ) 'CPFEM_cs: ' , shape ( CPFEM_cs )
write ( 6 , '(a32,1x,6(i8,1x))' ) 'CPFEM_dcsdE: ' , shape ( CPFEM_dcsdE )
write ( 6 , '(a32,1x,6(i8,1x),/)' ) 'CPFEM_dcsdE_knownGood: ' , shape ( CPFEM_dcsdE_knownGood )
2014-07-23 18:56:05 +05:30
write ( 6 , '(a32,l1)' ) 'symmetricSolver: ' , symmetricSolver
2013-03-01 17:18:29 +05:30
endif
2014-07-24 17:49:15 +05:30
#endif
2013-03-01 17:18:29 +05:30
flush ( 6 )
2009-06-15 18:41:21 +05:30
2012-03-09 01:55:28 +05:30
end subroutine CPFEM_init
2009-06-15 18:41:21 +05:30
2013-03-06 20:11:15 +05:30
!--------------------------------------------------------------------------------------------------
!> @brief perform initialization at first call, update variables and call the actual material model
!--------------------------------------------------------------------------------------------------
2014-07-24 17:49:15 +05:30
#if defined(Marc4DAMASK) || defined(Abaqus)
2013-10-16 18:34:59 +05:30
subroutine CPFEM_general ( mode , parallelExecution , ffn , ffn1 , temperature , dt , elFE , ip , cauchyStress , jacobian )
2014-07-24 17:49:15 +05:30
#else
subroutine CPFEM_general ( mode , ffn , ffn1 , temperature , dt , elFE , ip )
#endif
2013-10-16 18:34:59 +05:30
use numerics , only : &
defgradTolerance , &
iJacoStiffness
use debug , only : &
debug_level , &
debug_CPFEM , &
debug_levelBasic , &
debug_levelExtensive , &
debug_levelSelective , &
debug_e , &
debug_i , &
debug_stressMaxLocation , &
debug_stressMinLocation , &
debug_jacobianMaxLocation , &
debug_jacobianMinLocation , &
debug_stressMax , &
debug_stressMin , &
debug_jacobianMax , &
debug_jacobianMin
use FEsolving , only : &
outdatedFFN1 , &
terminallyIll , &
cycleCounter , &
theInc , &
theTime , &
theDelta , &
FEsolving_execElem , &
FEsolving_execIP , &
restartWrite
use math , only : &
math_identity2nd , &
math_mul33x33 , &
math_det33 , &
math_transpose33 , &
math_I3 , &
math_Mandel3333to66 , &
math_Mandel66to3333 , &
math_Mandel33to6 , &
math_Mandel6to33
use mesh , only : &
mesh_FEasCP , &
mesh_NcpElems , &
mesh_maxNips , &
mesh_element
use material , only : &
homogenization_maxNgrains , &
microstructure_elemhomo , &
2014-07-02 17:57:39 +05:30
plasticState , &
damageState , &
2014-08-21 23:18:20 +05:30
#ifdef NEWSTATE
homogState , &
mappingHomogenization , &
#endif
2014-07-02 17:57:39 +05:30
thermalState , &
mappingConstitutive , &
2014-07-24 14:08:52 +05:30
material_phase , &
phase_plasticity
2013-10-16 18:34:59 +05:30
use crystallite , only : &
crystallite_partionedF , &
crystallite_F0 , &
crystallite_Fp0 , &
crystallite_Fp , &
crystallite_Lp0 , &
crystallite_Lp , &
crystallite_dPdF0 , &
crystallite_dPdF , &
crystallite_Tstar0_v , &
crystallite_Tstar_v , &
crystallite_temperature
use homogenization , only : &
homogenization_sizeState , &
2014-08-21 23:18:20 +05:30
#ifndef NEWSTATE
2013-10-16 18:34:59 +05:30
homogenization_state , &
homogenization_state0 , &
2014-08-21 23:18:20 +05:30
#endif
2013-10-16 18:34:59 +05:30
materialpoint_F , &
materialpoint_F0 , &
materialpoint_P , &
materialpoint_dPdF , &
materialpoint_results , &
materialpoint_sizeResults , &
materialpoint_stressAndItsTangent , &
materialpoint_postResults
use IO , only : &
IO_write_jobRealFile , &
IO_warning
2013-03-01 17:18:29 +05:30
use DAMASK_interface
2013-10-16 18:34:59 +05:30
2013-03-01 17:18:29 +05:30
implicit none
2013-10-16 18:34:59 +05:30
integer ( pInt ) , intent ( in ) :: elFE , & !< FE element number
ip !< integration point number
real ( pReal ) , intent ( in ) :: temperature !< temperature
2013-04-09 15:38:00 +05:30
real ( pReal ) , intent ( in ) :: dt !< time increment
real ( pReal ) , dimension ( 3 , 3 ) , intent ( in ) :: ffn , & !< deformation gradient for t=t0
ffn1 !< deformation gradient for t=t1
integer ( pInt ) , intent ( in ) :: mode !< computation mode 1: regular computation plus aging of results
2014-07-24 17:49:15 +05:30
#if defined(Marc4DAMASK) || defined(Abaqus)
2013-08-02 16:50:11 +05:30
logical , intent ( in ) :: parallelExecution !< flag indicating parallel computation of requested IPs
2014-07-24 17:49:15 +05:30
real ( pReal ) , dimension ( 6 ) , intent ( out ) :: cauchyStress !< stress vector in Mandel notation
real ( pReal ) , dimension ( 6 , 6 ) , intent ( out ) :: jacobian !< jacobian in Mandel notation (Consistent tangent dcs/dE)
2013-04-09 15:38:00 +05:30
2014-07-24 17:49:15 +05:30
real ( pReal ) J_inverse , & ! inverse of Jacobian
2013-03-01 17:18:29 +05:30
rnd
2014-07-24 17:49:15 +05:30
real ( pReal ) , dimension ( 3 , 3 ) :: Kirchhoff , & ! Piola-Kirchhoff stress in Matrix notation
cauchyStress33 ! stress vector in Matrix notation
2013-03-01 17:18:29 +05:30
real ( pReal ) , dimension ( 3 , 3 , 3 , 3 ) :: H_sym , &
H , &
2014-07-24 17:49:15 +05:30
jacobian3333 ! jacobian in Matrix notation
#else
logical , parameter :: parallelExecution = . true .
#endif
integer ( pInt ) elCP , & ! crystal plasticity element number
2014-07-24 14:08:52 +05:30
i , j , k , l , m , n , ph
2014-07-24 17:49:15 +05:30
logical updateJaco ! flag indicating if JAcobian has to be updated
2013-03-01 17:18:29 +05:30
2014-07-24 17:49:15 +05:30
#if defined(Marc4DAMASK) || defined(Abaqus)
2013-10-16 18:34:59 +05:30
elCP = mesh_FEasCP ( 'elem' , elFE )
2014-07-24 17:49:15 +05:30
#else
elCP = elFE
#endif
2013-03-01 17:18:29 +05:30
2013-08-02 16:50:11 +05:30
if ( iand ( debug_level ( debug_CPFEM ) , debug_levelBasic ) / = 0_pInt &
2013-10-16 18:34:59 +05:30
. and . elCP == debug_e . and . ip == debug_i ) then
2014-05-16 19:31:27 +05:30
write ( 6 , '(/,a)' ) '#############################################'
write ( 6 , '(a1,a22,1x,i8,a13)' ) '#' , 'element' , elCP , '#'
write ( 6 , '(a1,a22,1x,i8,a13)' ) '#' , 'ip' , ip , '#'
write ( 6 , '(a1,a22,1x,f15.7,a6)' ) '#' , 'theTime' , theTime , '#'
write ( 6 , '(a1,a22,1x,f15.7,a6)' ) '#' , 'theDelta' , theDelta , '#'
write ( 6 , '(a1,a22,1x,i8,a13)' ) '#' , 'theInc' , theInc , '#'
write ( 6 , '(a1,a22,1x,i8,a13)' ) '#' , 'cycleCounter' , cycleCounter , '#'
write ( 6 , '(a1,a22,1x,i8,a13)' ) '#' , 'computationMode' , mode , '#'
write ( 6 , '(a,/)' ) '#############################################' ; flush ( 6 )
2013-03-01 17:18:29 +05:30
endif
2013-04-16 22:37:27 +05:30
2014-07-24 17:49:15 +05:30
#if defined(Marc4DAMASK) || defined(Abaqus)
2013-08-02 17:06:51 +05:30
!*** backup or restore jacobian
if ( iand ( mode , CPFEM_BACKUPJACOBIAN ) / = 0_pInt ) &
CPFEM_dcsde_knownGood = CPFEM_dcsde
if ( iand ( mode , CPFEM_RESTOREJACOBIAN ) / = 0_pInt ) &
CPFEM_dcsde = CPFEM_dcsde_knownGood
2014-07-24 17:49:15 +05:30
#endif
2013-08-02 17:06:51 +05:30
!*** age results and write restart data if requested
2013-04-16 22:37:27 +05:30
if ( iand ( mode , CPFEM_AGERESULTS ) / = 0_pInt ) then
crystallite_F0 = crystallite_partionedF ! crystallite deformation (_subF is perturbed...)
crystallite_Fp0 = crystallite_Fp ! crystallite plastic deformation
crystallite_Lp0 = crystallite_Lp ! crystallite plastic velocity
crystallite_dPdF0 = crystallite_dPdF ! crystallite stiffness
crystallite_Tstar0_v = crystallite_Tstar_v ! crystallite 2nd Piola Kirchhoff stress
2014-05-12 18:30:37 +05:30
2014-07-24 17:49:15 +05:30
forall ( i = 1 : size ( plasticState ) ) plasticState ( i ) % state0 = plasticState ( i ) % state ! copy state in this lenghty way because: A component cannot be an array if the encompassing structure is an array
forall ( i = 1 : size ( damageState ) ) damageState ( i ) % state0 = damageState ( i ) % state ! copy state in this lenghty way because: A component cannot be an array if the encompassing structure is an array
forall ( i = 1 : size ( thermalState ) ) thermalState ( i ) % state0 = thermalState ( i ) % state ! copy state in this lenghty way because: A component cannot be an array if the encompassing structure is an array
2014-08-04 23:20:01 +05:30
if ( iand ( debug_level ( debug_CPFEM ) , debug_levelBasic ) / = 0_pInt ) then
2014-05-16 19:31:27 +05:30
write ( 6 , '(a)' ) '<< CPFEM >> aging states'
if ( debug_e < = mesh_NcpElems . and . debug_i < = mesh_maxNips ) then
write ( 6 , '(a,1x,i8,1x,i2,1x,i4,/,(12x,6(e20.8,1x)),/)' ) &
'<< CPFEM >> aged state of elFE ip grain' , debug_e , debug_i , 1 , &
plasticState ( mappingConstitutive ( 2 , 1 , debug_i , debug_e ) ) % state ( : , mappingConstitutive ( 1 , 1 , debug_i , debug_e ) )
2014-05-12 18:30:37 +05:30
endif
endif
2014-07-02 17:57:39 +05:30
2013-08-01 21:40:56 +05:30
!$OMP PARALLEL DO
do k = 1 , mesh_NcpElems
do j = 1 , mesh_maxNips
if ( homogenization_sizeState ( j , k ) > 0_pInt ) &
2014-08-21 23:18:20 +05:30
#ifdef NEWSTATE
homogState ( mappingHomogenization ( 2 , j , k ) ) % state0 ( : , mappingHomogenization ( 1 , j , k ) ) = &
homogState ( mappingHomogenization ( 2 , j , k ) ) % state ( : , mappingHomogenization ( 1 , j , k ) ) ! internal state of homogenization scheme
#else
2013-08-01 21:40:56 +05:30
homogenization_state0 ( j , k ) % p = homogenization_state ( j , k ) % p ! internal state of homogenization scheme
2014-08-21 23:18:20 +05:30
#endif
2013-08-01 21:40:56 +05:30
enddo
enddo
!$OMP END PARALLEL DO
! * dump the last converged values of each essential variable to a binary file
if ( restartWrite ) then
2014-08-04 23:20:01 +05:30
if ( iand ( debug_level ( debug_CPFEM ) , debug_levelBasic ) / = 0_pInt ) &
2013-08-01 21:40:56 +05:30
write ( 6 , '(a)' ) '<< CPFEM >> writing state variables of last converged step to binary files'
2013-09-18 19:37:55 +05:30
call IO_write_jobRealFile ( 777 , 'recordedPhase' , size ( material_phase ) )
2013-08-01 21:40:56 +05:30
write ( 777 , rec = 1 ) material_phase
close ( 777 )
2013-03-01 17:18:29 +05:30
2013-09-18 19:37:55 +05:30
call IO_write_jobRealFile ( 777 , 'convergedF' , size ( crystallite_F0 ) )
2013-08-01 21:40:56 +05:30
write ( 777 , rec = 1 ) crystallite_F0
close ( 777 )
2013-03-01 17:18:29 +05:30
2013-09-18 19:37:55 +05:30
call IO_write_jobRealFile ( 777 , 'convergedFp' , size ( crystallite_Fp0 ) )
2013-08-01 21:40:56 +05:30
write ( 777 , rec = 1 ) crystallite_Fp0
close ( 777 )
2013-03-01 17:18:29 +05:30
2013-09-18 19:37:55 +05:30
call IO_write_jobRealFile ( 777 , 'convergedLp' , size ( crystallite_Lp0 ) )
2013-08-01 21:40:56 +05:30
write ( 777 , rec = 1 ) crystallite_Lp0
close ( 777 )
2013-03-01 17:18:29 +05:30
2013-09-18 19:37:55 +05:30
call IO_write_jobRealFile ( 777 , 'convergeddPdF' , size ( crystallite_dPdF0 ) )
2013-08-01 21:40:56 +05:30
write ( 777 , rec = 1 ) crystallite_dPdF0
close ( 777 )
2013-03-01 17:18:29 +05:30
2013-09-18 19:37:55 +05:30
call IO_write_jobRealFile ( 777 , 'convergedTstar' , size ( crystallite_Tstar0_v ) )
2013-08-01 21:40:56 +05:30
write ( 777 , rec = 1 ) crystallite_Tstar0_v
close ( 777 )
2014-07-24 14:08:52 +05:30
2013-09-18 19:37:55 +05:30
call IO_write_jobRealFile ( 777 , 'convergedStateConst' )
2013-08-01 21:40:56 +05:30
m = 0_pInt
2014-07-24 14:08:52 +05:30
writeInstances : do ph = 1_pInt , size ( phase_plasticity )
do k = 1_pInt , plasticState ( ph ) % sizeState
2014-08-04 23:20:01 +05:30
do l = 1 , size ( plasticState ( ph ) % state0 ( 1 , : ) )
2014-07-24 14:08:52 +05:30
m = m + 1_pInt
write ( 777 , rec = m ) plasticState ( ph ) % state0 ( k , l )
enddo ; enddo
enddo writeInstances
2013-08-01 21:40:56 +05:30
close ( 777 )
2014-07-23 18:56:05 +05:30
2013-09-18 19:37:55 +05:30
call IO_write_jobRealFile ( 777 , 'convergedStateHomog' )
2013-08-01 21:40:56 +05:30
m = 0_pInt
do k = 1 , mesh_NcpElems ; do j = 1 , mesh_maxNips
do l = 1 , homogenization_sizeState ( j , k )
m = m + 1_pInt
2014-08-21 23:18:20 +05:30
#ifdef NEWSTATE
write ( 777 , rec = m ) homogState ( mappingHomogenization ( 2 , j , k ) ) % state0 ( l , mappingHomogenization ( 1 , j , k ) )
#else
2013-08-01 21:40:56 +05:30
write ( 777 , rec = m ) homogenization_state0 ( j , k ) % p ( l )
2014-08-21 23:18:20 +05:30
#endif
2013-08-01 21:40:56 +05:30
enddo
enddo ; enddo
close ( 777 )
2014-07-24 17:49:15 +05:30
#if defined(Marc4DAMASK) || defined(Abaqus)
2013-09-18 19:37:55 +05:30
call IO_write_jobRealFile ( 777 , 'convergeddcsdE' , size ( CPFEM_dcsdE ) )
2013-08-01 21:40:56 +05:30
write ( 777 , rec = 1 ) CPFEM_dcsdE
close ( 777 )
2014-07-24 17:49:15 +05:30
#endif
2013-08-01 21:40:56 +05:30
endif
endif
2013-08-02 17:06:51 +05:30
!*** collection of FEM input with returning of randomize odd stress and jacobian
!* In case that no parallel execution is required, there is no need to collect FEM input
if ( . not . parallelExecution ) then
2013-10-16 18:34:59 +05:30
crystallite_temperature ( ip , elCP ) = temperature
materialpoint_F0 ( 1 : 3 , 1 : 3 , ip , elCP ) = ffn
materialpoint_F ( 1 : 3 , 1 : 3 , ip , elCP ) = ffn1
2013-08-02 17:06:51 +05:30
elseif ( iand ( mode , CPFEM_COLLECT ) / = 0_pInt ) then
2014-07-24 17:49:15 +05:30
#if defined(Marc4DAMASK) || defined(Abaqus)
2013-08-02 17:06:51 +05:30
call random_number ( rnd )
if ( rnd < 0.5_pReal ) rnd = rnd - 1.0_pReal
2014-07-24 17:49:15 +05:30
CPFEM_cs ( 1 : 6 , ip , elCP ) = rnd * CPFEM_odd_stress
CPFEM_dcsde ( 1 : 6 , 1 : 6 , ip , elCP ) = CPFEM_odd_jacobian * math_identity2nd ( 6 )
#endif
2013-10-16 18:34:59 +05:30
crystallite_temperature ( ip , elCP ) = temperature
materialpoint_F0 ( 1 : 3 , 1 : 3 , ip , elCP ) = ffn
materialpoint_F ( 1 : 3 , 1 : 3 , ip , elCP ) = ffn1
2013-08-02 17:06:51 +05:30
CPFEM_calc_done = . false .
2013-08-01 21:40:56 +05:30
endif
2013-08-02 17:06:51 +05:30
!*** calculation of stress and jacobian
2013-08-01 21:40:56 +05:30
if ( iand ( mode , CPFEM_CALCRESULTS ) / = 0_pInt ) then
!*** deformation gradient outdated or any actual deformation gradient differs more than relevantStrain from the stored one
2013-10-16 18:34:59 +05:30
if ( terminallyIll . or . outdatedFFN1 . or . any ( abs ( ffn1 - materialpoint_F ( 1 : 3 , 1 : 3 , ip , elCP ) ) > defgradTolerance ) ) then
if ( any ( abs ( ffn1 - materialpoint_F ( 1 : 3 , 1 : 3 , ip , elCP ) ) > defgradTolerance ) ) then
2013-08-01 21:40:56 +05:30
if ( iand ( debug_level ( debug_CPFEM ) , debug_levelExtensive ) / = 0_pInt ) then
2013-10-16 18:34:59 +05:30
write ( 6 , '(a,1x,i8,1x,i2)' ) '<< CPFEM >> OUTDATED at elFE ip' , elCP , ip
2013-10-19 00:27:28 +05:30
write ( 6 , '(a,/,3(12x,3(f10.6,1x),/))' ) '<< CPFEM >> FFN1 old:' , &
math_transpose33 ( materialpoint_F ( 1 : 3 , 1 : 3 , ip , elCP ) )
2013-08-01 21:40:56 +05:30
write ( 6 , '(a,/,3(12x,3(f10.6,1x),/))' ) '<< CPFEM >> FFN1 now:' , math_transpose33 ( ffn1 )
endif
outdatedFFN1 = . true .
endif
2014-07-24 17:49:15 +05:30
#if defined(Marc4DAMASK) || defined(Abaqus)
2013-08-01 21:40:56 +05:30
call random_number ( rnd )
if ( rnd < 0.5_pReal ) rnd = rnd - 1.0_pReal
2013-10-16 18:34:59 +05:30
CPFEM_cs ( 1 : 6 , ip , elCP ) = rnd * CPFEM_odd_stress
CPFEM_dcsde ( 1 : 6 , 1 : 6 , ip , elCP ) = CPFEM_odd_jacobian * math_identity2nd ( 6 )
2014-07-24 17:49:15 +05:30
#endif
2013-08-01 21:40:56 +05:30
!*** deformation gradient is not outdated
else
2014-07-24 17:49:15 +05:30
updateJaco = mod ( cycleCounter , iJacoStiffness ) == 0
2013-10-16 18:34:59 +05:30
!* no parallel computation, so we use just one single elFE and ip for computation
2013-08-01 21:40:56 +05:30
if ( . not . parallelExecution ) then
2013-10-16 18:34:59 +05:30
FEsolving_execElem ( 1 ) = elCP
FEsolving_execElem ( 2 ) = elCP
2014-03-12 13:03:51 +05:30
if ( . not . microstructure_elemhomo ( mesh_element ( 4 , elCP ) ) . or . & ! calculate unless homogeneous
( microstructure_elemhomo ( mesh_element ( 4 , elCP ) ) . and . ip == 1_pInt ) ) then ! and then only first ip
2013-10-16 18:34:59 +05:30
FEsolving_execIP ( 1 , elCP ) = ip
FEsolving_execIP ( 2 , elCP ) = ip
2014-05-16 19:31:27 +05:30
if ( iand ( debug_level ( debug_CPFEM ) , debug_levelExtensive ) / = 0_pInt ) &
write ( 6 , '(a,i8,1x,i2)' ) '<< CPFEM >> calculation for elFE ip ' , elCP , ip
2014-03-12 13:03:51 +05:30
call materialpoint_stressAndItsTangent ( updateJaco , dt ) ! calculate stress and its tangent
call materialpoint_postResults ( )
2013-03-01 17:18:29 +05:30
endif
2013-08-01 21:40:56 +05:30
!* parallel computation and calulation not yet done
elseif ( . not . CPFEM_calc_done ) then
2014-05-16 19:31:27 +05:30
if ( iand ( debug_level ( debug_CPFEM ) , debug_levelExtensive ) / = 0_pInt ) &
write ( 6 , '(a,i8,a,i8)' ) '<< CPFEM >> calculation for elements ' , FEsolving_execElem ( 1 ) , &
' to ' , FEsolving_execElem ( 2 )
2014-03-12 13:03:51 +05:30
call materialpoint_stressAndItsTangent ( updateJaco , dt ) ! calculate stress and its tangent (parallel execution inside)
call materialpoint_postResults ( )
2013-08-01 21:40:56 +05:30
CPFEM_calc_done = . true .
2013-03-01 17:18:29 +05:30
endif
2013-08-01 21:40:56 +05:30
!* map stress and stiffness (or return odd values if terminally ill)
2014-08-04 23:20:01 +05:30
#if defined(Marc4DAMASK) || defined(Abaqus)
2013-08-01 21:40:56 +05:30
if ( terminallyIll ) then
2013-03-01 17:18:29 +05:30
call random_number ( rnd )
if ( rnd < 0.5_pReal ) rnd = rnd - 1.0_pReal
2013-10-16 18:34:59 +05:30
CPFEM_cs ( 1 : 6 , ip , elCP ) = rnd * CPFEM_odd_stress
CPFEM_dcsde ( 1 : 6 , 1 : 6 , ip , elCP ) = CPFEM_odd_jacobian * math_identity2nd ( 6 )
2013-08-01 21:40:56 +05:30
else
2014-03-12 13:03:51 +05:30
if ( microstructure_elemhomo ( mesh_element ( 4 , elCP ) ) . and . ip > 1_pInt ) then ! me homogenous? --> copy from first ip
2013-10-16 18:34:59 +05:30
materialpoint_P ( 1 : 3 , 1 : 3 , ip , elCP ) = materialpoint_P ( 1 : 3 , 1 : 3 , 1 , elCP )
materialpoint_F ( 1 : 3 , 1 : 3 , ip , elCP ) = materialpoint_F ( 1 : 3 , 1 : 3 , 1 , elCP )
materialpoint_dPdF ( 1 : 3 , 1 : 3 , 1 : 3 , 1 : 3 , ip , elCP ) = materialpoint_dPdF ( 1 : 3 , 1 : 3 , 1 : 3 , 1 : 3 , 1 , elCP )
2013-10-19 00:27:28 +05:30
materialpoint_results ( 1 : materialpoint_sizeResults , ip , elCP ) = &
materialpoint_results ( 1 : materialpoint_sizeResults , 1 , elCP )
2013-08-01 21:40:56 +05:30
endif
! translate from P to CS
2013-10-16 18:34:59 +05:30
Kirchhoff = math_mul33x33 ( materialpoint_P ( 1 : 3 , 1 : 3 , ip , elCP ) , math_transpose33 ( materialpoint_F ( 1 : 3 , 1 : 3 , ip , elCP ) ) )
J_inverse = 1.0_pReal / math_det33 ( materialpoint_F ( 1 : 3 , 1 : 3 , ip , elCP ) )
CPFEM_cs ( 1 : 6 , ip , elCP ) = math_Mandel33to6 ( J_inverse * Kirchhoff )
2013-08-01 21:40:56 +05:30
! translate from dP/dF to dCS/dE
H = 0.0_pReal
do i = 1 , 3 ; do j = 1 , 3 ; do k = 1 , 3 ; do l = 1 , 3 ; do m = 1 , 3 ; do n = 1 , 3
H ( i , j , k , l ) = H ( i , j , k , l ) + &
2013-10-16 18:34:59 +05:30
materialpoint_F ( j , m , ip , elCP ) * &
materialpoint_F ( l , n , ip , elCP ) * &
materialpoint_dPdF ( i , m , k , n , ip , elCP ) - &
math_I3 ( j , l ) * materialpoint_F ( i , m , ip , elCP ) * materialpoint_P ( k , m , ip , elCP ) + &
2013-08-01 21:40:56 +05:30
0.5_pReal * ( math_I3 ( i , k ) * Kirchhoff ( j , l ) + math_I3 ( j , l ) * Kirchhoff ( i , k ) + &
math_I3 ( i , l ) * Kirchhoff ( j , k ) + math_I3 ( j , k ) * Kirchhoff ( i , l ) )
enddo ; enddo ; enddo ; enddo ; enddo ; enddo
2013-10-16 18:34:59 +05:30
forall ( i = 1 : 3 , j = 1 : 3 , k = 1 : 3 , l = 1 : 3 ) &
2013-08-01 21:40:56 +05:30
H_sym ( i , j , k , l ) = 0.25_pReal * ( H ( i , j , k , l ) + H ( j , i , k , l ) + H ( i , j , l , k ) + H ( j , i , l , k ) )
2013-10-16 18:34:59 +05:30
CPFEM_dcsde ( 1 : 6 , 1 : 6 , ip , elCP ) = math_Mandel3333to66 ( J_inverse * H_sym )
2013-03-01 17:18:29 +05:30
endif
2014-08-04 23:20:01 +05:30
#endif
2013-08-01 21:40:56 +05:30
endif
2014-07-24 17:49:15 +05:30
#if defined(Marc4DAMASK) || defined(Abaqus)
2013-08-01 21:57:37 +05:30
!* remember extreme values of stress and jacobian
2013-10-16 18:34:59 +05:30
cauchyStress33 = math_Mandel6to33 ( CPFEM_cs ( 1 : 6 , ip , elCP ) )
2013-08-01 21:57:37 +05:30
if ( maxval ( cauchyStress33 ) > debug_stressMax ) then
2013-10-16 18:34:59 +05:30
debug_stressMaxLocation = [ elCP , ip ]
2013-08-01 21:57:37 +05:30
debug_stressMax = maxval ( cauchyStress33 )
endif
if ( minval ( cauchyStress33 ) < debug_stressMin ) then
2013-10-16 18:34:59 +05:30
debug_stressMinLocation = [ elCP , ip ]
2013-08-01 21:57:37 +05:30
debug_stressMin = minval ( cauchyStress33 )
endif
2013-10-16 18:34:59 +05:30
jacobian3333 = math_Mandel66to3333 ( CPFEM_dcsdE ( 1 : 6 , 1 : 6 , ip , elCP ) )
2013-08-01 21:57:37 +05:30
if ( maxval ( jacobian3333 ) > debug_jacobianMax ) then
2013-10-16 18:34:59 +05:30
debug_jacobianMaxLocation = [ elCP , ip ]
2013-08-01 21:57:37 +05:30
debug_jacobianMax = maxval ( jacobian3333 )
endif
if ( minval ( jacobian3333 ) < debug_jacobianMin ) then
2013-10-16 18:34:59 +05:30
debug_jacobianMinLocation = [ elCP , ip ]
2013-08-01 21:57:37 +05:30
debug_jacobianMin = minval ( jacobian3333 )
endif
2013-08-01 21:40:56 +05:30
!* report stress and stiffness
if ( ( iand ( debug_level ( debug_CPFEM ) , debug_levelExtensive ) / = 0_pInt ) &
2013-10-16 18:34:59 +05:30
. and . ( ( debug_e == elCP . and . debug_i == ip ) &
2013-08-01 21:40:56 +05:30
. or . . not . iand ( debug_level ( debug_CPFEM ) , debug_levelSelective ) / = 0_pInt ) ) then
2013-10-16 18:34:59 +05:30
write ( 6 , '(a,i8,1x,i2,/,12x,6(f10.3,1x)/)' ) '<< CPFEM >> stress/MPa at elFE ip ' , &
elCP , ip , CPFEM_cs ( 1 : 6 , ip , elCP ) / 1.0e6_pReal
write ( 6 , '(a,i8,1x,i2,/,6(12x,6(f10.3,1x)/))' ) '<< CPFEM >> Jacobian/GPa at elFE ip ' , &
elCP , ip , transpose ( CPFEM_dcsdE ( 1 : 6 , 1 : 6 , ip , elCP ) ) / 1.0e9_pReal
2013-08-01 21:40:56 +05:30
flush ( 6 )
endif
2014-07-24 17:49:15 +05:30
#endif
2013-03-01 17:18:29 +05:30
endif
2014-07-24 17:49:15 +05:30
#if defined(Marc4DAMASK) || defined(Abaqus)
2013-08-01 21:40:56 +05:30
!*** warn if stiffness close to zero
2013-10-16 18:34:59 +05:30
if ( all ( abs ( CPFEM_dcsdE ( 1 : 6 , 1 : 6 , ip , elCP ) ) < 1e-10_pReal ) ) call IO_warning ( 601 , elCP , ip )
2014-08-04 23:20:01 +05:30
!*** copy to output if using commercial FEM solver
2014-07-24 17:49:15 +05:30
cauchyStress = CPFEM_cs ( 1 : 6 , ip , elCP )
jacobian = CPFEM_dcsdE ( 1 : 6 , 1 : 6 , ip , elCP )
#endif
2014-07-24 14:08:52 +05:30
2012-03-09 01:55:28 +05:30
end subroutine CPFEM_general
2009-03-04 19:31:36 +05:30
2013-03-01 17:18:29 +05:30
end module CPFEM