223031012
/
DAMASK_EICMD
mirror of https://github.com/achalhp/DAMASK_EICMD.git
1
0
Fork 0
Code Issues Projects Releases Wiki Activity

polishing

This commit is contained in:
Martin Diehl 2020-06-16 06:41:53 +02:00
parent 0a9902818c
commit 54aa5a67ff
1 changed files with 13 additions and 25 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

38
src/CPFEM.f90
View File

@ -48,15 +48,12 @@ module CPFEM
logical, public, protected :: & logical, public, protected :: &
CPFEM_init_done = .false. !< remember whether init has been done already 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 :: & integer(pInt), parameter, public :: &
CPFEM_COLLECT = 2_pInt**0_pInt, & CPFEM_CALCRESULTS = 2_pInt**0_pInt, &
CPFEM_CALCRESULTS = 2_pInt**1_pInt, & CPFEM_AGERESULTS = 2_pInt**1_pInt, &
CPFEM_AGERESULTS = 2_pInt**2_pInt, & CPFEM_BACKUPJACOBIAN = 2_pInt**2_pInt, &
CPFEM_BACKUPJACOBIAN = 2_pInt**3_pInt, & CPFEM_RESTOREJACOBIAN = 2_pInt**3_pInt
CPFEM_RESTOREJACOBIAN = 2_pInt**4_pInt
public :: & public :: &
CPFEM_general, & CPFEM_general, &
@ -134,7 +131,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
real(pReal) J_inverse, & ! inverse of Jacobian real(pReal) J_inverse, & ! inverse of Jacobian
rnd rnd
real(pReal), dimension (3,3) :: Kirchhoff ! Piola-Kirchhoff stress real(pReal), dimension (3,3) :: Kirchhoff ! Piola-Kirchhoff stress
real(pReal), dimension (3,3,3,3) :: H_sym, & real(pReal), dimension (3,3,3,3) :: H_sym, &
H H
@ -142,8 +139,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
i, j, k, l, m, n, ph, homog, mySource i, j, k, l, m, n, ph, homog, mySource
logical updateJaco ! flag indicating if Jacobian has to be updated logical updateJaco ! flag indicating if Jacobian has to be updated
real(pReal), parameter :: ODD_STRESS = 1e15_pReal, & !< return value for stress in case of ping pong dummy cycle real(pReal), parameter :: ODD_STRESS = 1e15_pReal, & !< return value for stress if terminallyIll
ODD_JACOBIAN = 1e50_pReal !< return value for jacobian in case of ping pong dummy cycle ODD_JACOBIAN = 1e50_pReal !< return value for jacobian if terminallyIll
elCP = mesh_FEM2DAMASK_elem(elFE) elCP = mesh_FEM2DAMASK_elem(elFE)
@ -167,10 +164,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
if (iand(mode, CPFEM_RESTOREJACOBIAN) /= 0_pInt) & if (iand(mode, CPFEM_RESTOREJACOBIAN) /= 0_pInt) &
CPFEM_dcsde = CPFEM_dcsde_knownGood CPFEM_dcsde = CPFEM_dcsde_knownGood
!*** age results
if (iand(mode, CPFEM_AGERESULTS) /= 0_pInt) call CPFEM_forward if (iand(mode, CPFEM_AGERESULTS) /= 0_pInt) call CPFEM_forward
chosenThermal1: select case (thermal_type(material_homogenizationAt(elCP))) chosenThermal1: select case (thermal_type(material_homogenizationAt(elCP)))
case (THERMAL_conduction_ID) chosenThermal1 case (THERMAL_conduction_ID) chosenThermal1
temperature(material_homogenizationAt(elCP))%p(thermalMapping(material_homogenizationAt(elCP))%p(ip,elCP)) = & temperature(material_homogenizationAt(elCP))%p(thermalMapping(material_homogenizationAt(elCP))%p(ip,elCP)) = &
@ -179,26 +174,22 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
materialpoint_F0(1:3,1:3,ip,elCP) = ffn materialpoint_F0(1:3,1:3,ip,elCP) = ffn
materialpoint_F(1:3,1:3,ip,elCP) = ffn1 materialpoint_F(1:3,1:3,ip,elCP) = ffn1
!*** calculation of stress and jacobian
if (iand(mode, CPFEM_CALCRESULTS) /= 0_pInt) then if (iand(mode, CPFEM_CALCRESULTS) /= 0_pInt) then
!*** deformation gradient outdated or any actual deformation gradient differs more than relevantStrain from the stored one
validCalculation: if (terminallyIll) then validCalculation: if (terminallyIll) then
call random_number(rnd) call random_number(rnd)
if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal if (rnd < 0.5_pReal) rnd = rnd - 1.0_pReal
CPFEM_cs(1:6,ip,elCP) = ODD_STRESS * rnd CPFEM_cs(1:6,ip,elCP) = ODD_STRESS * rnd
CPFEM_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_identity2nd(6) CPFEM_dcsde(1:6,1:6,ip,elCP) = ODD_JACOBIAN * math_identity2nd(6)
!*** deformation gradient is not outdated
else validCalculation else validCalculation
updateJaco = mod(cycleCounter,iJacoStiffness) == 0 updateJaco = mod(cycleCounter,iJacoStiffness) == 0
FEsolving_execElem = elCP FEsolving_execElem = elCP
FEsolving_execIP = ip FEsolving_execIP = ip
if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) & if (iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) &
write(6,'(a,i8,1x,i2)') '<< CPFEM >> calculation for elFE ip ',elFE,ip write(6,'(a,i8,1x,i2)') '<< CPFEM >> calculation for elFE ip ',elFE,ip
call materialpoint_stressAndItsTangent(updateJaco, dt) call materialpoint_stressAndItsTangent(updateJaco, dt)
!* map stress and stiffness (or return odd values if terminally ill)
terminalIllness: if (terminallyIll) then terminalIllness: if (terminallyIll) then
call random_number(rnd) call random_number(rnd)
@ -208,7 +199,7 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
else terminalIllness else terminalIllness
! translate from P to CS ! translate from P to sigma
Kirchhoff = matmul(materialpoint_P(1:3,1:3,ip,elCP), transpose(materialpoint_F(1:3,1:3,ip,elCP))) Kirchhoff = matmul(materialpoint_P(1:3,1:3,ip,elCP), transpose(materialpoint_F(1:3,1:3,ip,elCP)))
J_inverse = 1.0_pReal / math_det33(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_sym33to6(J_inverse * Kirchhoff,weighted=.false.) CPFEM_cs(1:6,ip,elCP) = math_sym33to6(J_inverse * Kirchhoff,weighted=.false.)
@ -232,7 +223,6 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
endif terminalIllness endif terminalIllness
endif validCalculation endif validCalculation
!* report stress and stiffness
if ((iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) & if ((iand(debug_level(debug_CPFEM), debug_levelExtensive) /= 0_pInt) &
.and. ((debug_e == elCP .and. debug_i == ip) & .and. ((debug_e == elCP .and. debug_i == ip) &
.or. .not. iand(debug_level(debug_CPFEM), debug_levelSelective) /= 0_pInt)) then .or. .not. iand(debug_level(debug_CPFEM), debug_levelSelective) /= 0_pInt)) then
@ -245,10 +235,8 @@ subroutine CPFEM_general(mode, ffn, ffn1, temperature_inp, dt, elFE, ip, cauchyS
endif endif
!*** warn if stiffness close to zero
if (all(abs(CPFEM_dcsdE(1:6,1:6,ip,elCP)) < 1e-10_pReal)) call IO_warning(601,elCP,ip) if (all(abs(CPFEM_dcsdE(1:6,1:6,ip,elCP)) < 1e-10_pReal)) call IO_warning(601,elCP,ip)
!*** copy to output if using commercial FEM solver
cauchyStress = CPFEM_cs (1:6, ip,elCP) cauchyStress = CPFEM_cs (1:6, ip,elCP)
jacobian = CPFEM_dcsdE(1:6,1:6,ip,elCP) jacobian = CPFEM_dcsdE(1:6,1:6,ip,elCP)