2013-03-06 20:11:15 +05:30
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
|
|
|
!> @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
|
2019-09-20 01:37:18 +05:30
|
|
|
use prec
|
|
|
|
|
use numerics
|
|
|
|
|
use debug
|
|
|
|
|
use FEsolving
|
|
|
|
|
use math
|
2019-09-22 19:23:03 +05:30
|
|
|
use rotations
|
2019-09-20 01:37:18 +05:30
|
|
|
use mesh
|
|
|
|
|
use material
|
|
|
|
|
use config
|
|
|
|
|
use crystallite
|
|
|
|
|
use homogenization
|
|
|
|
|
use IO
|
|
|
|
|
use discretization
|
|
|
|
|
use DAMASK_interface
|
|
|
|
|
use numerics
|
|
|
|
|
use HDF5_utilities
|
|
|
|
|
use results
|
|
|
|
|
use lattice
|
|
|
|
|
use constitutive
|
|
|
|
|
|
|
|
|
|
implicit none
|
|
|
|
|
private
|
|
|
|
|
|
|
|
|
|
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
|
|
|
|
|
integer(pInt), public :: &
|
|
|
|
|
cycleCounter = 0_pInt, & !< needs description
|
|
|
|
|
theInc = -1_pInt, & !< needs description
|
|
|
|
|
lastLovl = 0_pInt, & !< lovl in previous call to marc hypela2
|
|
|
|
|
lastStep = 0_pInt !< kstep in previous call to abaqus umat
|
|
|
|
|
real(pReal), public :: &
|
|
|
|
|
theTime = 0.0_pReal, & !< needs description
|
|
|
|
|
theDelta = 0.0_pReal
|
|
|
|
|
logical, public :: &
|
|
|
|
|
outdatedFFN1 = .false., & !< needs description
|
|
|
|
|
lastIncConverged = .false., & !< needs description
|
|
|
|
|
outdatedByNewInc = .false. !< needs description
|
|
|
|
|
|
|
|
|
|
logical, public, protected :: &
|
|
|
|
|
CPFEM_init_done = .false. !< remember whether init has been done already
|
|
|
|
|
logical, private :: &
|
|
|
|
|
CPFEM_calc_done = .false. !< remember whether first ip has already calced the results
|
|
|
|
|
|
|
|
|
|
integer(pInt), parameter, public :: &
|
|
|
|
|
CPFEM_COLLECT = 2_pInt**0_pInt, &
|
|
|
|
|
CPFEM_CALCRESULTS = 2_pInt**1_pInt, &
|
|
|
|
|
CPFEM_AGERESULTS = 2_pInt**2_pInt, &
|
|
|
|
|
CPFEM_BACKUPJACOBIAN = 2_pInt**3_pInt, &
|
|
|
|
|
CPFEM_RESTOREJACOBIAN = 2_pInt**4_pInt
|
|
|
|
|
|
|
|
|
|
public :: &
|
|
|
|
|
CPFEM_general, &
|
|
|
|
|
CPFEM_initAll, &
|
|
|
|
|
CPFEM_results
|
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
|
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
2015-07-24 20:27:29 +05:30
|
|
|
subroutine CPFEM_initAll(el,ip)
|
2015-03-15 21:00:14 +05:30
|
|
|
integer(pInt), intent(in) :: el, & !< FE el number
|
|
|
|
|
ip !< FE integration point number
|
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
|
|
|
call DAMASK_interface_init ! Spectral and FEM interface to commandline
|
2013-03-01 17:18:29 +05:30
|
|
|
call prec_init
|
|
|
|
|
call IO_init
|
|
|
|
|
call numerics_init
|
|
|
|
|
call debug_init
|
2018-06-14 10:09:49 +05:30
|
|
|
call config_init
|
2013-03-01 17:18:29 +05:30
|
|
|
call math_init
|
2019-09-22 19:23:03 +05:30
|
|
|
call rotations_init
|
2013-03-01 17:18:29 +05:30
|
|
|
call FE_init
|
2019-04-05 20:23:41 +05:30
|
|
|
#ifdef DAMASK_HDF5
|
|
|
|
|
call HDF5_utilities_init
|
|
|
|
|
call results_init
|
|
|
|
|
#endif
|
2019-06-04 23:36:08 +05:30
|
|
|
call mesh_init(ip, el)
|
|
|
|
|
call lattice_init
|
2015-05-28 22:32:23 +05:30
|
|
|
call material_init
|
|
|
|
|
call constitutive_init
|
2014-10-10 17:58:57 +05:30
|
|
|
call crystallite_init
|
2015-07-24 20:27:29 +05:30
|
|
|
call homogenization_init
|
2013-03-01 17:18:29 +05:30
|
|
|
call CPFEM_init
|
|
|
|
|
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
|
|
|
|
2019-09-20 01:37:18 +05:30
|
|
|
write(6,'(/,a)') ' <<<+- CPFEM init -+>>>'
|
|
|
|
|
flush(6)
|
|
|
|
|
|
|
|
|
|
allocate(CPFEM_cs( 6,discretization_nIP,discretization_nElem), source= 0.0_pReal)
|
|
|
|
|
allocate(CPFEM_dcsdE( 6,6,discretization_nIP,discretization_nElem), source= 0.0_pReal)
|
|
|
|
|
allocate(CPFEM_dcsdE_knownGood(6,6,discretization_nIP,discretization_nElem), source= 0.0_pReal)
|
|
|
|
|
|
|
|
|
|
if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0) then
|
|
|
|
|
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)
|
|
|
|
|
write(6,'(a32,l1)') 'symmetricSolver: ', symmetricSolver
|
|
|
|
|
flush(6)
|
|
|
|
|
endif
|
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
|
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
2015-05-28 22:32:23 +05:30
|
|
|
subroutine CPFEM_general(mode, parallelExecution, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyStress, jacobian)
|
2013-10-16 18:34:59 +05:30
|
|
|
|
|
|
|
|
integer(pInt), intent(in) :: elFE, & !< FE element number
|
|
|
|
|
ip !< integration point number
|
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
|
2015-07-24 20:27:29 +05:30
|
|
|
real(pReal), intent(in) :: temperature_inp !< temperature
|
2013-08-02 16:50:11 +05:30
|
|
|
logical, intent(in) :: parallelExecution !< flag indicating parallel computation of requested IPs
|
2019-01-14 15:09:38 +05:30
|
|
|
real(pReal), dimension(6), intent(out) :: cauchyStress !< stress as 6 vector
|
|
|
|
|
real(pReal), dimension(6,6), intent(out) :: jacobian !< jacobian as 66 tensor (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
|
|
|
|
|
|
|
|
|
|
integer(pInt) elCP, & ! crystal plasticity element number
|
2015-05-28 22:32:23 +05:30
|
|
|
i, j, k, l, m, n, ph, homog, mySource
|
2014-07-24 17:49:15 +05:30
|
|
|
logical updateJaco ! flag indicating if JAcobian has to be updated
|
2014-08-26 19:51:24 +05:30
|
|
|
|
2013-10-16 18:34:59 +05:30
|
|
|
elCP = mesh_FEasCP('elem',elFE)
|
2014-08-26 19:51:24 +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, '#'
|
2014-08-27 21:24:50 +05:30
|
|
|
if (terminallyIll) &
|
|
|
|
|
write(6,'(a,/)') '# --- terminallyIll --- #'
|
2014-05-16 19:31:27 +05:30
|
|
|
write(6,'(a,/)') '#############################################'; flush (6)
|
2013-03-01 17:18:29 +05:30
|
|
|
endif
|
|
|
|
|
|
2013-08-02 17:06:51 +05:30
|
|
|
if (iand(mode, CPFEM_BACKUPJACOBIAN) /= 0_pInt) &
|
|
|
|
|
CPFEM_dcsde_knownGood = CPFEM_dcsde
|
|
|
|
|
if (iand(mode, CPFEM_RESTOREJACOBIAN) /= 0_pInt) &
|
|
|
|
|
CPFEM_dcsde = CPFEM_dcsde_knownGood
|
|
|
|
|
|
2019-10-24 00:47:56 +05:30
|
|
|
!*** age results
|
2013-04-16 22:37:27 +05:30
|
|
|
if (iand(mode, CPFEM_AGERESULTS) /= 0_pInt) then
|
2019-10-24 00:47:56 +05:30
|
|
|
crystallite_F0 = crystallite_partionedF ! crystallite deformation
|
2019-03-09 05:03:11 +05:30
|
|
|
crystallite_Fp0 = crystallite_Fp ! crystallite plastic deformation
|
|
|
|
|
crystallite_Lp0 = crystallite_Lp ! crystallite plastic velocity
|
|
|
|
|
crystallite_Fi0 = crystallite_Fi ! crystallite intermediate deformation
|
|
|
|
|
crystallite_Li0 = crystallite_Li ! crystallite intermediate velocity
|
|
|
|
|
crystallite_S0 = crystallite_S ! crystallite 2nd Piola Kirchhoff stress
|
2014-05-12 18:30:37 +05:30
|
|
|
|
2019-10-24 00:47:56 +05:30
|
|
|
forall (i = 1:size(plasticState)) plasticState(i)%state0 = plasticState(i)%state
|
2015-05-28 22:32:23 +05:30
|
|
|
do i = 1, size(sourceState)
|
|
|
|
|
do mySource = 1,phase_Nsources(i)
|
2019-10-24 00:47:56 +05:30
|
|
|
sourceState(i)%p(mySource)%state0 = sourceState(i)%p(mySource)%state
|
2015-05-28 22:32:23 +05:30
|
|
|
enddo; enddo
|
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'
|
2019-06-07 11:19:45 +05:30
|
|
|
if (debug_e <= discretization_nElem .and. debug_i <=discretization_nIP) then
|
2014-05-16 19:31:27 +05:30
|
|
|
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, &
|
2019-06-14 12:32:28 +05:30
|
|
|
plasticState(material_phaseAt(1,debug_e))%state(:,material_phasememberAt(1,debug_i,debug_e))
|
2014-05-12 18:30:37 +05:30
|
|
|
endif
|
2014-08-26 19:51:24 +05:30
|
|
|
endif
|
2014-07-02 17:57:39 +05:30
|
|
|
|
2014-09-19 23:29:06 +05:30
|
|
|
do homog = 1_pInt, material_Nhomogenization
|
2015-05-28 22:32:23 +05:30
|
|
|
homogState (homog)%state0 = homogState (homog)%state
|
|
|
|
|
thermalState (homog)%state0 = thermalState (homog)%state
|
|
|
|
|
damageState (homog)%state0 = damageState (homog)%state
|
2014-09-19 23:29:06 +05:30
|
|
|
enddo
|
2019-03-09 03:46:08 +05:30
|
|
|
endif
|
2013-08-01 21:40:56 +05:30
|
|
|
|
2013-08-02 17:06:51 +05:30
|
|
|
|
|
|
|
|
|
|
|
|
|
!*** collection of FEM input with returning of randomize odd stress and jacobian
|
2014-08-27 21:24:50 +05:30
|
|
|
!* If no parallel execution is required, there is no need to collect FEM input
|
2014-08-26 19:51:24 +05:30
|
|
|
|
2013-08-02 17:06:51 +05:30
|
|
|
if (.not. parallelExecution) then
|
2019-06-07 09:48:42 +05:30
|
|
|
chosenThermal1: select case (thermal_type(material_homogenizationAt(elCP)))
|
2017-04-28 16:09:01 +05:30
|
|
|
case (THERMAL_conduction_ID) chosenThermal1
|
2018-10-04 10:09:03 +05:30
|
|
|
temperature(material_homogenizationAt(elCP))%p(thermalMapping(material_homogenizationAt(elCP))%p(ip,elCP)) = &
|
2017-04-24 12:40:37 +05:30
|
|
|
temperature_inp
|
2017-04-28 16:09:01 +05:30
|
|
|
end select chosenThermal1
|
2013-10-16 18:34:59 +05:30
|
|
|
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
|
|
|
|
|
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)
|
2019-06-07 09:48:42 +05:30
|
|
|
chosenThermal2: select case (thermal_type(material_homogenizationAt(elCP)))
|
2017-04-28 16:09:01 +05:30
|
|
|
case (THERMAL_conduction_ID) chosenThermal2
|
2018-10-04 10:09:03 +05:30
|
|
|
temperature(material_homogenizationAt(elCP))%p(thermalMapping(material_homogenizationAt(elCP))%p(ip,elCP)) = &
|
2017-04-24 12:40:37 +05:30
|
|
|
temperature_inp
|
2017-04-28 16:09:01 +05:30
|
|
|
end select chosenThermal2
|
2013-10-16 18:34:59 +05:30
|
|
|
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.
|
2014-08-27 21:24:50 +05:30
|
|
|
endif ! collection
|
2013-08-01 21:40:56 +05:30
|
|
|
|
2014-08-26 19:51:24 +05:30
|
|
|
|
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
|
2014-08-26 19:51:24 +05:30
|
|
|
|
2013-08-01 21:40:56 +05:30
|
|
|
!*** deformation gradient outdated or any actual deformation gradient differs more than relevantStrain from the stored one
|
2014-08-27 21:24:50 +05:30
|
|
|
validCalculation: if (terminallyIll &
|
|
|
|
|
.or. outdatedFFN1 &
|
|
|
|
|
.or. any(abs(ffn1 - materialpoint_F(1:3,1:3,ip,elCP)) > defgradTolerance)) then
|
2014-08-26 19:51:24 +05:30
|
|
|
if (any(abs(ffn1 - materialpoint_F(1:3,1:3,ip,elCP)) > defgradTolerance)) then
|
|
|
|
|
if (iand(debug_level(debug_CPFEM), debug_levelBasic) /= 0_pInt) then
|
2014-08-27 21:24:50 +05:30
|
|
|
write(6,'(a,1x,i8,1x,i2)') '<< CPFEM >> OUTDATED at elFE ip',elFE,ip
|
2013-10-19 00:27:28 +05:30
|
|
|
write(6,'(a,/,3(12x,3(f10.6,1x),/))') '<< CPFEM >> FFN1 old:',&
|
2019-01-14 15:09:38 +05:30
|
|
|
transpose(materialpoint_F(1:3,1:3,ip,elCP))
|
|
|
|
|
write(6,'(a,/,3(12x,3(f10.6,1x),/))') '<< CPFEM >> FFN1 now:',transpose(ffn1)
|
2013-08-01 21:40:56 +05:30
|
|
|
endif
|
|
|
|
|
outdatedFFN1 = .true.
|
|
|
|
|
endif
|
|
|
|
|
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-08-26 19:51:24 +05:30
|
|
|
|
2013-08-01 21:40:56 +05:30
|
|
|
!*** deformation gradient is not outdated
|
2014-08-26 19:51:24 +05:30
|
|
|
|
2014-08-27 21:24:50 +05:30
|
|
|
else validCalculation
|
2014-08-26 19:51:24 +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
|
2014-08-26 19:51:24 +05:30
|
|
|
|
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
|
2019-06-05 00:46:18 +05:30
|
|
|
FEsolving_execIP(1,elCP) = ip
|
|
|
|
|
FEsolving_execIP(2,elCP) = ip
|
|
|
|
|
if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) &
|
|
|
|
|
write(6,'(a,i8,1x,i2)') '<< CPFEM >> calculation for elFE ip ',elFE,ip
|
|
|
|
|
call materialpoint_stressAndItsTangent(updateJaco, dt) ! calculate stress and its tangent
|
2014-08-26 19:51:24 +05:30
|
|
|
|
2013-08-01 21:40:56 +05:30
|
|
|
!* parallel computation and calulation not yet done
|
2014-08-26 19:51:24 +05:30
|
|
|
|
2013-08-01 21:40:56 +05:30
|
|
|
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)
|
2013-08-01 21:40:56 +05:30
|
|
|
CPFEM_calc_done = .true.
|
2013-03-01 17:18:29 +05:30
|
|
|
endif
|
2014-08-26 19:51:24 +05:30
|
|
|
|
2013-08-01 21:40:56 +05:30
|
|
|
!* map stress and stiffness (or return odd values if terminally ill)
|
2014-08-27 21:24:50 +05:30
|
|
|
terminalIllness: 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)
|
2014-08-27 21:24:50 +05:30
|
|
|
|
|
|
|
|
else terminalIllness
|
|
|
|
|
|
|
|
|
|
|
2013-08-01 21:40:56 +05:30
|
|
|
! translate from P to CS
|
2019-04-03 11:52:04 +05:30
|
|
|
Kirchhoff = matmul(materialpoint_P(1:3,1:3,ip,elCP), transpose(materialpoint_F(1:3,1:3,ip,elCP)))
|
2013-10-16 18:34:59 +05:30
|
|
|
J_inverse = 1.0_pReal / math_det33(materialpoint_F(1:3,1:3,ip,elCP))
|
2019-01-14 21:06:08 +05:30
|
|
|
CPFEM_cs(1:6,ip,elCP) = math_sym33to6(J_inverse * Kirchhoff,weighted=.false.)
|
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
|
2019-01-14 15:09:38 +05:30
|
|
|
H(i,j,k,l) = H(i,j,k,l) &
|
|
|
|
|
+ materialpoint_F(j,m,ip,elCP) * materialpoint_F(l,n,ip,elCP) &
|
|
|
|
|
* materialpoint_dPdF(i,m,k,n,ip,elCP) &
|
|
|
|
|
- math_delta(j,l) * materialpoint_F(i,m,ip,elCP) * materialpoint_P(k,m,ip,elCP) &
|
|
|
|
|
+ 0.5_pReal * ( Kirchhoff(j,l)*math_delta(i,k) + Kirchhoff(i,k)*math_delta(j,l) &
|
|
|
|
|
+ Kirchhoff(j,k)*math_delta(i,l) + Kirchhoff(i,l)*math_delta(j,k))
|
2013-08-01 21:40:56 +05:30
|
|
|
enddo; enddo; enddo; enddo; enddo; enddo
|
2014-08-27 21:24:50 +05:30
|
|
|
|
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))
|
2014-08-27 21:24:50 +05:30
|
|
|
|
2019-01-14 21:06:08 +05:30
|
|
|
CPFEM_dcsde(1:6,1:6,ip,elCP) = math_sym3333to66(J_inverse * H_sym,weighted=.false.)
|
2013-08-01 21:40:56 +05:30
|
|
|
|
2014-08-27 21:24:50 +05:30
|
|
|
endif terminalIllness
|
|
|
|
|
endif validCalculation
|
2013-08-01 21:57:37 +05:30
|
|
|
|
2014-08-26 19:51:24 +05:30
|
|
|
!* report stress and stiffness
|
2013-08-01 21:40:56 +05:30
|
|
|
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
|
2014-08-27 21:24:50 +05:30
|
|
|
write(6,'(a,i8,1x,i2,/,12x,6(f10.3,1x)/)') &
|
|
|
|
|
'<< CPFEM >> stress/MPa at elFE ip ', elFE, ip, CPFEM_cs(1:6,ip,elCP)*1.0e-6_pReal
|
|
|
|
|
write(6,'(a,i8,1x,i2,/,6(12x,6(f10.3,1x)/))') &
|
|
|
|
|
'<< CPFEM >> Jacobian/GPa at elFE ip ', elFE, ip, transpose(CPFEM_dcsdE(1:6,1:6,ip,elCP))*1.0e-9_pReal
|
2013-08-01 21:40:56 +05:30
|
|
|
flush(6)
|
|
|
|
|
endif
|
2014-08-27 21:24:50 +05:30
|
|
|
|
2013-03-01 17:18:29 +05:30
|
|
|
endif
|
|
|
|
|
|
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-27 21:24:50 +05:30
|
|
|
|
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)
|
2014-08-27 21:24:50 +05:30
|
|
|
|
|
|
|
|
|
|
|
|
|
!*** remember extreme values of stress ...
|
2019-01-14 21:06:08 +05:30
|
|
|
cauchyStress33 = math_6toSym33(CPFEM_cs(1:6,ip,elCP),weighted=.false.)
|
2014-08-27 21:24:50 +05:30
|
|
|
if (maxval(cauchyStress33) > debug_stressMax) then
|
|
|
|
|
debug_stressMaxLocation = [elCP, ip]
|
|
|
|
|
debug_stressMax = maxval(cauchyStress33)
|
|
|
|
|
endif
|
|
|
|
|
if (minval(cauchyStress33) < debug_stressMin) then
|
|
|
|
|
debug_stressMinLocation = [elCP, ip]
|
|
|
|
|
debug_stressMin = minval(cauchyStress33)
|
|
|
|
|
endif
|
|
|
|
|
!*** ... and Jacobian
|
2019-01-14 21:06:08 +05:30
|
|
|
jacobian3333 = math_66toSym3333(CPFEM_dcsdE(1:6,1:6,ip,elCP),weighted=.false.)
|
2014-08-27 21:24:50 +05:30
|
|
|
if (maxval(jacobian3333) > debug_jacobianMax) then
|
|
|
|
|
debug_jacobianMaxLocation = [elCP, ip]
|
|
|
|
|
debug_jacobianMax = maxval(jacobian3333)
|
|
|
|
|
endif
|
|
|
|
|
if (minval(jacobian3333) < debug_jacobianMin) then
|
|
|
|
|
debug_jacobianMinLocation = [elCP, ip]
|
|
|
|
|
debug_jacobianMin = minval(jacobian3333)
|
|
|
|
|
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
|
|
|
|
2019-05-05 15:36:55 +05:30
|
|
|
|
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
|
|
|
!> @brief triggers writing of the results
|
|
|
|
|
!--------------------------------------------------------------------------------------------------
|
|
|
|
|
subroutine CPFEM_results(inc,time)
|
|
|
|
|
|
2019-09-20 01:37:18 +05:30
|
|
|
integer(pInt), intent(in) :: inc
|
|
|
|
|
real(pReal), intent(in) :: time
|
2019-05-05 15:36:55 +05:30
|
|
|
|
|
|
|
|
#ifdef DAMASK_HDF5
|
2019-09-20 01:37:18 +05:30
|
|
|
call results_openJobFile
|
|
|
|
|
call results_addIncrement(inc,time)
|
|
|
|
|
call constitutive_results
|
|
|
|
|
call crystallite_results
|
2019-10-12 19:20:10 +05:30
|
|
|
call homogenization_results
|
2019-09-20 01:37:18 +05:30
|
|
|
call results_removeLink('current') ! ToDo: put this into closeJobFile
|
|
|
|
|
call results_closeJobFile
|
2019-05-05 15:36:55 +05:30
|
|
|
#endif
|
|
|
|
|
|
|
|
|
|
end subroutine CPFEM_results
|
|
|
|
|
|
2013-03-01 17:18:29 +05:30
|
|
|
end module CPFEM
|
