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

# separated CPFEM_Crystallite from CPFEM. Now in separate module crystallite.f90 as "SingleCrystallite" 

# improved SingleCrystallite to advance by true cutbacking (instead of improving guess and integrating always from t_0)
# module "crystal" renamed to "lattice" together with its prefix for variables
# extension of "computationMode" to deal with cutbacks (CPFEM_general).
# cutback and new inc detection for MARC is based on common block variable cptim (and inc), not incsub anymore!
# generalized GrainInterAction as new homogenization scheme

# two symbolic links are required: constitutive.f90 and CPFEM.f90
This commit is contained in:
Denny Tjahjanto 2008-04-07 14:54:29 +00:00
parent c063ce5bc1
commit 12dfbaf6b4
11 changed files with 1887 additions and 61 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

798
trunk/CPFEM_GIA8.f90 Normal file
View File

@ -0,0 +1,798 @@
!##############################################################
MODULE CPFEM
!##############################################################
! *** CPFEM engine ***
!
use prec, only: pReal,pInt
implicit none
!
! ****************************************************************
! *** General variables for the material behaviour calculation ***
! ****************************************************************
real(pReal), dimension (:,:), allocatable :: CPFEM_Temperature
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_ffn_bar
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_ffn1_bar
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_PK1_bar
real(pReal), dimension (:,:,:,:,:,:),allocatable :: CPFEM_dPdF_bar
real(pReal), dimension (:,:,:), allocatable :: CPFEM_stress_bar
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_jaco_bar
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_jaco_knownGood
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_results
real(pReal), dimension (:,:,:,:,:), allocatable :: CPFEM_Fp_old
real(pReal), dimension (:,:,:,:,:), allocatable :: CPFEM_Fp_new
real(pReal), parameter :: CPFEM_odd_stress = 1e15_pReal, CPFEM_odd_jacobian = 1e50_pReal
integer(pInt) :: CPFEM_Nresults = 4_pInt ! three Euler angles plus volume fraction
logical :: CPFEM_init_done = .false. ! remember if init has been done already
logical :: CPFEM_calc_done = .false. ! remember if first IP has already calced the results
!
real(pReal), dimension (:,:,:,:), allocatable :: GIA_rVect_new ! boundary relaxation vectors
real(pReal), dimension (:,:,:,:), allocatable :: GIA_rVect_old ! boundary relaxation vectors
real(pReal), dimension (:,:), allocatable :: GIA_bNorm ! grain boundary normals
!
CONTAINS
!
!*********************************************************
!*** allocate the arrays defined in module CPFEM ***
!*** and initialize them ***
!*********************************************************
SUBROUTINE CPFEM_init(Temperature)
!
use prec
use math, only: math_EulertoR, math_I3, math_identity2nd
use mesh
use constitutive
!
implicit none
!
real(pReal) Temperature
integer(pInt) e,i,g,b
!
! *** mpie.marc parameters ***
allocate(CPFEM_Temperature (mesh_maxNips,mesh_NcpElems)) ; CPFEM_Temperature = Temperature
allocate(CPFEM_ffn_bar (3,3,mesh_maxNips,mesh_NcpElems))
forall(e=1:mesh_NcpElems,i=1:mesh_maxNips) CPFEM_ffn_bar(:,:,i,e) = math_I3
allocate(CPFEM_ffn1_bar (3,3,mesh_maxNips,mesh_NcpElems)) ; CPFEM_ffn1_bar = CPFEM_ffn_bar
allocate(CPFEM_PK1_bar (3,3,mesh_maxNips,mesh_NcpElems)) ; CPFEM_PK1_bar = 0.0_pReal
allocate(CPFEM_dPdF_bar(3,3,3,3,mesh_maxNips,mesh_NcpElems)) ; CPFEM_dPdF_bar = 0.0_pReal
allocate(CPFEM_stress_bar(6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_stress_bar = 0.0_pReal
allocate(CPFEM_jaco_bar(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_jaco_bar = 0.0_pReal
allocate(CPFEM_jaco_knownGood(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_jaco_knownGood = 0.0_pReal
!
! *** User defined results !!! MISSING incorporate consti_Nresults ***
allocate(CPFEM_results(CPFEM_Nresults+constitutive_maxNresults,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems))
CPFEM_results = 0.0_pReal
!
! *** Plastic deformation gradient at (t=t0) and (t=t1) ***
allocate(CPFEM_Fp_new(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_Fp_new = 0.0_pReal
allocate(CPFEM_Fp_old(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems))
forall (e=1:mesh_NcpElems,i=1:mesh_maxNips,g=1:constitutive_maxNgrains) &
CPFEM_Fp_old(:,:,g,i,e) = math_EulerToR(constitutive_EulerAngles(:,g,i,e)) ! plastic def gradient reflects init orientation
!
allocate(GIA_rVect_new(3,12,mesh_maxNips,mesh_NcpElems)) ; GIA_rVect_new = 0.0_pReal
allocate(GIA_rVect_old(3,12,mesh_maxNips,mesh_NcpElems)) ; GIA_rVect_old = 0.0_pReal
allocate(GIA_bNorm(3,12)) ; GIA_bNorm = 0.0_pReal
do b = 1,4
GIA_bNorm(1,b) = 1.0_pReal
GIA_bNorm(2,b+4) = 1.0_pReal
GIA_bNorm(3,b+8) = 1.0_pReal
enddo
!
! *** Output to MARC output file ***
write(6,*)
write(6,*) 'CPFEM Initialization'
write(6,*)
write(6,*) 'CPFEM_Temperature: ', shape(CPFEM_Temperature)
write(6,*) 'CPFEM_ffn_bar: ', shape(CPFEM_ffn_bar)
write(6,*) 'CPFEM_ffn1_bar: ', shape(CPFEM_ffn1_bar)
write(6,*) 'CPFEM_PK1_bar: ', shape(CPFEM_PK1_bar)
write(6,*) 'CPFEM_dPdF_bar: ', shape(CPFEM_dPdF_bar)
write(6,*) 'CPFEM_stress_bar: ', shape(CPFEM_stress_bar)
write(6,*) 'CPFEM_jaco_bar: ', shape(CPFEM_jaco_bar)
write(6,*) 'CPFEM_jaco_knownGood: ', shape(CPFEM_jaco_knownGood)
write(6,*) 'CPFEM_results: ', shape(CPFEM_results)
write(6,*) 'CPFEM_Fp_old: ', shape(CPFEM_Fp_old)
write(6,*) 'CPFEM_Fp_new: ', shape(CPFEM_Fp_new)
!
write(6,*) 'GIA_rVect_new: ', shape(GIA_rVect_new)
write(6,*) 'GIA_rVect_old: ', shape(GIA_rVect_old)
write(6,*) 'GIA_bNorm: ', shape(GIA_bNorm)
write(6,*)
call flush(6)
return
!
END SUBROUTINE
!
!
!***********************************************************************
!*** perform initialization at first call, update variables and ***
!*** call the actual material model ***
!
! CPFEM_mode computation mode (regular, collection, recycle)
! ffn deformation gradient for t=t0
! ffn1 deformation gradient for t=t1
! Temperature temperature
! CPFEM_dt time increment
! CPFEM_en element number
! CPFEM_in intergration point number
! CPFEM_stress stress vector in Mandel notation
! CPFEM_updateJaco flag to initiate computation of Jacobian
! CPFEM_jaco jacobian in Mandel notation
! CPFEM_ngens size of stress strain law
!***********************************************************************
SUBROUTINE CPFEM_general(CPFEM_mode, ffn, ffn1, Temperature, CPFEM_dt,&
CPFEM_en, CPFEM_in, CPFEM_stress, CPFEM_updateJaco, CPFEM_jaco, CPFEM_ngens)
! note: CPFEM_stress = Cauchy stress cs(6) and CPFEM_jaco = Consistent tangent dcs/de
!
use prec, only: pReal,pInt
use FEsolving
use debug
use math, only: math_init, invnrmMandel, math_identity2nd, math_Mandel3333to66,math_Mandel33to6,math_Mandel6to33,math_det3x3,math_I3
use mesh, only: mesh_init,mesh_FEasCP, mesh_NcpElems, FE_Nips, FE_mapElemtype, mesh_element
use lattice, only: lattice_init
use constitutive, only: constitutive_init,constitutive_state_old,constitutive_state_new,material_Cslip_66
implicit none
!
integer(pInt) CPFEM_en, CPFEM_in, cp_en, CPFEM_ngens, i,j,k,l,m,n, e
real(pReal), dimension (3,3) :: ffn,ffn1,Kirchhoff_bar
real(pReal), dimension (3,3,3,3) :: H_bar
real(pReal), dimension(CPFEM_ngens) :: CPFEM_stress
real(pReal), dimension(CPFEM_ngens,CPFEM_ngens) :: CPFEM_jaco
real(pReal) Temperature,CPFEM_dt,J_inverse
integer(pInt) CPFEM_mode ! 1: regular computation with aged results&
! 2: regular computation&
! 3: collection of FEM data&
! 4: recycling of former results (MARC speciality)&
! 5: record tangent from former converged inc&
! 6: restore tangent from former converged inc
logical CPFEM_updateJaco
!
if (.not. CPFEM_init_done) then ! initialization step (three dimensional stress state check missing?)
call math_init()
call mesh_init()
call lattice_init()
call constitutive_init()
call CPFEM_init(Temperature)
CPFEM_init_done = .true.
endif
!
cp_en = mesh_FEasCP('elem',CPFEM_en)
if (cp_en == 1 .and. CPFEM_in == 1) &
write(6,'(a6,x,i4,x,a4,x,i4,x,a10,x,f8.4,x,a10,x,i2,x,a10,x,i2,x,a10,x,i2,x,a10,x,i2)') &
'elem',cp_en,'IP',CPFEM_in,&
'theTime',theTime,'theInc',theInc,'theCycle',theCycle,'theLovl',theLovl,&
'mode',CPFEM_mode
!
select case (CPFEM_mode)
case (2,1) ! regular computation (with aging of results)
if (.not. CPFEM_calc_done) then ! puuh, me needs doing all the work...
write (6,*) 'puuh me needs doing all the work', cp_en
if (CPFEM_mode == 1) then ! age results at start of new increment
CPFEM_Fp_old = CPFEM_Fp_new
constitutive_state_old = constitutive_state_new
GIA_rVect_old = GIA_rVect_new
write (6,*) '#### aged results'
endif
debug_cutbackDistribution = 0_pInt ! initialize debugging data
debug_InnerLoopDistribution = 0_pInt
debug_OuterLoopDistribution = 0_pInt
!
do e=1,mesh_NcpElems ! ## this shall be done in a parallel loop in the future ##
do i=1,FE_Nips(mesh_element(2,e)) ! iterate over all IPs of this element's type
debugger = (e==1 .and. i==1) ! switch on debugging for first IP in first element
call CPFEM_MaterialPoint(CPFEM_updateJaco, CPFEM_dt, i, e)
enddo
enddo
call debug_info() ! output of debugging/performance statistics
CPFEM_calc_done = .true. ! now calc is done
endif
! translate from P and dP/dF to CS and dCS/dE
Kirchhoff_bar = matmul(CPFEM_PK1_bar(:,:,CPFEM_in, cp_en),transpose(CPFEM_ffn1_bar(:,:,CPFEM_in, cp_en)))
J_inverse = 1.0_pReal/math_det3x3(CPFEM_ffn1_bar(:,:,CPFEM_in, cp_en))
CPFEM_stress_bar(1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel33to6(J_inverse*Kirchhoff_bar)
!
H_bar = 0.0_pReal
forall(i=1:3,j=1:3,k=1:3,l=1:3,m=1:3,n=1:3) &
H_bar(i,j,k,l) = H_bar(i,j,k,l) + &
(CPFEM_ffn1_bar(j,m,CPFEM_in,cp_en)*CPFEM_ffn1_bar(l,n,CPFEM_in,cp_en)*CPFEM_dPdF_bar(i,m,k,n,CPFEM_in,cp_en) - &
math_I3(j,l)*CPFEM_ffn1_bar(i,m,CPFEM_in,cp_en)*CPFEM_PK1_bar(k,m,CPFEM_in,cp_en)) + &
0.5_pReal*(math_I3(i,k)*Kirchhoff_bar(j,l) + math_I3(j,l)*Kirchhoff_bar(i,k) + &
math_I3(i,l)*Kirchhoff_bar(j,k) + math_I3(j,k)*Kirchhoff_bar(i,l))
CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel3333to66(J_inverse*H_bar)
!
case (3) ! collect and return odd result
CPFEM_Temperature(CPFEM_in,cp_en) = Temperature
CPFEM_ffn_bar(:,:,CPFEM_in,cp_en) = ffn
CPFEM_ffn1_bar(:,:,CPFEM_in,cp_en) = ffn1
CPFEM_stress_bar(1:CPFEM_ngens,CPFEM_in,cp_en) = CPFEM_odd_stress
CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = CPFEM_odd_jacobian*math_identity2nd(CPFEM_ngens)
CPFEM_calc_done = .false.
case (4) ! do nothing since we can recycle the former results (MARC specialty)
case (5) ! record consistent tangent at beginning of new increment
CPFEM_jaco_knownGood = CPFEM_jaco_bar
case (6) ! restore consistent tangent after cutback
CPFEM_jaco_bar = CPFEM_jaco_knownGood
end select
!
! return the local stress and the jacobian from storage
CPFEM_stress(1:CPFEM_ngens) = CPFEM_stress_bar(1:CPFEM_ngens,CPFEM_in,cp_en)
CPFEM_jaco(1:CPFEM_ngens,1:CPFEM_ngens) = CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en)
if (cp_en == 1 .and. CPFEM_in == 1) write (6,*) 'stress',CPFEM_stress
if (cp_en == 1 .and. CPFEM_in == 1 .and. CPFEM_updateJaco) write (6,*) 'stiffness',CPFEM_jaco
!
return
!
END SUBROUTINE
!
!**********************************************************
!*** calculate the material point behaviour ***
!**********************************************************
SUBROUTINE CPFEM_MaterialPoint(&
updateJaco,& ! flag to initiate Jacobian updating
CPFEM_dt,& ! Time increment (dt)
CPFEM_in,& ! Integration point number
cp_en) ! Element number
!
use prec
use debug
use math, only: math_pDecomposition,math_RtoEuler,inDeg,math_I3,math_invert3x3,math_permut,math_invert,math_delta
use IO, only: IO_error
use mesh, only: mesh_element
use crystallite
use constitutive
implicit none
!
character(len=128) msg
integer(pInt) cp_en,CPFEM_in,grain,max_cutbacks,i,j,k,l,m,n,iBoun,NRiter,dummy,ii,jj,kk,ll
logical updateJaco,error,NRconvergent,failed
real(pReal) CPFEM_dt,volfrac,dTime,shMod,C_kb,resNorm,resMax,subStep,subFrac,temp1,temp2
real(pReal), dimension(3,3) :: F0_bar,F1_bar,dF_bar
real(pReal), dimension(3,3) :: U,R
real(pReal), dimension(3,3,8) :: PK1,Fp0,Fp1,Fe1,F1,F0,dFgrain,dFg_cor
real(pReal), dimension(3,3,12) :: GPK1,GF1,Nye,GRB1
real(pReal), dimension(3,3,3,3,8) :: dPdF
real(pReal), dimension(3,3,3,3,12) :: dRdX1
real(pReal), dimension(36) :: var,res
real(pReal), dimension(36,36) :: dresdvar,dvardres
real(pReal), dimension(3,12) :: rx,rVect
real(pReal), dimension(12) :: NyeNorm
real(pReal), dimension(constitutive_maxNstatevars,8) :: state0,state1
!
if (texture_Ngrains(mesh_element(4,cp_en)) /= 8_pInt) then
call IO_error(800)
return
endif
!
CPFEM_PK1_bar(:,:,CPFEM_in,cp_en) = 0.0_pReal ! zero out average first PK stress
if (updateJaco) CPFEM_dPdF_bar(:,:,:,:,CPFEM_in,cp_en) = 0.0_pReal ! zero out average consistent tangent
!
! ------------- GIA loop --------------------
!
! collect information
shMod = 0.2_pReal*(material_C11(1) - material_C12(1)) + 0.3_pReal*material_C44(1) ! equivalent shear modulus
C_kb = material_bg(1)*shMod/material_GrainSize(1) ! equivalent boundary stiffness
!
F0_bar = CPFEM_ffn_bar(:,:,CPFEM_in,cp_en) ! effective deformation gradient at t_n
state0 = constitutive_state_old(:,:,CPFEM_in,cp_en) ! state variables at t_n
Fp0 = CPFEM_Fp_old(:,:,:,CPFEM_in,cp_en) ! grain plastic def. gradient at t_n
rVect = GIA_rVect_old(:,:,CPFEM_in,cp_en) ! relaxation vectors from previous convergent step
!
dF_bar = CPFEM_ffn1_bar(:,:,CPFEM_in,cp_en) - CPFEM_ffn_bar(:,:,CPFEM_in,cp_en) ! deformation gradient increment
subFrac = 0.0_pReal
subStep = 1.0_pReal
!
! Substepping procedure to improve N-R iteration
SubStepping: do
dTime = subStep*CPFEM_dt
call GIA_RelaxedDeformation(F0,F0_bar,rVect) ! def. gradient of indiv. grains at t_n
F1_bar = F0_bar + subStep*dF_bar ! effective def. gradient at t_n+1
forall (iBoun=1:12,i=1:3) var(3_pInt*(iBoun-1_pInt)+i) = rVect(i,iBoun) ! primary variable: relaxation vector
!
! Newton-Raphson iteration block
NRiter = 1_pInt
NRIteration: do
forall (iBoun=1:12,i=1:3) rx(i,iBoun) = var(3_pInt*(iBoun-1_pInt)+i) ! relaxation vectors (guess)
!
! deformation gradients of grains at t_n+1 (guess)
call GIA_RelaxedDeformation(F1,F1_bar,rx)
!
! -------------- grain loop -----------------
do grain = 1,8
call SingleCrystallite(msg,PK1(:,:,grain),dPdF(:,:,:,:,grain),&
CPFEM_results(5:4+constitutive_Nresults(grain,CPFEM_in,cp_en),grain,CPFEM_in,cp_en),&
Fp1(:,:,grain),Fe1(:,:,grain),state1(:,grain),& ! output up to here
dTime,cp_en,CPFEM_in,grain,.true.,&
CPFEM_Temperature(CPFEM_in,cp_en),F1(:,:,grain),F0(:,:,grain),Fp0(:,:,grain),state0(:,grain))
if (msg /= 'ok') then ! solution not reached --> exit NRIteration
write(6,*) 'GIA: grain loop failed to converge within allowable step-size'
NRconvergent = .false.
exit NRiteration
endif
enddo ! grain loop
!
! calculate the deformation jump and stress jump across the boundaries
dFgrain = F1 - F0
call GIA_BoundaryJump(GF1,F1)
call GIA_BoundaryJump(GPK1,PK1)
!
! compute the Nye tensor at the boundary
Nye = 0.0_pReal
NyeNorm = 0.0_pReal
do iBoun = 1,12
do i = 1,3
do j = 1,3
do k = 1,3
do l = 1,3
Nye(i,j,iBoun) = Nye(i,j,iBoun) - 0.5_pReal*math_permut(j,k,l)*GIA_bNorm(k,iBoun)*GF1(i,l,iBoun)
enddo
enddo
NyeNorm(iBoun) = NyeNorm(iBoun) + Nye(i,j,iBoun)*Nye(i,j,iBoun)
enddo
enddo
NyeNorm(iBoun) = sqrt(NyeNorm(iBoun))
if (NyeNorm(iBoun) > 1.0e-8_pReal) Nye(:,:,iBoun) = Nye(:,:,iBoun)/NyeNorm(iBoun)
enddo
!
! compute the stress-like penalty at the boundary
GRB1 = 0.0_pReal
do iBoun = 1,12
do i = 1,3
do j = 1,3
do k = 1,3
do l = 1,3
GRB1(i,j,iBoun) = GRB1(i,j,iBoun) + Nye(i,k,iBoun)*GIA_bNorm(l,iBoun)*math_permut(k,l,j)
enddo
enddo
enddo
enddo
GRB1(:,:,iBoun) = 0.5_pReal*(C_kb + C_kb)*GRB1(:,:,iBoun)
enddo
!
! compute the resiudal of stress at the boundary
res = 0.0_pReal
resNorm = 0.0_pReal
do iBoun = 1,12
do j = 1,3
do i = 1,3
res(3_pInt*(iBoun-1_pInt)+j) = res(3_pInt*(iBoun-1_pInt)+j) - &
GIA_bNorm(i,iBoun)*(GPK1(i,j,iBoun) - GRB1(i,j,iBoun))
enddo
resNorm = resNorm + res(3_pInt*(iBoun-1_pInt)+j)*res(3_pInt*(iBoun-1_pInt)+j)
enddo
enddo
resNorm = sqrt(resNorm)
!
! write(6,'(x,a,i3,a,i3,a,i3,a,e10.4)')'EL = ',cp_en,':IP = ',CPFEM_in,':Iter = ',NRiter,':RNorm = ',resNorm
if (NRiter == 1_pInt) resMax = resNorm
if ((resNorm < resToler*resMax) .or. (resNorm < resAbsol)) then ! resNorm < tolerance ===> convergent
NRconvergent = .true.
exit NRiteration
elseif ((NRiter > NRiterMax) .or. (resNorm > resBound*resMax)) then ! resNorm > up. bound ===> substepping
NRconvergent = .false.
exit NRiteration
else ! update the residual
dRdX1 = 0.0_pReal
do iBoun = 1,12
if (NyeNorm(iBoun) < 1.0e-8_pReal) NyeNorm(iBoun) = 1.0e-8_pReal
do i = 1,3
do j = 1,3
do k = 1,3
do l = 1,3
temp1 = 0.0_pReal
temp2 = 0.0_pReal
do ii = 1,3
do jj = 1,3
do kk = 1,3
temp1 = temp1 + GIA_bNorm(jj,iBoun)*math_permut(ii,jj,j)*math_delta(i,k)* &
GIA_bNorm(kk,iBoun)*math_permut(ii,kk,l)
do ll = 1,3
temp2 = temp2 + Nye(i,ii,iBoun)*GIA_bNorm(jj,iBoun)*math_permut(ii,jj,j)* &
Nye(k,kk,iBoun)*GIA_bNorm(ll,iBoun)*math_permut(kk,ll,l)
enddo
enddo
enddo
enddo
dRdX1(i,j,k,l,iBoun) = 0.25_pReal*(C_kb + C_kb)*(temp1 - temp2)/NyeNorm(iBoun)
enddo
enddo
enddo
enddo
enddo
call GIA_JacobianMatrix(dresdvar,dPdF,dRdX1)
dvardres = 0.0_pReal
call math_invert(36,dresdvar,dvardres,dummy,failed)
if (failed) then
write(6,*) 'GIA: failed to invert the Jacobian'
NRconvergent = .false.
exit NRiteration
endif
forall (i=1:36,j=1:36) var(i) = var(i) - dvardres(i,j)*res(j)
endif
!
NRiter = NRiter + 1_pInt
enddo NRIteration ! End of N-R iteration blok
!
if (.not. NRconvergent) then
subStep = 0.5_pReal*subStep
else
subFrac = subFrac + subStep
subStep = 1.0_pReal - subFrac
Fp0 = Fp1
F0_bar = F1_bar
state0 = state1
rVect = rx
endif
!
if (subStep < subStepMin) exit SubStepping
enddo SubStepping ! End of substepping blok
!
! ------------- GIA loop (end) --------------
!
! return to the general subroutine when convergence is not reached
if (.not. NRconvergent) then
write(6,'(x,a)') 'GIA: convergence is not reached within allowable step-size'
write(6,'(x,a,i3,a,i3)') 'Element: ',cp_en,' @ IP: ',CPFEM_in
call IO_error(600)
return
endif
!
! updates all variables, deformation gradients, and vectors
GIA_rVect_new(:,:,CPFEM_in,cp_en) = rVect
CPFEM_Fp_new(:,:,:,CPFEM_in,cp_en) = Fp1
constitutive_state_new(:,:,CPFEM_in,cp_en) = state1
!
! compute the effective stress and consistent tangent
call GIA_TangentCorrection(dFgrain,dFg_cor)
do grain = 1,texture_Ngrains(mesh_element(4,cp_en))
volfrac = constitutive_matVolFrac(grain,CPFEM_in,cp_en)*constitutive_texVolFrac(grain,CPFEM_in,cp_en)
CPFEM_PK1_bar(:,:,CPFEM_in,cp_en) = CPFEM_PK1_bar(:,:,CPFEM_in,cp_en) + &
volfrac*PK1(:,:,grain) ! average Cauchy stress
if (updateJaco) then ! consistent tangent
do i = 1,3
do j = 1,3
CPFEM_dPdF_bar(:,:,i,j,CPFEM_in,cp_en) = CPFEM_dPdF_bar(:,:,i,j,CPFEM_in,cp_en) + &
volfrac*dPdF(:,:,i,j,grain)*dFg_cor(i,j,grain)
enddo
enddo
endif
!
! update results plotted in MENTAT
call math_pDecomposition(Fe1(:,:,grain),U,R,error) ! polar decomposition
if (error) then
write(6,*) Fe1(:,:,grain)
write(6,*) 'polar decomposition'
write(6,*) 'Grain: ',grain
write(6,*) 'Integration point: ',CPFEM_in
write(6,*) 'Element: ',mesh_element(1,cp_en)
call IO_error(650)
return
endif
CPFEM_results(1:3,grain,CPFEM_in,cp_en) = math_RtoEuler(transpose(R))*inDeg ! orientation
CPFEM_results(4 ,grain,CPFEM_in,cp_en) = volfrac ! volume fraction of orientation
enddo
!
return
!
END SUBROUTINE
!
!
!********************************************************************
! Calculates the relaxed deformation gradients of grains
!********************************************************************
subroutine GIA_RelaxedDeformation(&
F,& ! relaxed deformation gradient of grains
F_bar,& ! effective deformation gradient
r) ! relaxation vectors at boundary
!
implicit none
!
real(pReal), dimension(3,3) :: F_bar
real(pReal), dimension(3,3,8) :: F
real(pReal), dimension(3,12) :: r,n
integer(pInt) i,j,iBoun,grain
!
n = GIA_bNorm
do i = 1,3
do j = 1,3
F(i,j,1) = F_bar(i,j) + n(i, 1)*r(j, 1) + n(i, 5)*r(j, 5) + n(i, 9)*r(j, 9)
F(i,j,2) = F_bar(i,j) - n(i, 1)*r(j, 1) + n(i, 6)*r(j, 6) + n(i,10)*r(j,10)
F(i,j,3) = F_bar(i,j) + n(i, 2)*r(j, 2) - n(i, 5)*r(j, 5) + n(i,11)*r(j,11)
F(i,j,4) = F_bar(i,j) - n(i, 2)*r(j, 2) - n(i, 6)*r(j, 6) + n(i,12)*r(j,12)
F(i,j,5) = F_bar(i,j) + n(i, 3)*r(j, 3) + n(i, 7)*r(j, 7) - n(i, 9)*r(j, 9)
F(i,j,6) = F_bar(i,j) - n(i, 3)*r(j, 3) + n(i, 8)*r(j, 8) - n(i,10)*r(j,10)
F(i,j,7) = F_bar(i,j) + n(i, 4)*r(j, 4) - n(i, 7)*r(j, 7) - n(i,11)*r(j,11)
F(i,j,8) = F_bar(i,j) - n(i, 4)*r(j, 4) - n(i, 8)*r(j, 8) - n(i,12)*r(j,12)
enddo
enddo
!
return
!
END SUBROUTINE
!
!
!********************************************************************
! Calculates the jump of tensors across the grain boundary
!********************************************************************
subroutine GIA_BoundaryJump(&
F_boun,& ! tensor jump across the boundary
F_bulk) ! bulk tensor
!
implicit none
!
real(pReal), dimension(3,3,12) :: F_boun
real(pReal), dimension(3,3,8) :: F_bulk
integer(pInt) i,j,iBoun,grain
!
F_boun(:,:, 1) = F_bulk(:,:,2) - F_bulk(:,:,1)
F_boun(:,:, 2) = F_bulk(:,:,4) - F_bulk(:,:,3)
F_boun(:,:, 3) = F_bulk(:,:,6) - F_bulk(:,:,5)
F_boun(:,:, 4) = F_bulk(:,:,8) - F_bulk(:,:,7)
F_boun(:,:, 5) = F_bulk(:,:,3) - F_bulk(:,:,1)
F_boun(:,:, 6) = F_bulk(:,:,4) - F_bulk(:,:,2)
F_boun(:,:, 7) = F_bulk(:,:,7) - F_bulk(:,:,5)
F_boun(:,:, 8) = F_bulk(:,:,8) - F_bulk(:,:,6)
F_boun(:,:, 9) = F_bulk(:,:,5) - F_bulk(:,:,1)
F_boun(:,:,10) = F_bulk(:,:,6) - F_bulk(:,:,2)
F_boun(:,:,11) = F_bulk(:,:,7) - F_bulk(:,:,3)
F_boun(:,:,12) = F_bulk(:,:,8) - F_bulk(:,:,4)
!
return
!
END SUBROUTINE
!
!
!********************************************************************
! Calculates the jump of tensors across the grain boundary
!********************************************************************
subroutine GIA_JacobianMatrix(&
dresdvar,& ! Jacobian matrix
dPdF,& ! stress consistent tangent of bulk
dRdX) ! stress-like penalty tangent at boundary
!
implicit none
!
real(pReal), dimension(3,3,3,3,8) :: dPdF
real(pReal), dimension(3,3,3,3,12) :: dRdX
real(pReal), dimension(36,36) :: dresdvar
real(pReal), dimension(3,12) :: n
integer(pInt) i,j,k,l
!
n = GIA_bNorm
dresdvar = 0.0_pReal
do i = 1,3
do k = 1,3
do l = 1,3
do j = 1,3
!
! at boundary 1, influenced by boundary +5, -6, +9, -10
dresdvar(( 1-1)*3 + j,( 1-1)*3 + l) = dresdvar(( 1-1)*3 + j,( 1-1)*3 + l) &
+ (dPdF(i,j,k,l, 1) + dPdF(i,j,k,l, 2))*n(i, 1)*n(k, 1) &
+ (dRdX(i,j,k,l, 1) + dRdX(i,j,k,l, 1))*n(i, 1)*n(k, 1)
dresdvar(( 1-1)*3 + j,( 5-1)*3 + l) = dresdvar(( 1-1)*3 + j,( 5-1)*3 + l) + dPdF(i,j,k,l, 1)*n(i, 1)*n(k, 5) &
+ dRdX(i,j,k,l, 1)*n(i, 1)*n(k, 5)
dresdvar(( 1-1)*3 + j,( 6-1)*3 + l) = dresdvar(( 1-1)*3 + j,( 6-1)*3 + l) - dPdF(i,j,k,l, 2)*n(i, 1)*n(k, 6) &
- dRdX(i,j,k,l, 1)*n(i, 1)*n(k, 6)
dresdvar(( 1-1)*3 + j,( 9-1)*3 + l) = dresdvar(( 1-1)*3 + j,( 9-1)*3 + l) + dPdF(i,j,k,l, 1)*n(i, 1)*n(k, 9) &
+ dRdX(i,j,k,l, 1)*n(i, 1)*n(k, 9)
dresdvar(( 1-1)*3 + j,(10-1)*3 + l) = dresdvar(( 1-1)*3 + j,(10-1)*3 + l) - dPdF(i,j,k,l, 2)*n(i, 1)*n(k,10) &
- dRdX(i,j,k,l, 1)*n(i, 1)*n(k,10)
!
! at boundary 2, influenced by boundary -5, +6, +11, -12
dresdvar(( 2-1)*3 + j,( 2-1)*3 + l) = dresdvar(( 2-1)*3 + j,( 2-1)*3 + l) &
+ (dPdF(i,j,k,l, 3) + dPdF(i,j,k,l, 4))*n(i, 2)*n(k, 2) &
+ (dRdX(i,j,k,l, 2) + dRdX(i,j,k,l, 2))*n(i, 2)*n(k, 2)
dresdvar(( 2-1)*3 + j,( 5-1)*3 + l) = dresdvar(( 2-1)*3 + j,( 5-1)*3 + l) - dPdF(i,j,k,l, 3)*n(i, 2)*n(k, 5) &
- dRdX(i,j,k,l, 2)*n(i, 2)*n(k, 5)
dresdvar(( 2-1)*3 + j,( 6-1)*3 + l) = dresdvar(( 2-1)*3 + j,( 6-1)*3 + l) + dPdF(i,j,k,l, 4)*n(i, 2)*n(k, 6) &
+ dRdX(i,j,k,l, 2)*n(i, 2)*n(k, 6)
dresdvar(( 2-1)*3 + j,(11-1)*3 + l) = dresdvar(( 2-1)*3 + j,(11-1)*3 + l) + dPdF(i,j,k,l, 3)*n(i, 2)*n(k,11) &
+ dRdX(i,j,k,l, 2)*n(i, 2)*n(k,11)
dresdvar(( 2-1)*3 + j,(12-1)*3 + l) = dresdvar(( 2-1)*3 + j,(12-1)*3 + l) - dPdF(i,j,k,l, 4)*n(i, 2)*n(k,12) &
- dRdX(i,j,k,l, 2)*n(i, 2)*n(k,12)
!
! at boundary 3, influenced by boundary +7, -8, -9, +10
dresdvar(( 3-1)*3 + j,( 3-1)*3 + l) = dresdvar(( 3-1)*3 + j,( 3-1)*3 + l) &
+ (dPdF(i,j,k,l, 5) + dPdF(i,j,k,l, 6))*n(i, 3)*n(k, 3) &
+ (dRdX(i,j,k,l, 3) + dRdX(i,j,k,l, 3))*n(i, 3)*n(k, 3)
dresdvar(( 3-1)*3 + j,( 7-1)*3 + l) = dresdvar(( 3-1)*3 + j,( 7-1)*3 + l) + dPdF(i,j,k,l, 5)*n(i, 3)*n(k, 7) &
+ dRdX(i,j,k,l, 3)*n(i, 3)*n(k, 7)
dresdvar(( 3-1)*3 + j,( 8-1)*3 + l) = dresdvar(( 3-1)*3 + j,( 8-1)*3 + l) - dPdF(i,j,k,l, 6)*n(i, 3)*n(k, 8) &
- dRdX(i,j,k,l, 3)*n(i, 3)*n(k, 8)
dresdvar(( 3-1)*3 + j,( 9-1)*3 + l) = dresdvar(( 3-1)*3 + j,( 9-1)*3 + l) - dPdF(i,j,k,l, 5)*n(i, 3)*n(k, 9) &
- dRdX(i,j,k,l, 3)*n(i, 3)*n(k, 9)
dresdvar(( 3-1)*3 + j,(10-1)*3 + l) = dresdvar(( 3-1)*3 + j,(10-1)*3 + l) + dPdF(i,j,k,l, 6)*n(i, 3)*n(k,10) &
+ dRdX(i,j,k,l, 3)*n(i, 3)*n(k,10)
!
! at boundary 4, influenced by boundary -7, +8, -11, +12
dresdvar(( 4-1)*3 + j,( 4-1)*3 + l) = dresdvar(( 4-1)*3 + j,( 4-1)*3 + l) &
+ (dPdF(i,j,k,l, 7) + dPdF(i,j,k,l, 8))*n(i, 4)*n(k, 4) &
+ (dRdX(i,j,k,l, 4) + dRdX(i,j,k,l, 4))*n(i, 4)*n(k, 4)
dresdvar(( 4-1)*3 + j,( 7-1)*3 + l) = dresdvar(( 4-1)*3 + j,( 7-1)*3 + l) - dPdF(i,j,k,l, 7)*n(i, 4)*n(k, 7) &
- dRdX(i,j,k,l, 4)*n(i, 4)*n(k, 7)
dresdvar(( 4-1)*3 + j,( 8-1)*3 + l) = dresdvar(( 4-1)*3 + j,( 8-1)*3 + l) + dPdF(i,j,k,l, 8)*n(i, 4)*n(k, 8) &
+ dRdX(i,j,k,l, 4)*n(i, 4)*n(k, 8)
dresdvar(( 4-1)*3 + j,(11-1)*3 + l) = dresdvar(( 4-1)*3 + j,(11-1)*3 + l) - dPdF(i,j,k,l, 7)*n(i, 4)*n(k,11) &
- dRdX(i,j,k,l, 4)*n(i, 4)*n(k,11)
dresdvar(( 4-1)*3 + j,(12-1)*3 + l) = dresdvar(( 4-1)*3 + j,(12-1)*3 + l) + dPdF(i,j,k,l, 8)*n(i, 4)*n(k,12) &
+ dRdX(i,j,k,l, 4)*n(i, 4)*n(k,12)
!
! at boundary 5, influenced by boundary +1, -2, +9, -11
dresdvar(( 5-1)*3 + j,( 5-1)*3 + l) = dresdvar(( 5-1)*3 + j,( 5-1)*3 + l) &
+ (dPdF(i,j,k,l, 1) + dPdF(i,j,k,l, 3))*n(i, 5)*n(k, 5) &
+ (dRdX(i,j,k,l, 5) + dRdX(i,j,k,l, 5))*n(i, 5)*n(k, 5)
dresdvar(( 5-1)*3 + j,( 1-1)*3 + l) = dresdvar(( 5-1)*3 + j,( 1-1)*3 + l) + dPdF(i,j,k,l, 1)*n(i, 5)*n(k, 1) &
+ dRdX(i,j,k,l, 5)*n(i, 5)*n(k, 1)
dresdvar(( 5-1)*3 + j,( 2-1)*3 + l) = dresdvar(( 5-1)*3 + j,( 2-1)*3 + l) - dPdF(i,j,k,l, 3)*n(i, 5)*n(k, 2) &
- dRdX(i,j,k,l, 5)*n(i, 5)*n(k, 2)
dresdvar(( 5-1)*3 + j,( 9-1)*3 + l) = dresdvar(( 5-1)*3 + j,( 9-1)*3 + l) + dPdF(i,j,k,l, 1)*n(i, 5)*n(k, 9) &
+ dRdX(i,j,k,l, 5)*n(i, 5)*n(k, 9)
dresdvar(( 5-1)*3 + j,(11-1)*3 + l) = dresdvar(( 5-1)*3 + j,(11-1)*3 + l) - dPdF(i,j,k,l, 3)*n(i, 5)*n(k,11) &
- dRdX(i,j,k,l, 5)*n(i, 5)*n(k,11)
!
! at boundary 6, influenced by boundary -1, +2, +10, -12
dresdvar(( 6-1)*3 + j,( 6-1)*3 + l) = dresdvar(( 6-1)*3 + j,( 6-1)*3 + l) &
+ (dPdF(i,j,k,l, 2) + dPdF(i,j,k,l, 4))*n(i, 6)*n(k, 6) &
+ (dRdX(i,j,k,l, 6) + dRdX(i,j,k,l, 6))*n(i, 6)*n(k, 6)
dresdvar(( 6-1)*3 + j,( 1-1)*3 + l) = dresdvar(( 6-1)*3 + j,( 1-1)*3 + l) - dPdF(i,j,k,l, 2)*n(i, 6)*n(k, 1) &
- dRdX(i,j,k,l, 6)*n(i, 6)*n(k, 1)
dresdvar(( 6-1)*3 + j,( 2-1)*3 + l) = dresdvar(( 6-1)*3 + j,( 2-1)*3 + l) + dPdF(i,j,k,l, 4)*n(i, 6)*n(k, 2) &
+ dRdX(i,j,k,l, 6)*n(i, 6)*n(k, 2)
dresdvar(( 6-1)*3 + j,(10-1)*3 + l) = dresdvar(( 6-1)*3 + j,(10-1)*3 + l) + dPdF(i,j,k,l, 2)*n(i, 6)*n(k,10) &
+ dRdX(i,j,k,l, 6)*n(i, 6)*n(k,10)
dresdvar(( 6-1)*3 + j,(12-1)*3 + l) = dresdvar(( 6-1)*3 + j,(12-1)*3 + l) - dPdF(i,j,k,l, 4)*n(i, 6)*n(k,12) &
- dRdX(i,j,k,l, 6)*n(i, 6)*n(k,12)
!
! at boundary 7, influenced by boundary +3, -4, -9, +11
dresdvar(( 7-1)*3 + j,( 7-1)*3 + l) = dresdvar(( 7-1)*3 + j,( 7-1)*3 + l) &
+ (dPdF(i,j,k,l, 5) + dPdF(i,j,k,l, 7))*n(i, 7)*n(k, 7) &
+ (dRdX(i,j,k,l, 7) + dRdX(i,j,k,l, 7))*n(i, 7)*n(k, 7)
dresdvar(( 7-1)*3 + j,( 3-1)*3 + l) = dresdvar(( 7-1)*3 + j,( 3-1)*3 + l) + dPdF(i,j,k,l, 5)*n(i, 7)*n(k, 3) &
+ dRdX(i,j,k,l, 7)*n(i, 7)*n(k, 3)
dresdvar(( 7-1)*3 + j,( 4-1)*3 + l) = dresdvar(( 7-1)*3 + j,( 4-1)*3 + l) - dPdF(i,j,k,l, 7)*n(i, 7)*n(k, 4) &
- dRdX(i,j,k,l, 7)*n(i, 7)*n(k, 4)
dresdvar(( 7-1)*3 + j,( 9-1)*3 + l) = dresdvar(( 7-1)*3 + j,( 9-1)*3 + l) - dPdF(i,j,k,l, 5)*n(i, 7)*n(k, 9) &
- dRdX(i,j,k,l, 7)*n(i, 7)*n(k, 9)
dresdvar(( 7-1)*3 + j,(11-1)*3 + l) = dresdvar(( 7-1)*3 + j,(11-1)*3 + l) + dPdF(i,j,k,l, 7)*n(i, 7)*n(k,11) &
+ dRdX(i,j,k,l, 7)*n(i, 7)*n(k,11)
!
! at boundary 8, influenced by boundary -3, +4, -10, +12
dresdvar(( 8-1)*3 + j,( 8-1)*3 + l) = dresdvar(( 8-1)*3 + j,( 8-1)*3 + l) &
+ (dPdF(i,j,k,l, 6) + dPdF(i,j,k,l, 8))*n(i, 8)*n(k, 8) &
+ (dRdX(i,j,k,l, 8) + dRdX(i,j,k,l, 8))*n(i, 8)*n(k, 8)
dresdvar(( 8-1)*3 + j,( 3-1)*3 + l) = dresdvar(( 8-1)*3 + j,( 3-1)*3 + l) - dPdF(i,j,k,l, 6)*n(i, 8)*n(k, 3) &
- dRdX(i,j,k,l, 8)*n(i, 8)*n(k, 3)
dresdvar(( 8-1)*3 + j,( 4-1)*3 + l) = dresdvar(( 8-1)*3 + j,( 4-1)*3 + l) + dPdF(i,j,k,l, 8)*n(i, 8)*n(k, 4) &
+ dRdX(i,j,k,l, 8)*n(i, 8)*n(k, 4)
dresdvar(( 8-1)*3 + j,(10-1)*3 + l) = dresdvar(( 8-1)*3 + j,(10-1)*3 + l) - dPdF(i,j,k,l, 6)*n(i, 8)*n(k,10) &
- dRdX(i,j,k,l, 8)*n(i, 8)*n(k,10)
dresdvar(( 8-1)*3 + j,(12-1)*3 + l) = dresdvar(( 8-1)*3 + j,(12-1)*3 + l) + dPdF(i,j,k,l, 8)*n(i, 8)*n(k,12) &
+ dRdX(i,j,k,l, 8)*n(i, 8)*n(k,12)
!
! at boundary 9, influenced by boundary +1, -3, +5, -7
dresdvar(( 9-1)*3 + j,( 9-1)*3 + l) = dresdvar(( 9-1)*3 + j,( 9-1)*3 + l) &
+ (dPdF(i,j,k,l, 1) + dPdF(i,j,k,l, 5))*n(i, 9)*n(k, 9) &
+ (dRdX(i,j,k,l, 9) + dRdX(i,j,k,l, 9))*n(i, 9)*n(k, 9)
dresdvar(( 9-1)*3 + j,( 1-1)*3 + l) = dresdvar(( 9-1)*3 + j,( 1-1)*3 + l) + dPdF(i,j,k,l, 1)*n(i, 9)*n(k, 1) &
+ dRdX(i,j,k,l, 9)*n(i, 9)*n(k, 1)
dresdvar(( 9-1)*3 + j,( 3-1)*3 + l) = dresdvar(( 9-1)*3 + j,( 3-1)*3 + l) - dPdF(i,j,k,l, 5)*n(i, 9)*n(k, 3) &
- dRdX(i,j,k,l, 9)*n(i, 9)*n(k, 3)
dresdvar(( 9-1)*3 + j,( 5-1)*3 + l) = dresdvar(( 9-1)*3 + j,( 5-1)*3 + l) + dPdF(i,j,k,l, 1)*n(i, 9)*n(k, 5) &
+ dRdX(i,j,k,l, 9)*n(i, 9)*n(k, 5)
dresdvar(( 9-1)*3 + j,( 7-1)*3 + l) = dresdvar(( 9-1)*3 + j,( 7-1)*3 + l) - dPdF(i,j,k,l, 5)*n(i, 9)*n(k, 7) &
- dRdX(i,j,k,l, 9)*n(i, 9)*n(k, 7)
!
! at boundary 10, influenced by boundary -1, +3, +6, -8
dresdvar((10-1)*3 + j,(10-1)*3 + l) = dresdvar((10-1)*3 + j,(10-1)*3 + l) &
+ (dPdF(i,j,k,l, 2) + dPdF(i,j,k,l, 6))*n(i,10)*n(k,10) &
+ (dRdX(i,j,k,l,10) + dRdX(i,j,k,l,10))*n(i,10)*n(k,10)
dresdvar((10-1)*3 + j,( 1-1)*3 + l) = dresdvar((10-1)*3 + j,( 1-1)*3 + l) - dPdF(i,j,k,l, 2)*n(i,10)*n(k, 1) &
- dRdX(i,j,k,l,10)*n(i,10)*n(k, 1)
dresdvar((10-1)*3 + j,( 3-1)*3 + l) = dresdvar((10-1)*3 + j,( 3-1)*3 + l) + dPdF(i,j,k,l, 6)*n(i,10)*n(k, 3) &
+ dRdX(i,j,k,l,10)*n(i,10)*n(k, 3)
dresdvar((10-1)*3 + j,( 6-1)*3 + l) = dresdvar((10-1)*3 + j,( 6-1)*3 + l) + dPdF(i,j,k,l, 2)*n(i,10)*n(k, 6) &
+ dRdX(i,j,k,l,10)*n(i,10)*n(k, 6)
dresdvar((10-1)*3 + j,( 8-1)*3 + l) = dresdvar((10-1)*3 + j,( 8-1)*3 + l) - dPdF(i,j,k,l, 6)*n(i,10)*n(k, 8) &
- dRdX(i,j,k,l,10)*n(i,10)*n(k, 8)
!
! at boundary 11, influenced by boundary +2, -4, -5, +7
dresdvar((11-1)*3 + j,(11-1)*3 + l) = dresdvar((11-1)*3 + j,(11-1)*3 + l) &
+ (dPdF(i,j,k,l, 3) + dPdF(i,j,k,l, 7))*n(i,11)*n(k,11) &
+ (dRdX(i,j,k,l,11) + dRdX(i,j,k,l,11))*n(i,11)*n(k,11)
dresdvar((11-1)*3 + j,( 2-1)*3 + l) = dresdvar((11-1)*3 + j,( 2-1)*3 + l) + dPdF(i,j,k,l, 3)*n(i,11)*n(k, 2) &
+ dRdX(i,j,k,l,11)*n(i,11)*n(k, 2)
dresdvar((11-1)*3 + j,( 4-1)*3 + l) = dresdvar((11-1)*3 + j,( 4-1)*3 + l) - dPdF(i,j,k,l, 7)*n(i,11)*n(k, 4) &
- dRdX(i,j,k,l,11)*n(i,11)*n(k, 4)
dresdvar((11-1)*3 + j,( 5-1)*3 + l) = dresdvar((11-1)*3 + j,( 5-1)*3 + l) - dPdF(i,j,k,l, 3)*n(i,11)*n(k, 5) &
- dRdX(i,j,k,l,11)*n(i,11)*n(k, 5)
dresdvar((11-1)*3 + j,( 7-1)*3 + l) = dresdvar((11-1)*3 + j,( 7-1)*3 + l) + dPdF(i,j,k,l, 7)*n(i,11)*n(k, 7) &
+ dRdX(i,j,k,l,11)*n(i,11)*n(k, 7)
!
! at boundary 12, influenced by boundary -2, +4, -6, +8
dresdvar((12-1)*3 + j,(12-1)*3 + l) = dresdvar((12-1)*3 + j,(12-1)*3 + l) &
+ (dPdF(i,j,k,l, 4) + dPdF(i,j,k,l, 8))*n(i,12)*n(k,12) &
+ (dRdX(i,j,k,l,12) + dRdX(i,j,k,l,12))*n(i,12)*n(k,12)
dresdvar((12-1)*3 + j,( 2-1)*3 + l) = dresdvar((12-1)*3 + j,( 2-1)*3 + l) - dPdF(i,j,k,l, 4)*n(i,12)*n(k, 2) &
- dRdX(i,j,k,l,12)*n(i,12)*n(k, 2)
dresdvar((12-1)*3 + j,( 4-1)*3 + l) = dresdvar((12-1)*3 + j,( 4-1)*3 + l) + dPdF(i,j,k,l, 8)*n(i,12)*n(k, 4) &
+ dRdX(i,j,k,l,12)*n(i,12)*n(k, 4)
dresdvar((12-1)*3 + j,( 6-1)*3 + l) = dresdvar((12-1)*3 + j,( 6-1)*3 + l) - dPdF(i,j,k,l, 4)*n(i,12)*n(k, 6) &
- dRdX(i,j,k,l,12)*n(i,12)*n(k, 6)
dresdvar((12-1)*3 + j,( 8-1)*3 + l) = dresdvar((12-1)*3 + j,( 8-1)*3 + l) + dPdF(i,j,k,l, 8)*n(i,12)*n(k, 8) &
+ dRdX(i,j,k,l,12)*n(i,12)*n(k, 8)
!
enddo
enddo
enddo
enddo
!
return
!
END SUBROUTINE
!
!
!********************************************************************
! Calculate the correction for the effective consisten tangent
!********************************************************************
subroutine GIA_TangentCorrection(&
dF_grain,& ! deformation gradient increment of grains
ddF_corr) !
!
implicit none
!
real(pReal), dimension(3,3,8) :: dF_grain,ddF_corr
integer(pInt) i,j,iBoun,grain
!
ddF_corr = 0.0_pReal
do j = 1,3
!
! first relaxation direction
ddF_corr(1,j,1) = 2.0_pReal*abs(dF_grain(1,j,1))/(abs(dF_grain(1,j,1)) + abs(dF_grain(1,j,2)))
if (abs(dF_grain(1,j,1)) < 1.0e-8_pReal) ddF_corr(1,j,1) = 1.0_pReal
ddF_corr(1,j,2) = 2.0_pReal - ddF_corr(1,j,1)
ddF_corr(1,j,3) = 2.0_pReal*abs(dF_grain(1,j,3))/(abs(dF_grain(1,j,3)) + abs(dF_grain(1,j,4)))
if (abs(dF_grain(1,j,3)) < 1.0e-8_pReal) ddF_corr(1,j,3) = 1.0_pReal
ddF_corr(1,j,4) = 2.0_pReal - ddF_corr(1,j,3)
ddF_corr(1,j,5) = 2.0_pReal*abs(dF_grain(1,j,5))/(abs(dF_grain(1,j,5)) + abs(dF_grain(1,j,6)))
if (abs(dF_grain(1,j,5)) < 1.0e-8_pReal) ddF_corr(1,j,5) = 1.0_pReal
ddF_corr(1,j,6) = 2.0_pReal - ddF_corr(1,j,5)
ddF_corr(1,j,7) = 2.0_pReal*abs(dF_grain(1,j,7))/(abs(dF_grain(1,j,7)) + abs(dF_grain(1,j,8)))
if (abs(dF_grain(1,j,7)) < 1.0e-8_pReal) ddF_corr(1,j,7) = 1.0_pReal
ddF_corr(1,j,8) = 2.0_pReal - ddF_corr(1,j,7)
!
! second relaxation direction
ddF_corr(2,j,1) = 2.0_pReal*abs(dF_grain(2,j,1))/(abs(dF_grain(2,j,1)) + abs(dF_grain(2,j,3)))
if (abs(dF_grain(2,j,1)) < 1.0e-8_pReal) ddF_corr(2,j,1) = 1.0_pReal
ddF_corr(2,j,2) = 2.0_pReal*abs(dF_grain(2,j,2))/(abs(dF_grain(2,j,2)) + abs(dF_grain(2,j,4)))
if (abs(dF_grain(2,j,2)) < 1.0e-8_pReal) ddF_corr(2,j,2) = 1.0_pReal
ddF_corr(2,j,3) = 2.0_pReal - ddF_corr(2,j,1)
ddF_corr(2,j,4) = 2.0_pReal - ddF_corr(2,j,2)
ddF_corr(2,j,5) = 2.0_pReal*abs(dF_grain(2,j,5))/(abs(dF_grain(2,j,5)) + abs(dF_grain(2,j,7)))
if (abs(dF_grain(2,j,5)) < 1.0e-8_pReal) ddF_corr(2,j,5) = 1.0_pReal
ddF_corr(2,j,6) = 2.0_pReal*abs(dF_grain(2,j,6))/(abs(dF_grain(2,j,6)) + abs(dF_grain(2,j,8)))
if (abs(dF_grain(2,j,6)) < 1.0e-8_pReal) ddF_corr(2,j,6) = 1.0_pReal
ddF_corr(2,j,7) = 2.0_pReal - ddF_corr(2,j,5)
ddF_corr(2,j,8) = 2.0_pReal - ddF_corr(2,j,6)
!
! third relaxation direction
ddF_corr(3,j,1) = 2.0_pReal*abs(dF_grain(3,j,1))/(abs(dF_grain(3,j,1)) + abs(dF_grain(3,j,5)))
if (abs(dF_grain(3,j,1)) < 1.0e-8_pReal) ddF_corr(3,j,1) = 1.0_pReal
ddF_corr(3,j,2) = 2.0_pReal*abs(dF_grain(3,j,2))/(abs(dF_grain(3,j,2)) + abs(dF_grain(3,j,6)))
if (abs(dF_grain(3,j,2)) < 1.0e-8_pReal) ddF_corr(3,j,2) = 1.0_pReal
ddF_corr(3,j,3) = 2.0_pReal*abs(dF_grain(3,j,3))/(abs(dF_grain(3,j,3)) + abs(dF_grain(3,j,7)))
if (abs(dF_grain(3,j,3)) < 1.0e-8_pReal) ddF_corr(3,j,3) = 1.0_pReal
ddF_corr(3,j,4) = 2.0_pReal*abs(dF_grain(3,j,4))/(abs(dF_grain(3,j,4)) + abs(dF_grain(3,j,8)))
if (abs(dF_grain(3,j,4)) < 1.0e-8_pReal) ddF_corr(3,j,4) = 1.0_pReal
ddF_corr(3,j,5) = 2.0_pReal - ddF_corr(3,j,1)
ddF_corr(3,j,6) = 2.0_pReal - ddF_corr(3,j,2)
ddF_corr(3,j,7) = 2.0_pReal - ddF_corr(3,j,3)
ddF_corr(3,j,8) = 2.0_pReal - ddF_corr(3,j,4)
!
enddo
!
return
!
END SUBROUTINE
!
END MODULE
!##############################################################

290
trunk/CPFEM_Taylor.f90 Normal file
View File

@ -0,0 +1,290 @@
!##############################################################
MODULE CPFEM
!##############################################################
! *** CPFEM engine ***
!
use prec, only: pReal,pInt
implicit none
!
! ****************************************************************
! *** General variables for the material behaviour calculation ***
! ****************************************************************
real(pReal), dimension (:,:), allocatable :: CPFEM_Temperature
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_ffn_bar
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_ffn1_bar
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_PK1_bar
real(pReal), dimension (:,:,:,:,:,:),allocatable :: CPFEM_dPdF_bar
real(pReal), dimension (:,:,:), allocatable :: CPFEM_stress_bar
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_jaco_bar
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_jaco_knownGood
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_results
real(pReal), dimension (:,:,:,:,:), allocatable :: CPFEM_Fp_old
real(pReal), dimension (:,:,:,:,:), allocatable :: CPFEM_Fp_new
real(pReal), parameter :: CPFEM_odd_stress = 1e15_pReal, CPFEM_odd_jacobian = 1e50_pReal
integer(pInt) :: CPFEM_Nresults = 4_pInt ! three Euler angles plus volume fraction
logical :: CPFEM_init_done = .false. ! remember if init has been done already
logical :: CPFEM_calc_done = .false. ! remember if first IP has already calced the results
!
CONTAINS
!
!*********************************************************
!*** allocate the arrays defined in module CPFEM ***
!*** and initialize them ***
!*********************************************************
SUBROUTINE CPFEM_init(Temperature)
!
use prec
use math, only: math_EulertoR, math_I3, math_identity2nd
use mesh
use constitutive
!
implicit none
!
real(pReal) Temperature
integer(pInt) e,i,g
!
! *** mpie.marc parameters ***
allocate(CPFEM_Temperature (mesh_maxNips,mesh_NcpElems)) ; CPFEM_Temperature = Temperature
allocate(CPFEM_ffn_bar (3,3,mesh_maxNips,mesh_NcpElems))
forall(e=1:mesh_NcpElems,i=1:mesh_maxNips) CPFEM_ffn_bar(:,:,i,e) = math_I3
allocate(CPFEM_ffn1_bar (3,3,mesh_maxNips,mesh_NcpElems)) ; CPFEM_ffn1_bar = CPFEM_ffn_bar
allocate(CPFEM_PK1_bar (3,3,mesh_maxNips,mesh_NcpElems)) ; CPFEM_PK1_bar = 0.0_pReal
allocate(CPFEM_dPdF_bar(3,3,3,3,mesh_maxNips,mesh_NcpElems)) ; CPFEM_dPdF_bar = 0.0_pReal
allocate(CPFEM_stress_bar(6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_stress_bar = 0.0_pReal
allocate(CPFEM_jaco_bar(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_jaco_bar = 0.0_pReal
allocate(CPFEM_jaco_knownGood(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_jaco_knownGood = 0.0_pReal
!
! *** User defined results !!! MISSING incorporate consti_Nresults ***
allocate(CPFEM_results(CPFEM_Nresults+constitutive_maxNresults,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems))
CPFEM_results = 0.0_pReal
!
! *** Plastic deformation gradient at (t=t0) and (t=t1) ***
allocate(CPFEM_Fp_new(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) ; CPFEM_Fp_new = 0.0_pReal
allocate(CPFEM_Fp_old(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems))
forall (e=1:mesh_NcpElems,i=1:mesh_maxNips,g=1:constitutive_maxNgrains) &
CPFEM_Fp_old(:,:,g,i,e) = math_EulerToR(constitutive_EulerAngles(:,g,i,e)) ! plastic def gradient reflects init orientation
!
! *** Output to MARC output file ***
write(6,*)
write(6,*) 'CPFEM Initialization'
write(6,*)
write(6,*) 'CPFEM_Temperature: ', shape(CPFEM_Temperature)
write(6,*) 'CPFEM_ffn_bar: ', shape(CPFEM_ffn_bar)
write(6,*) 'CPFEM_ffn1_bar: ', shape(CPFEM_ffn1_bar)
write(6,*) 'CPFEM_PK1_bar: ', shape(CPFEM_PK1_bar)
write(6,*) 'CPFEM_dPdF_bar: ', shape(CPFEM_dPdF_bar)
write(6,*) 'CPFEM_stress_bar: ', shape(CPFEM_stress_bar)
write(6,*) 'CPFEM_jaco_bar: ', shape(CPFEM_jaco_bar)
write(6,*) 'CPFEM_jaco_knownGood: ', shape(CPFEM_jaco_knownGood)
write(6,*) 'CPFEM_results: ', shape(CPFEM_results)
write(6,*) 'CPFEM_Fp_old: ', shape(CPFEM_Fp_old)
write(6,*) 'CPFEM_Fp_new: ', shape(CPFEM_Fp_new)
write(6,*)
call flush(6)
return
!
END SUBROUTINE
!
!
!***********************************************************************
!*** perform initialization at first call, update variables and ***
!*** call the actual material model ***
!
! CPFEM_mode computation mode (regular, collection, recycle)
! ffn deformation gradient for t=t0
! ffn1 deformation gradient for t=t1
! Temperature temperature
! CPFEM_dt time increment
! CPFEM_en element number
! CPFEM_in intergration point number
! CPFEM_stress stress vector in Mandel notation
! CPFEM_updateJaco flag to initiate computation of Jacobian
! CPFEM_jaco jacobian in Mandel notation
! CPFEM_ngens size of stress strain law
!***********************************************************************
SUBROUTINE CPFEM_general(CPFEM_mode, ffn, ffn1, Temperature, CPFEM_dt,&
CPFEM_en, CPFEM_in, CPFEM_stress, CPFEM_updateJaco, CPFEM_jaco, CPFEM_ngens)
! note: CPFEM_stress = Cauchy stress cs(6) and CPFEM_jaco = Consistent tangent dcs/de
!
use prec, only: pReal,pInt
use FEsolving
use debug
use math, only: math_init, invnrmMandel, math_identity2nd, math_Mandel3333to66,math_Mandel33to6,math_Mandel6to33,math_det3x3,math_I3
use mesh, only: mesh_init,mesh_FEasCP, mesh_NcpElems, FE_Nips, FE_mapElemtype, mesh_element
use lattice, only: lattice_init
use constitutive, only: constitutive_init,constitutive_state_old,constitutive_state_new,material_Cslip_66
implicit none
!
integer(pInt) CPFEM_en, CPFEM_in, cp_en, CPFEM_ngens, i,j,k,l,m,n, e
real(pReal), dimension (3,3) :: ffn,ffn1,Kirchhoff_bar
real(pReal), dimension (3,3,3,3) :: H_bar
real(pReal), dimension(CPFEM_ngens) :: CPFEM_stress
real(pReal), dimension(CPFEM_ngens,CPFEM_ngens) :: CPFEM_jaco
real(pReal) Temperature,CPFEM_dt,J_inverse
integer(pInt) CPFEM_mode ! 1: regular computation with aged results&
! 2: regular computation&
! 3: collection of FEM data&
! 4: recycling of former results (MARC speciality)&
! 5: record tangent from former converged inc&
! 6: restore tangent from former converged inc
logical CPFEM_updateJaco
!
if (.not. CPFEM_init_done) then ! initialization step (three dimensional stress state check missing?)
call math_init()
call mesh_init()
call lattice_init()
call constitutive_init()
call CPFEM_init(Temperature)
CPFEM_init_done = .true.
endif
!
cp_en = mesh_FEasCP('elem',CPFEM_en)
if (cp_en == 1 .and. CPFEM_in == 1) &
write(6,'(a6,x,i4,x,a4,x,i4,x,a10,x,f8.4,x,a10,x,i2,x,a10,x,i2,x,a10,x,i2,x,a10,x,i2)') &
'elem',cp_en,'IP',CPFEM_in,&
'theTime',theTime,'theInc',theInc,'theCycle',theCycle,'theLovl',theLovl,&
'mode',CPFEM_mode
!
select case (CPFEM_mode)
case (2,1) ! regular computation (with aging of results)
if (.not. CPFEM_calc_done) then ! puuh, me needs doing all the work...
write (6,*) 'puuh me needs doing all the work', cp_en
if (CPFEM_mode == 1) then ! age results at start of new increment
CPFEM_Fp_old = CPFEM_Fp_new
constitutive_state_old = constitutive_state_new
write (6,*) '#### aged results'
endif
debug_cutbackDistribution = 0_pInt ! initialize debugging data
debug_InnerLoopDistribution = 0_pInt
debug_OuterLoopDistribution = 0_pInt
!
do e=1,mesh_NcpElems ! ## this shall be done in a parallel loop in the future ##
do i=1,FE_Nips(mesh_element(2,e)) ! iterate over all IPs of this element's type
debugger = (e==1 .and. i==1) ! switch on debugging for first IP in first element
call CPFEM_MaterialPoint(CPFEM_updateJaco, CPFEM_dt, i, e)
enddo
enddo
call debug_info() ! output of debugging/performance statistics
CPFEM_calc_done = .true. ! now calc is done
endif
! translate from P and dP/dF to CS and dCS/dE
Kirchhoff_bar = matmul(CPFEM_PK1_bar(:,:,CPFEM_in, cp_en),transpose(CPFEM_ffn1_bar(:,:,CPFEM_in, cp_en)))
J_inverse = 1.0_pReal/math_det3x3(CPFEM_ffn1_bar(:,:,CPFEM_in, cp_en))
CPFEM_stress_bar(1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel33to6(J_inverse*Kirchhoff_bar)
!
H_bar = 0.0_pReal
forall(i=1:3,j=1:3,k=1:3,l=1:3,m=1:3,n=1:3) &
H_bar(i,j,k,l) = H_bar(i,j,k,l) + &
(CPFEM_ffn1_bar(j,m,CPFEM_in,cp_en)*CPFEM_ffn1_bar(l,n,CPFEM_in,cp_en)*CPFEM_dPdF_bar(i,m,k,n,CPFEM_in,cp_en) - &
math_I3(j,l)*CPFEM_ffn1_bar(i,m,CPFEM_in,cp_en)*CPFEM_PK1_bar(k,m,CPFEM_in,cp_en)) + &
0.5_pReal*(math_I3(i,k)*Kirchhoff_bar(j,l) + math_I3(j,l)*Kirchhoff_bar(i,k) + &
math_I3(i,l)*Kirchhoff_bar(j,k) + math_I3(j,k)*Kirchhoff_bar(i,l))
CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = math_Mandel3333to66(J_inverse*H_bar)
!
case (3) ! collect and return odd result
CPFEM_Temperature(CPFEM_in,cp_en) = Temperature
CPFEM_ffn_bar(:,:,CPFEM_in,cp_en) = ffn
CPFEM_ffn1_bar(:,:,CPFEM_in,cp_en) = ffn1
CPFEM_stress_bar(1:CPFEM_ngens,CPFEM_in,cp_en) = CPFEM_odd_stress
CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en) = CPFEM_odd_jacobian*math_identity2nd(CPFEM_ngens)
CPFEM_calc_done = .false.
case (4) ! do nothing since we can recycle the former results (MARC specialty)
case (5) ! record consistent tangent at beginning of new increment
CPFEM_jaco_knownGood = CPFEM_jaco_bar
case (6) ! restore consistent tangent after cutback
CPFEM_jaco_bar = CPFEM_jaco_knownGood
end select
!
! return the local stress and the jacobian from storage
CPFEM_stress(1:CPFEM_ngens) = CPFEM_stress_bar(1:CPFEM_ngens,CPFEM_in,cp_en)
CPFEM_jaco(1:CPFEM_ngens,1:CPFEM_ngens) = CPFEM_jaco_bar(1:CPFEM_ngens,1:CPFEM_ngens,CPFEM_in,cp_en)
if (cp_en == 1 .and. CPFEM_in == 1) write (6,*) 'stress',CPFEM_stress
if (cp_en == 1 .and. CPFEM_in == 1 .and. CPFEM_updateJaco) write (6,*) 'stiffness',CPFEM_jaco
!
return
!
END SUBROUTINE
!
!**********************************************************
!*** calculate the material point behaviour ***
!**********************************************************
SUBROUTINE CPFEM_MaterialPoint(&
updateJaco,& ! flag to initiate Jacobian updating
CPFEM_dt,& ! Time increment (dt)
CPFEM_in,& ! Integration point number
cp_en) ! Element number
!
use prec
use debug
use math, only: math_pDecomposition,math_RtoEuler,inDeg,math_I3,math_invert3x3,math_permut,math_invert,math_delta
use IO, only: IO_error
use mesh, only: mesh_element
use crystallite
use constitutive
implicit none
!
character(len=128) msg
integer(pInt) cp_en,CPFEM_in,grain,max_cutbacks,i,j,k,l,m,n
logical updateJaco,error
real(pReal) CPFEM_dt,volfrac
real(pReal), dimension(3,3) :: F0_bar,F1_bar,dF_bar
real(pReal), dimension(3,3) :: U,R
real(pReal), dimension(3,3,8) :: PK1,Fp0,Fp1,Fe1
real(pReal), dimension(3,3,3,3,8) :: dPdF
real(pReal), dimension(constitutive_maxNstatevars,8) :: state0,state1
!
CPFEM_PK1_bar(:,:,CPFEM_in,cp_en) = 0.0_pReal ! zero out average first PK stress
if (updateJaco) CPFEM_dPdF_bar(:,:,:,:,CPFEM_in,cp_en) = 0.0_pReal ! zero out average consistent tangent
!
F0_bar = CPFEM_ffn_bar(:,:,CPFEM_in,cp_en) ! effective deformation gradient at t_n
F1_bar = CPFEM_ffn1_bar(:,:,CPFEM_in,cp_en) ! effective deformation gradient at t_n+1
state0 = constitutive_state_old(:,:,CPFEM_in,cp_en) ! state variables at t_n
Fp0 = CPFEM_Fp_old(:,:,:,CPFEM_in,cp_en) ! grain plastic def. gradient at t_n
!
! -------------- grain loop -----------------
do grain = 1,texture_Ngrains(mesh_element(4,cp_en))
call SingleCrystallite(msg,PK1(:,:,grain),dPdF(:,:,:,:,grain),&
CPFEM_results(5:4+constitutive_Nresults(grain,CPFEM_in,cp_en),grain,CPFEM_in,cp_en),&
Fp1(:,:,grain),Fe1(:,:,grain),state1(:,grain),& ! output up to here
CPFEM_dt,cp_en,CPFEM_in,grain,.true.,&
CPFEM_Temperature(CPFEM_in,cp_en),F1_bar,F0_bar,Fp0(:,:,grain),state0(:,grain))
if (msg /= 'ok') then ! solution not reached
call IO_error(600)
return
endif
! update results plotted in MENTAT
call math_pDecomposition(Fe1(:,:,grain),U,R,error) ! polar decomposition
if (error) then
write(6,*) Fe1(:,:,grain)
write(6,*) 'polar decomposition'
write(6,*) 'Grain: ',grain
write(6,*) 'Integration point: ',CPFEM_in
write(6,*) 'Element: ',mesh_element(1,cp_en)
call IO_error(650)
return
endif
!
volfrac = constitutive_matVolFrac(grain,CPFEM_in,cp_en)*constitutive_texVolFrac(grain,CPFEM_in,cp_en)
CPFEM_PK1_bar(:,:,CPFEM_in,cp_en) = CPFEM_PK1_bar(:,:,CPFEM_in,cp_en) + &
volfrac*PK1(:,:,grain) ! average Cauchy stress
if (updateJaco) then ! consistent tangent
CPFEM_dPdF_bar(:,:,:,:,CPFEM_in,cp_en) = CPFEM_dPdF_bar(:,:,:,:,CPFEM_in,cp_en) + &
volfrac*dPdF(:,:,:,:,grain)
endif
CPFEM_results(1:3,grain,CPFEM_in,cp_en) = math_RtoEuler(transpose(R))*inDeg ! orientation
CPFEM_results(4 ,grain,CPFEM_in,cp_en) = volfrac ! volume fraction of orientation
enddo ! grain loop
!
! updates all variables, deformation gradients, and vectors
CPFEM_Fp_new(:,:,:,CPFEM_in,cp_en) = Fp1
constitutive_state_new(:,:,CPFEM_in,cp_en) = state1
!
return
!
END SUBROUTINE
!
END MODULE
!##############################################################

7
trunk/FEsolving.f90
View File

@ -3,13 +3,12 @@
MODULE FEsolving MODULE FEsolving
!############################################################## !##############################################################
use prec, only: pInt use prec, only: pInt,pReal
implicit none implicit none
integer(pInt) cycleCounter integer(pInt) cycleCounter
integer(pInt) theInc,theCycle,theLovl integer(pInt) theInc,theCycle,theLovl
logical :: outdatedByNewInc = .false. real(pReal) theTime
logical :: lastIncConverged = .false.,outdatedByNewInc = .false.
END MODULE FEsolving END MODULE FEsolving

2
trunk/IO.f90
View File

@ -595,6 +595,8 @@
msg='Polar decomposition failed' msg='Polar decomposition failed'
case (700) case (700)
msg='Singular matrix in stress iteration' msg='Singular matrix in stress iteration'
case (800)
msg='GIA requires 8 grains per IP (bonehead, you!)'
case default case default
msg='Unknown error number' msg='Unknown error number'
end select end select

30
trunk/constitutive_pheno.f90
View File

@ -271,7 +271,7 @@ do while(.true.)
else else
if (section>0) then if (section>0) then
select case(tag) select case(tag)
case ('crystal_structure') case ('lattice_structure')
material_CrystalStructure(section)=IO_intValue(line,positions,2) material_CrystalStructure(section)=IO_intValue(line,positions,2)
case ('nslip') case ('nslip')
material_Nslip(section)=IO_intValue(line,positions,2) material_Nslip(section)=IO_intValue(line,positions,2)
@ -413,7 +413,7 @@ subroutine constitutive_Parse_MatTexDat(filename)
use prec, only: pReal,pInt use prec, only: pReal,pInt
use IO, only: IO_error, IO_open_file use IO, only: IO_error, IO_open_file
use math, only: math_Mandel3333to66, math_Voigt66to3333 use math, only: math_Mandel3333to66, math_Voigt66to3333
use crystal, only: crystal_SlipIntType use lattice, only: lattice_SlipIntType
implicit none implicit none
!* Definition of variables !* Definition of variables
@ -454,7 +454,7 @@ allocate(material_n_slip(material_maxN)) ; material_n_slip=0.0_pReal
allocate(material_h0(material_maxN)) ; material_h0=0.0_pReal allocate(material_h0(material_maxN)) ; material_h0=0.0_pReal
allocate(material_s_sat(material_maxN)) ; material_s_sat=0.0_pReal allocate(material_s_sat(material_maxN)) ; material_s_sat=0.0_pReal
allocate(material_w0(material_maxN)) ; material_w0=0.0_pReal allocate(material_w0(material_maxN)) ; material_w0=0.0_pReal
allocate(material_SlipIntCoeff(maxval(crystal_SlipIntType),material_maxN)) ; material_SlipIntCoeff=0.0_pReal allocate(material_SlipIntCoeff(maxval(lattice_SlipIntType),material_maxN)) ; material_SlipIntCoeff=0.0_pReal
allocate(material_GrainSize(material_maxN)) ; material_GrainSize=0.0_pReal allocate(material_GrainSize(material_maxN)) ; material_GrainSize=0.0_pReal
allocate(material_bg(material_maxN)) ; material_bg=0.0_pReal allocate(material_bg(material_maxN)) ; material_bg=0.0_pReal
allocate(texture_ODFfile(texture_maxN)) ; texture_ODFfile='' allocate(texture_ODFfile(texture_maxN)) ; texture_ODFfile=''
@ -546,7 +546,7 @@ use prec, only: pReal,pInt
use math, only: math_sampleGaussOri,math_sampleFiberOri,math_sampleRandomOri,math_symmetricEulers,math_EulerToR use math, only: math_sampleGaussOri,math_sampleFiberOri,math_sampleRandomOri,math_symmetricEulers,math_EulerToR
use mesh, only: mesh_NcpElems,FE_Nips,mesh_maxNips,mesh_element use mesh, only: mesh_NcpElems,FE_Nips,mesh_maxNips,mesh_element
use IO, only: IO_hybridIA use IO, only: IO_hybridIA
use crystal, only: crystal_SlipIntType use lattice, only: lattice_SlipIntType
implicit none implicit none
@ -688,7 +688,7 @@ do i=1,material_maxN
do j=1,material_Nslip(i) do j=1,material_Nslip(i)
do k=1,material_Nslip(i) do k=1,material_Nslip(i)
!* min function is used to distinguish self hardening from latent hardening !* min function is used to distinguish self hardening from latent hardening
constitutive_HardeningMatrix(k,j,i) = material_SlipIntCoeff(max(2,min(3,crystal_SlipIntType(k,j,i)))-1,i) ! 1,2,3,4,5 --> 1,1,2,2,2 constitutive_HardeningMatrix(k,j,i) = material_SlipIntCoeff(max(2,min(3,lattice_SlipIntType(k,j,i)))-1,i) ! 1,2,3,4,5 --> 1,1,2,2,2
enddo enddo
enddo enddo
enddo enddo
@ -759,7 +759,7 @@ subroutine constitutive_LpAndItsTangent(Lp,dLp_dTstar,Tstar_v,state,Temperature,
!* - dLp_dTstar : derivative of Lp (4th-order tensor) * !* - dLp_dTstar : derivative of Lp (4th-order tensor) *
!********************************************************************* !*********************************************************************
use prec, only: pReal,pInt use prec, only: pReal,pInt
use crystal, only: crystal_Sslip,crystal_Sslip_v use lattice, only: lattice_Sslip,lattice_Sslip_v
use math, only: math_Plain3333to99 use math, only: math_Plain3333to99
use debug use debug
@ -782,10 +782,10 @@ matID = constitutive_matID(ipc,ip,el)
!* Calculation of Lp !* Calculation of Lp
Lp = 0.0_pReal Lp = 0.0_pReal
do i=1,material_Nslip(matID) do i=1,material_Nslip(matID)
tau_slip(i)=dot_product(Tstar_v,crystal_Sslip_v(:,i,material_CrystalStructure(matID))) tau_slip(i)=dot_product(Tstar_v,lattice_Sslip_v(:,i,material_CrystalStructure(matID)))
gdot_slip(i)=material_gdot0_slip(matID)*(abs(tau_slip(i))/state(i))**& gdot_slip(i)=material_gdot0_slip(matID)*(abs(tau_slip(i))/state(i))**&
material_n_slip(matID)*sign(1.0_pReal,tau_slip(i)) material_n_slip(matID)*sign(1.0_pReal,tau_slip(i))
Lp=Lp+gdot_slip(i)*crystal_Sslip(:,:,i,material_CrystalStructure(matID)) Lp=Lp+gdot_slip(i)*lattice_Sslip(:,:,i,material_CrystalStructure(matID))
enddo enddo
!* Calculation of the tangent of Lp !* Calculation of the tangent of Lp
@ -796,8 +796,8 @@ do i=1,material_Nslip(matID)
(material_n_slip(matID)-1.0_pReal)*material_n_slip(matID)/state(i) (material_n_slip(matID)-1.0_pReal)*material_n_slip(matID)/state(i)
forall (k=1:3,l=1:3,m=1:3,n=1:3) & forall (k=1:3,l=1:3,m=1:3,n=1:3) &
dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + & dLp_dTstar3333(k,l,m,n) = dLp_dTstar3333(k,l,m,n) + &
dgdot_dtauslip(i)*crystal_Sslip(k,l,i,material_CrystalStructure(matID))* & dgdot_dtauslip(i)*lattice_Sslip(k,l,i,material_CrystalStructure(matID))* &
crystal_Sslip(m,n,i,material_CrystalStructure(matID)) lattice_Sslip(m,n,i,material_CrystalStructure(matID))
enddo enddo
dLp_dTstar = math_Plain3333to99(dLp_dTstar3333) dLp_dTstar = math_Plain3333to99(dLp_dTstar3333)
@ -819,7 +819,7 @@ function constitutive_dotState(Tstar_v,state,Temperature,ipc,ip,el)
!* - constitutive_dotState : evolution of state variable * !* - constitutive_dotState : evolution of state variable *
!********************************************************************* !*********************************************************************
use prec, only: pReal,pInt use prec, only: pReal,pInt
use crystal, only: crystal_Sslip_v use lattice, only: lattice_Sslip_v
implicit none implicit none
!* Definition of variables !* Definition of variables
@ -834,7 +834,7 @@ matID = constitutive_matID(ipc,ip,el)
!* Self-Hardening of each system !* Self-Hardening of each system
do i=1,constitutive_Nstatevars(ipc,ip,el) do i=1,constitutive_Nstatevars(ipc,ip,el)
tau_slip = dot_product(Tstar_v,crystal_Sslip_v(:,i,material_CrystalStructure(matID))) tau_slip = dot_product(Tstar_v,lattice_Sslip_v(:,i,material_CrystalStructure(matID)))
gdot_slip = material_gdot0_slip(matID)*(abs(tau_slip)/state(i))**& gdot_slip = material_gdot0_slip(matID)*(abs(tau_slip)/state(i))**&
material_n_slip(matID)*sign(1.0_pReal,tau_slip) material_n_slip(matID)*sign(1.0_pReal,tau_slip)
self_hardening(i)=material_h0(matID)*(1.0_pReal-state(i)/& self_hardening(i)=material_h0(matID)*(1.0_pReal-state(i)/&
@ -848,7 +848,7 @@ return
end function end function
function constitutive_post_results(Tstar_v,state,dt,Temperature,ipc,ip,el) function constitutive_post_results(Tstar_v,state,Temperature,dt,ipc,ip,el)
!********************************************************************* !*********************************************************************
!* return array of constitutive results * !* return array of constitutive results *
!* INPUT: * !* INPUT: *
@ -860,7 +860,7 @@ function constitutive_post_results(Tstar_v,state,dt,Temperature,ipc,ip,el)
!* - el : current element * !* - el : current element *
!********************************************************************* !*********************************************************************
use prec, only: pReal,pInt use prec, only: pReal,pInt
use crystal, only: crystal_Sslip_v use lattice, only: lattice_Sslip_v
implicit none implicit none
!* Definition of variables !* Definition of variables
@ -878,7 +878,7 @@ if(constitutive_Nresults(ipc,ip,el)==0) return
do i=1,material_Nslip(matID) do i=1,material_Nslip(matID)
constitutive_post_results(i) = state(i) constitutive_post_results(i) = state(i)
tau_slip=dot_product(Tstar_v,crystal_Sslip_v(:,i,material_CrystalStructure(matID))) tau_slip=dot_product(Tstar_v,lattice_Sslip_v(:,i,material_CrystalStructure(matID)))
constitutive_post_results(i+material_Nslip(matID)) = & constitutive_post_results(i+material_Nslip(matID)) = &
dt*material_gdot0_slip(matID)*(abs(tau_slip)/state(i))**material_n_slip(matID)*sign(1.0_pReal,tau_slip) dt*material_gdot0_slip(matID)*(abs(tau_slip)/state(i))**material_n_slip(matID)*sign(1.0_pReal,tau_slip)
enddo enddo

307
trunk/crystallite.f90 Normal file
View File

@ -0,0 +1,307 @@
!##############################################################
MODULE crystallite
!##############################################################
! *** Solution at single crystallite level ***
!
CONTAINS
!
!
!********************************************************************
! Calculates the stress for a single component
!********************************************************************
!***********************************************************************
!*** calculation of stress (P), stiffness (dPdF), ***
!*** and announcment of any ***
!*** acceleration of the Newton-Raphson correction ***
!***********************************************************************
subroutine SingleCrystallite(&
msg,& ! return message
P,& ! first PK stress
dPdF,& ! consistent tangent
post_results,& ! plot results from constitutive model
Fp_new,& ! new plastic deformation gradient
Fe_new,& ! new "elastic" deformation gradient
state_new,& ! new state variable array
!
dt,& ! time increment
cp_en,& ! element number
ip,& ! integration point number
grain,& ! grain number
updateJaco,& ! update of Jacobian required
Temperature,& ! temperature of crystallite
Fg_new,& ! new global deformation gradient
Fg_old,& ! old global deformation gradient
Fp_old,& ! old plastic deformation gradient
state_old) ! old state variable array
!
use prec, only: pReal,pInt,pert_Fg,subStepMin
use debug
use constitutive, only: constitutive_Nstatevars,constitutive_Nresults
use mesh, only: mesh_element
use math, only: math_Mandel6to33,math_Mandel33to6,math_Mandel3333to66,&
math_I3,math_det3x3,math_invert3x3
implicit none
!
character(len=*) msg
logical updateJaco,error,cuttedBack,guessNew
integer(pInt) cp_en,ip,grain,i,j,k,l,m,n, nCutbacks
real(pReal) Temperature
real(pReal) dt,dt_aim,subFrac,subStep,invJ,det
real(pReal), dimension(3,3) :: Lp,Lp_pert,inv
real(pReal), dimension(3,3) :: Fg_old,Fg_current,Fg_aim,Fg_new,Fg_pert,deltaFg
real(pReal), dimension(3,3) :: Fp_old,Fp_current,Fp_new,Fp_pert
real(pReal), dimension(3,3) :: Fe_old,Fe_current,Fe_new,Fe_pert
real(pReal), dimension(3,3) :: Tstar,tau,P,P_pert
real(pReal), dimension(3,3,3,3) :: dPdF
real(pReal), dimension(constitutive_Nstatevars(grain,ip,cp_en)) :: state_old,state_current,state_new,state_pert
real(pReal), dimension(constitutive_Nresults(grain,ip,cp_en)) :: post_results
!
deltaFg = Fg_new-Fg_old
subFrac = 0.0_pReal
subStep = 1.0_pReal
nCutbacks = 0_pInt
!
Fg_aim = Fg_old ! make "new", "aim" a synonym for "old"
Fp_new = Fp_old
call math_invert3x3(Fp_old,inv,det,error)
Fe_new = matmul(Fg_old,inv)
state_new = state_old
!
cuttedBack = .false.
guessNew = .true.
!
! begin the cutback loop
do while (subStep > subStepMin) ! continue until finished or too much cut backing
if (.not. cuttedBack) then
Fg_current = Fg_aim ! wind forward
Fp_current = Fp_new
Fe_current = Fe_new
state_current = state_new
endif
!
Fg_aim = Fg_current + subStep*deltaFg ! aim for Fg
dt_aim = subStep*dt ! aim for dt
msg = '' ! error free so far
if (guessNew) then ! calculate new Lp guess when cutted back
if (dt_aim /= 0.0_pReal) then
call math_invert3x3(Fg_aim,inv,det,error)
Lp = (math_I3-matmul(Fp_current,matmul(inv,Fe_current)))/dt ! fully plastic initial guess
else
Lp = 0.0_pReal ! fully elastic guess
endif
endif
call TimeIntegration(msg,Lp,Fp_new,Fe_new,P,state_new,post_results,.true., & ! def gradients and PK2 at end of time step
dt_aim,cp_en,ip,grain,Temperature,Fg_aim,Fp_current,state_current)
!
if (msg == 'ok') then
cuttedBack = .false. ! no cut back required
guessNew = .false. ! keep the Lp
subFrac = subFrac + subStep
subStep = 1.0_pReal - subFrac ! try one go
else
nCutbacks = nCutbacks + 1 ! record additional cutback
cuttedBack = .true. ! encountered problems -->
guessNew = .true. ! redo plastic Lp guess
subStep = subStep / 2.0_pReal ! cut time step in half
endif
enddo
!
debug_cutbackDistribution(nCutbacks+1) = debug_cutbackDistribution(nCutbacks+1)+1
!
if (msg /= 'ok') return ! solution not reached --> report back
if (updateJaco) then ! consistent tangent using
do k=1,3
do l=1,3
Fg_pert = Fg_new ! initialize perturbed Fg
Fg_pert(k,l) = Fg_pert(k,l) + pert_Fg ! perturb single component
Lp_pert = Lp
state_pert = state_new ! initial guess from end of time step
call TimeIntegration(msg,Lp,Fp_pert,Fe_pert,P_pert,state_pert,post_results,.false., & ! def gradients and PK2 at end of time step
dt_aim,cp_en,ip,grain,Temperature,Fg_pert,Fp_current,state_current)
if (msg /= 'ok') then
msg = 'consistent tangent --> '//msg
return
endif
dPdF(:,:,k,l) = (P_pert-P)/pert_Fg ! constructing the tangent dP_ij/dFg_kl from forward differences
enddo
enddo
endif
!
return
!
END SUBROUTINE
!
!***********************************************************************
!*** fully-implicit two-level time integration ***
!*** based on a residuum in Lp and intermediate ***
!*** acceleration of the Newton-Raphson correction ***
!***********************************************************************
SUBROUTINE TimeIntegration(&
msg,& ! return message
Lpguess,& ! guess of plastic velocity gradient
Fp_new,& ! new plastic deformation gradient
Fe_new,& ! new "elastic" deformation gradient
P,& ! 1nd PK stress (taken as initial guess if /= 0)
state,& ! current microstructure at end of time inc (taken as guess if /= 0)
results,& ! post results from constitutive
wantsConstitutiveResults,& ! its flag
!
dt,& ! time increment
cp_en,& ! element number
ip,& ! integration point number
grain,& ! grain number
Temperature,& ! temperature
Fg_new,& ! new total def gradient
Fp_old,& ! former plastic def gradient
state_old) ! former microstructure
!
use prec
use debug
use mesh, only: mesh_element
use constitutive, only: constitutive_Nstatevars,&
constitutive_homogenizedC,constitutive_dotState,constitutive_LpAndItsTangent,&
constitutive_Nresults,constitutive_Microstructure,constitutive_post_results
use math
implicit none
!
character(len=*) msg
logical failed,wantsConstitutiveResults
integer(pInt) cp_en, ip, grain
integer(pInt) iOuter,iInner,dummy, i,j,k,l,m,n
real(pReal) dt, Temperature, det, p_hydro, leapfrog,maxleap
real(pReal), dimension(6) :: Tstar_v
real(pReal), dimension(9,9) :: dLp,dTdLp,dRdLp,invdRdLp,eye2
real(pReal), dimension(6,6) :: C_66
real(pReal), dimension(3,3) :: Fg_new,Fp_new,invFp_new,Fp_old,invFp_old,Fe_new,Fe_old
real(pReal), dimension(3,3) :: P,Tstar
real(pReal), dimension(3,3) :: Lp,Lpguess,Lpguess_old,Rinner,Rinner_old,A,B,BT,AB,BTA
real(pReal), dimension(3,3,3,3) :: C
real(pReal), dimension(constitutive_Nstatevars(grain, ip, cp_en)) :: state_old,state,ROuter
real(pReal), dimension(constitutive_Nresults(grain,ip,cp_en)) :: results
!
msg = 'ok' ! error-free so far
eye2 = math_identity2nd(9)
call math_invert3x3(Fp_old,invFp_old,det,failed) ! inversion of Fp_old
if (failed) then
msg = 'inversion Fp_old'
return
endif
A = matmul(transpose(invFp_old), matmul(transpose(Fg_new),matmul(Fg_new,invFp_old)))
!
if (all(state == 0.0_pReal)) state = state_old ! former state guessed, if none specified
iOuter = 0_pInt ! outer counter
!
!
Outer: do ! outer iteration: State
iOuter = iOuter+1
if (iOuter > nOuter) then
msg = 'limit Outer iteration'
debug_OuterLoopDistribution(nOuter) = debug_OuterLoopDistribution(nOuter)+1
return
endif
call constitutive_Microstructure(state,Temperature,grain,ip,cp_en)
C_66 = constitutive_HomogenizedC(state, grain, ip, cp_en)
C = math_Mandel66to3333(C_66) ! 4th rank elasticity tensor
!
iInner = 0_pInt
leapfrog = 1.0_pReal ! correction as suggested by invdRdLp-step
maxleap = 1024.0_pReal ! preassign maximum acceleration level
!
Inner: do ! inner iteration: Lp
iInner = iInner+1
if (iInner > nInner) then ! too many loops required
msg = 'limit Inner iteration'
debug_InnerLoopDistribution(nInner) = debug_InnerLoopDistribution(nInner)+1
return
endif
B = math_i3 - dt*Lpguess
BT = transpose(B)
AB = matmul(A,B)
BTA = matmul(BT,A)
Tstar_v = 0.5_pReal*matmul(C_66,math_mandel33to6(matmul(BT,AB)-math_I3))
Tstar = math_Mandel6to33(Tstar_v)
p_hydro=(Tstar_v(1)+Tstar_v(2)+Tstar_v(3))/3.0_pReal
forall(i=1:3) Tstar_v(i) = Tstar_v(i)-p_hydro ! subtract hydrostatic pressure
call constitutive_LpAndItsTangent(Lp,dLp, &
Tstar_v,state,Temperature,grain,ip,cp_en)
Rinner = Lpguess - Lp ! update current residuum
if ((maxval(abs(Rinner)) < abstol_Inner) .or. &
(any(abs(dt*Lpguess) > relevantStrain) .and. &
maxval(abs(Rinner/Lpguess),abs(dt*Lpguess) > relevantStrain) < reltol_Inner)) &
exit Inner
!
! check for acceleration/deceleration in Newton--Raphson correction
if (leapfrog > 1.0_pReal .and. &
(sum(Rinner*Rinner) > sum(Rinner_old*Rinner_old) .or. & ! worse residuum
sum(Rinner*Rinner_old) < 0.0_pReal)) then ! residuum changed sign (overshoot)
maxleap = 0.5_pReal * leapfrog ! limit next acceleration
leapfrog = 1.0_pReal ! grinding halt
else ! better residuum
dTdLp = 0.0_pReal ! calc dT/dLp
forall (i=1:3,j=1:3,k=1:3,l=1:3,m=1:3,n=1:3) &
dTdLp(3*(i-1)+j,3*(k-1)+l) = dTdLp(3*(i-1)+j,3*(k-1)+l) + &
C(i,j,l,n)*AB(k,n)+C(i,j,m,l)*BTA(m,k)
dTdLp = -0.5_pReal*dt*dTdLp
dRdLp = eye2 - matmul(dLp,dTdLp) ! calc dR/dLp
invdRdLp = 0.0_pReal
call math_invert(9,dRdLp,invdRdLp,dummy,failed) ! invert dR/dLp --> dLp/dR
if (failed) then
msg = 'inversion dR/dLp'
if (debugger) then
write (6,*) msg
write (6,'(a,/,9(9(e9.3,x)/))') 'dRdLp', dRdLp(1:9,:)
write (6,*) 'state',state
write (6,'(a,/,3(3(f12.7,x)/))') 'Lpguess',Lpguess(1:3,:)
write (6,*) 'Tstar',Tstar_v
endif
return
endif
!
Rinner_old = Rinner ! remember current residuum
Lpguess_old = Lpguess ! remember current Lp guess
if (iInner > 1 .and. leapfrog < maxleap) leapfrog = 2.0_pReal * leapfrog ! accelerate
endif
!
Lpguess = Lpguess_old ! start from current guess
Rinner = Rinner_old ! use current residuum
forall (i=1:3,j=1:3,k=1:3,l=1:3) & ! leapfrog to updated Lpguess
Lpguess(i,j) = Lpguess(i,j) - leapfrog*invdRdLp(3*(i-1)+j,3*(k-1)+l)*Rinner(k,l)
enddo Inner
!
debug_InnerLoopDistribution(iInner) = debug_InnerLoopDistribution(iInner)+1
ROuter = state - state_old - &
dt*constitutive_dotState(Tstar_v,state,Temperature,&
grain,ip,cp_en) ! residuum from evolution of microstructure
state = state - ROuter ! update of microstructure
if (maxval(abs(Router/state),state /= 0.0_pReal) < reltol_Outer) exit Outer
enddo Outer
!
debug_OuterLoopDistribution(iOuter) = debug_OuterLoopDistribution(iOuter)+1
invFp_new = matmul(invFp_old,B)
call math_invert3x3(invFp_new,Fp_new,det,failed)
if (failed) then
msg = 'inversion Fp_new^-1'
return
endif
!
if (wantsConstitutiveResults) then ! get the post_results upon request
results = 0.0_pReal
results = constitutive_post_results(Tstar_v,state,Temperature,dt,grain,ip,cp_en)
endif
!
Fp_new = Fp_new*det**(1.0_pReal/3.0_pReal) ! regularize Fp by det = det(InvFp_new) !!
Fe_new = matmul(Fg_new,invFp_new) ! calc resulting Fe
forall (i=1:3) Tstar_v(i) = Tstar_v(i)+p_hydro ! add hydrostatic component back
P = matmul(Fe_new,matmul(Tstar,transpose(invFp_new))) ! first PK stress
!
return
!
END SUBROUTINE
!!
END MODULE
!##############################################################

370
trunk/lattice.f90 Normal file
View File

@ -0,0 +1,370 @@
!************************************
!* Module: LATTICE *
!************************************
!* contains: *
!* - Lattice structure definition *
!* - Slip system definition *
!* - Schmid matrices calculation *
!************************************
MODULE lattice
!*** Include other modules ***
use prec, only: pReal,pInt
implicit none
!************************************
!* Lattice structures *
!************************************
!* Number of lattice structures (1-FCC,2-BCC,3-HCP)
integer(pInt), parameter :: lattice_MaxLatticeStructure = 3
!* Total number of slip systems per lattice structure
!* (has to be changed according the definition of slip systems)
integer(pInt), dimension(lattice_MaxLatticeStructure), parameter :: lattice_MaxNslipOfStructure = &
reshape((/12,48,12/),(/lattice_MaxLatticeStructure/))
!* Total number of twin systems per lattice structure
!* (has to be changed according the definition of twin systems)
integer(pInt), dimension(lattice_MaxLatticeStructure), parameter :: lattice_MaxNtwinOfStructure = &
reshape((/12,12,6/),(/lattice_MaxLatticeStructure/))
!* Maximum number of slip systems over lattice structures
integer(pInt), parameter :: lattice_MaxMaxNslipOfStructure = 48
!* Maximum number of twin systems over lattice structures
integer(pInt), parameter :: lattice_MaxMaxNtwinOfStructure = 12
!* Slip direction, slip normales and Schmid matrices
real(pReal), dimension(3,3,lattice_MaxMaxNslipOfStructure,lattice_MaxLatticeStructure) :: lattice_Sslip
real(pReal), dimension(6,lattice_MaxMaxNslipOfStructure,lattice_MaxLatticeStructure) :: lattice_Sslip_v
real(pReal), dimension(3,lattice_MaxMaxNslipOfStructure,lattice_MaxLatticeStructure) :: lattice_sn
real(pReal), dimension(3,lattice_MaxMaxNslipOfStructure,lattice_MaxLatticeStructure) :: lattice_sd
real(pReal), dimension(3,lattice_MaxMaxNslipOfStructure,lattice_MaxLatticeStructure) :: lattice_st
!* twin direction, twin normales, Schmid matrices and transformation matrices
real(pReal), dimension(3,3,lattice_MaxMaxNtwinOfStructure,lattice_MaxLatticeStructure) :: lattice_Stwin
real(pReal), dimension(6,lattice_MaxMaxNtwinOfStructure,lattice_MaxLatticeStructure) :: lattice_Stwin_v
real(pReal), dimension(3,lattice_MaxMaxNtwinOfStructure,lattice_MaxLatticeStructure) :: lattice_tn
real(pReal), dimension(3,lattice_MaxMaxNtwinOfStructure,lattice_MaxLatticeStructure) :: lattice_td
real(pReal), dimension(3,lattice_MaxMaxNtwinOfStructure,lattice_MaxLatticeStructure) :: lattice_tt
real(pReal), dimension(3,3,lattice_MaxMaxNtwinOfStructure,lattice_MaxLatticeStructure) :: lattice_Qtwin
real(pReal), dimension(lattice_MaxLatticeStructure), parameter :: lattice_TwinShear = &
reshape((/0.7071067812,0.7071067812,0.7071067812/),(/lattice_MaxLatticeStructure/)) ! Depends surely on c/a ratio for HCP
!* Slip_slip interaction matrices
integer(pInt), dimension(lattice_MaxMaxNslipOfStructure,lattice_MaxMaxNslipOfStructure,lattice_MaxLatticeStructure) :: &
lattice_SlipIntType
!* Twin-twin interaction matrices
integer(pInt), dimension(lattice_MaxMaxNtwinOfStructure,lattice_MaxMaxNtwinOfStructure,lattice_MaxLatticeStructure) :: &
lattice_TwinIntType
!*** Slip systems for FCC structures (1) ***
!* System {111}<110> Sort according Eisenlohr&Hantcherli
data lattice_sd(:, 1,1)/ 0, 1,-1/ ; data lattice_sn(:, 1,1)/ 1, 1, 1/
data lattice_sd(:, 2,1)/-1, 0, 1/ ; data lattice_sn(:, 2,1)/ 1, 1, 1/
data lattice_sd(:, 3,1)/ 1,-1, 0/ ; data lattice_sn(:, 3,1)/ 1, 1, 1/
data lattice_sd(:, 4,1)/ 0,-1,-1/ ; data lattice_sn(:, 4,1)/-1,-1, 1/
data lattice_sd(:, 5,1)/ 1, 0, 1/ ; data lattice_sn(:, 5,1)/-1,-1, 1/
data lattice_sd(:, 6,1)/-1, 1, 0/ ; data lattice_sn(:, 6,1)/-1,-1, 1/
data lattice_sd(:, 7,1)/ 0,-1, 1/ ; data lattice_sn(:, 7,1)/ 1,-1,-1/
data lattice_sd(:, 8,1)/-1, 0,-1/ ; data lattice_sn(:, 8,1)/ 1,-1,-1/
data lattice_sd(:, 9,1)/ 1, 1, 0/ ; data lattice_sn(:, 9,1)/ 1,-1,-1/
data lattice_sd(:,10,1)/ 0, 1, 1/ ; data lattice_sn(:,10,1)/-1, 1,-1/
data lattice_sd(:,11,1)/ 1, 0,-1/ ; data lattice_sn(:,11,1)/-1, 1,-1/
data lattice_sd(:,12,1)/-1,-1, 0/ ; data lattice_sn(:,12,1)/-1, 1,-1/
!*** Twin systems for FCC structures (1) ***
!* System {111}<112> Sort according Eisenlohr&Hantcherli
data lattice_td(:, 1,1)/-2, 1, 1/ ; data lattice_tn(:, 1,1)/ 1, 1, 1/
data lattice_td(:, 2,1)/ 1,-2, 1/ ; data lattice_tn(:, 2,1)/ 1, 1, 1/
data lattice_td(:, 3,1)/ 1, 1,-2/ ; data lattice_tn(:, 3,1)/ 1, 1, 1/
data lattice_td(:, 4,1)/ 2,-1, 1/ ; data lattice_tn(:, 4,1)/-1,-1, 1/
data lattice_td(:, 5,1)/-1, 2, 1/ ; data lattice_tn(:, 5,1)/-1,-1, 1/
data lattice_td(:, 6,1)/-1,-1,-2/ ; data lattice_tn(:, 6,1)/-1,-1, 1/
data lattice_td(:, 7,1)/-2,-1,-1/ ; data lattice_tn(:, 7,1)/ 1,-1,-1/
data lattice_td(:, 8,1)/ 1, 2,-1/ ; data lattice_tn(:, 8,1)/ 1,-1,-1/
data lattice_td(:, 9,1)/ 1,-1, 2/ ; data lattice_tn(:, 9,1)/ 1,-1,-1/
data lattice_td(:,10,1)/ 2, 1,-1/ ; data lattice_tn(:,10,1)/-1, 1,-1/
data lattice_td(:,11,1)/-1,-2,-1/ ; data lattice_tn(:,11,1)/-1, 1,-1/
data lattice_td(:,12,1)/-1, 1, 2/ ; data lattice_tn(:,12,1)/-1, 1,-1/
!*** Slip-Slip interactions for FCC structures (1) ***
data lattice_SlipIntType( 1,1:lattice_MaxNslipOfStructure(1),1)/1,2,2,4,6,5,3,5,5,4,5,6/
data lattice_SlipIntType( 2,1:lattice_MaxNslipOfStructure(1),1)/2,1,2,6,4,5,5,4,6,5,3,5/
data lattice_SlipIntType( 3,1:lattice_MaxNslipOfStructure(1),1)/2,2,1,5,5,3,5,6,4,6,5,4/
data lattice_SlipIntType( 4,1:lattice_MaxNslipOfStructure(1),1)/4,6,5,1,2,2,4,5,6,3,5,5/
data lattice_SlipIntType( 5,1:lattice_MaxNslipOfStructure(1),1)/6,4,5,2,1,2,5,3,5,5,4,6/
data lattice_SlipIntType( 6,1:lattice_MaxNslipOfStructure(1),1)/5,5,3,2,2,1,6,5,4,5,6,4/
data lattice_SlipIntType( 7,1:lattice_MaxNslipOfStructure(1),1)/3,5,5,4,5,6,1,2,2,4,6,5/
data lattice_SlipIntType( 8,1:lattice_MaxNslipOfStructure(1),1)/5,4,6,5,3,5,2,1,2,6,4,5/
data lattice_SlipIntType( 9,1:lattice_MaxNslipOfStructure(1),1)/5,6,4,6,5,4,2,2,1,5,5,3/
data lattice_SlipIntType(10,1:lattice_MaxNslipOfStructure(1),1)/4,5,6,3,5,5,4,6,5,1,2,2/
data lattice_SlipIntType(11,1:lattice_MaxNslipOfStructure(1),1)/5,3,5,5,4,6,6,4,5,2,1,2/
data lattice_SlipIntType(12,1:lattice_MaxNslipOfStructure(1),1)/6,5,4,5,6,4,5,5,3,2,2,1/
!*** Twin-Twin interactions for FCC structures (1) ***
data lattice_TwinIntType( 1,1:lattice_MaxNtwinOfStructure(1),1)/0,0,0,1,1,1,1,1,1,1,1,1/
data lattice_TwinIntType( 2,1:lattice_MaxNtwinOfStructure(1),1)/0,0,0,1,1,1,1,1,1,1,1,1/
data lattice_TwinIntType( 3,1:lattice_MaxNtwinOfStructure(1),1)/0,0,0,1,1,1,1,1,1,1,1,1/
data lattice_TwinIntType( 4,1:lattice_MaxNtwinOfStructure(1),1)/1,1,1,0,0,0,1,1,1,1,1,1/
data lattice_TwinIntType( 5,1:lattice_MaxNtwinOfStructure(1),1)/1,1,1,0,0,0,1,1,1,1,1,1/
data lattice_TwinIntType( 6,1:lattice_MaxNtwinOfStructure(1),1)/1,1,1,0,0,0,1,1,1,1,1,1/
data lattice_TwinIntType( 7,1:lattice_MaxNtwinOfStructure(1),1)/1,1,1,1,1,1,0,0,0,1,1,1/
data lattice_TwinIntType( 8,1:lattice_MaxNtwinOfStructure(1),1)/1,1,1,1,1,1,0,0,0,1,1,1/
data lattice_TwinIntType( 9,1:lattice_MaxNtwinOfStructure(1),1)/1,1,1,1,1,1,0,0,0,1,1,1/
data lattice_TwinIntType(10,1:lattice_MaxNtwinOfStructure(1),1)/1,1,1,1,1,1,1,1,1,0,0,0/
data lattice_TwinIntType(11,1:lattice_MaxNtwinOfStructure(1),1)/1,1,1,1,1,1,1,1,1,0,0,0/
data lattice_TwinIntType(12,1:lattice_MaxNtwinOfStructure(1),1)/1,1,1,1,1,1,1,1,1,0,0,0/
!*** Slip systems for BCC structures (2) ***
!* System {110}<111>
!* Sort?
data lattice_sd(:, 1,2)/ 1,-1, 1/ ; data lattice_sn(:, 1,2)/ 0, 1, 1/
data lattice_sd(:, 2,2)/-1,-1, 1/ ; data lattice_sn(:, 2,2)/ 0, 1, 1/
data lattice_sd(:, 3,2)/ 1, 1, 1/ ; data lattice_sn(:, 3,2)/ 0,-1, 1/
data lattice_sd(:, 4,2)/-1, 1, 1/ ; data lattice_sn(:, 4,2)/ 0,-1, 1/
data lattice_sd(:, 5,2)/-1, 1, 1/ ; data lattice_sn(:, 5,2)/ 1, 0, 1/
data lattice_sd(:, 6,2)/-1,-1, 1/ ; data lattice_sn(:, 6,2)/ 1, 0, 1/
data lattice_sd(:, 7,2)/ 1, 1, 1/ ; data lattice_sn(:, 7,2)/-1, 0, 1/
data lattice_sd(:, 8,2)/ 1,-1, 1/ ; data lattice_sn(:, 8,2)/-1, 0, 1/
data lattice_sd(:, 9,2)/-1, 1, 1/ ; data lattice_sn(:, 9,2)/ 1, 1, 0/
data lattice_sd(:,10,2)/-1, 1,-1/ ; data lattice_sn(:,10,2)/ 1, 1, 0/
data lattice_sd(:,11,2)/ 1, 1, 1/ ; data lattice_sn(:,11,2)/-1, 1, 0/
data lattice_sd(:,12,2)/ 1, 1,-1/ ; data lattice_sn(:,12,2)/-1, 1, 0/
!* System {112}<111>
!* Sort?
data lattice_sd(:,13,2)/-1, 1, 1/ ; data lattice_sn(:,13,2)/ 2, 1, 1/
data lattice_sd(:,14,2)/ 1, 1, 1/ ; data lattice_sn(:,14,2)/-2, 1, 1/
data lattice_sd(:,15,2)/ 1, 1,-1/ ; data lattice_sn(:,15,2)/ 2,-1, 1/
data lattice_sd(:,16,2)/ 1,-1, 1/ ; data lattice_sn(:,16,2)/ 2, 1,-1/
data lattice_sd(:,17,2)/ 1,-1, 1/ ; data lattice_sn(:,17,2)/ 1, 2, 1/
data lattice_sd(:,18,2)/ 1, 1,-1/ ; data lattice_sn(:,18,2)/-1, 2, 1/
data lattice_sd(:,19,2)/ 1, 1, 1/ ; data lattice_sn(:,19,2)/ 1,-2, 1/
data lattice_sd(:,20,2)/-1, 1, 1/ ; data lattice_sn(:,20,2)/ 1, 2,-1/
data lattice_sd(:,21,2)/ 1, 1,-1/ ; data lattice_sn(:,21,2)/ 1, 1, 2/
data lattice_sd(:,22,2)/ 1,-1, 1/ ; data lattice_sn(:,22,2)/-1, 1, 2/
data lattice_sd(:,23,2)/-1, 1, 1/ ; data lattice_sn(:,23,2)/ 1,-1, 2/
data lattice_sd(:,24,2)/ 1, 1, 1/ ; data lattice_sn(:,24,2)/ 1, 1,-2/
!* System {123}<111>
!* Sort?
data lattice_sd(:,25,2)/ 1, 1,-1/ ; data lattice_sn(:,25,2)/ 1, 2, 3/
data lattice_sd(:,26,2)/ 1,-1, 1/ ; data lattice_sn(:,26,2)/-1, 2, 3/
data lattice_sd(:,27,2)/-1, 1, 1/ ; data lattice_sn(:,27,2)/ 1,-2, 3/
data lattice_sd(:,28,2)/ 1, 1, 1/ ; data lattice_sn(:,28,2)/ 1, 2,-3/
data lattice_sd(:,29,2)/ 1,-1, 1/ ; data lattice_sn(:,29,2)/ 1, 3, 2/
data lattice_sd(:,30,2)/ 1, 1,-1/ ; data lattice_sn(:,30,2)/-1, 3, 2/
data lattice_sd(:,31,2)/ 1, 1, 1/ ; data lattice_sn(:,31,2)/ 1,-3, 2/
data lattice_sd(:,32,2)/-1, 1, 1/ ; data lattice_sn(:,32,2)/ 1, 3,-2/
data lattice_sd(:,33,2)/ 1, 1,-1/ ; data lattice_sn(:,33,2)/ 2, 1, 3/
data lattice_sd(:,34,2)/ 1,-1, 1/ ; data lattice_sn(:,34,2)/-2, 1, 3/
data lattice_sd(:,35,2)/-1, 1, 1/ ; data lattice_sn(:,35,2)/ 2,-1, 3/
data lattice_sd(:,36,2)/ 1, 1, 1/ ; data lattice_sn(:,36,2)/ 2, 1,-3/
data lattice_sd(:,37,2)/ 1,-1, 1/ ; data lattice_sn(:,37,2)/ 2, 3, 1/
data lattice_sd(:,38,2)/ 1, 1,-1/ ; data lattice_sn(:,38,2)/-2, 3, 1/
data lattice_sd(:,39,2)/ 1, 1, 1/ ; data lattice_sn(:,39,2)/ 2,-3, 1/
data lattice_sd(:,40,2)/-1, 1, 1/ ; data lattice_sn(:,40,2)/ 2, 3,-1/
data lattice_sd(:,41,2)/-1, 1, 1/ ; data lattice_sn(:,41,2)/ 3, 1, 2/
data lattice_sd(:,42,2)/ 1, 1, 1/ ; data lattice_sn(:,42,2)/-3, 1, 2/
data lattice_sd(:,43,2)/ 1, 1,-1/ ; data lattice_sn(:,43,2)/ 3,-1, 2/
data lattice_sd(:,44,2)/ 1,-1, 1/ ; data lattice_sn(:,44,2)/ 3, 1,-2/
data lattice_sd(:,45,2)/-1, 1, 1/ ; data lattice_sn(:,45,2)/ 3, 2, 1/
data lattice_sd(:,46,2)/ 1, 1, 1/ ; data lattice_sn(:,46,2)/-3, 2, 1/
data lattice_sd(:,47,2)/ 1, 1,-1/ ; data lattice_sn(:,47,2)/ 3,-2, 1/
data lattice_sd(:,48,2)/ 1,-1, 1/ ; data lattice_sn(:,48,2)/ 3, 2,-1/
!*** Twin systems for BCC structures (2) ***
!* System {112}<111>
!* Sort?
!* MISSING: not implemented yet
!*** Slip-Slip interactions for BCC structures (2) ***
data lattice_SlipIntType( 1,:,2)/1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType( 2,:,2)/2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType( 3,:,2)/2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType( 4,:,2)/2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType( 5,:,2)/2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType( 6,:,2)/2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType( 7,:,2)/2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType( 8,:,2)/2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType( 9,:,2)/2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType(10,:,2)/2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType(11,:,2)/2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType(12,:,2)/2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType(13,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType(14,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType(15,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType(16,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType(17,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType(18,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType(19,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType(20,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType(21,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType(22,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType(23,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType(24,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType(25,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType(26,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType(27,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType(28,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType(29,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType(30,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType(31,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType(32,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType(33,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType(34,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType(35,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType(36,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType(37,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType(38,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType(39,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType(40,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2,2/
data lattice_SlipIntType(41,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2,2/
data lattice_SlipIntType(42,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2,2/
data lattice_SlipIntType(43,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2,2/
data lattice_SlipIntType(44,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2,2/
data lattice_SlipIntType(45,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2,2/
data lattice_SlipIntType(46,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2,2/
data lattice_SlipIntType(47,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1,2/
data lattice_SlipIntType(48,:,2)/2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,2,1/
!*** Twin-twin interactions for BCC structures (2) ***
! MISSING: not implemented yet
!*** Slip systems for HCP structures (3) ***
!* Basal systems {0001}<1120> (independent of c/a-ratio)
!* 1- (0 0 0 1)[-2 1 1 0]
!* 2- (0 0 0 1)[ 1 -2 1 0]
!* 3- (0 0 0 1)[ 1 1 -2 0]
!* Plane (hkil)->(hkl)
!* Direction [uvtw]->[(u-t) (v-t) w]
!* Automatical transformation from Bravais to Miller
!* not done for the moment
!* Sort?
data lattice_sd(:, 1,3)/-1, 0, 0/ ; data lattice_sn(:, 1,3)/ 0, 0, 1/
data lattice_sd(:, 2,3)/ 0,-1, 0/ ; data lattice_sn(:, 2,3)/ 0, 0, 1/
data lattice_sd(:, 3,3)/ 1, 1, 0/ ; data lattice_sn(:, 3,3)/ 0, 0, 1/
!* 1st type prismatic systems {1010}<1120> (independent of c/a-ratio)
!* 1- ( 0 1 -1 0)[-2 1 1 0]
!* 2- ( 1 0 -1 0)[ 1 -2 1 0]
!* 3- (-1 1 0 0)[ 1 1 -2 0]
!* Sort?
data lattice_sd(:, 4,3)/-1, 0, 0/ ; data lattice_sn(:, 4,3)/ 0, 1, 0/
data lattice_sd(:, 5,3)/ 0,-1, 0/ ; data lattice_sn(:, 5,3)/ 1, 0, 0/
data lattice_sd(:, 6,3)/ 1, 1, 0/ ; data lattice_sn(:, 6,3)/-1, 1, 0/
!* 1st type 1st order pyramidal systems {1011}<1120>
!* plane normales depend on the c/a-ratio
!* 1- ( 0 -1 1 1)[-2 1 1 0]
!* 2- ( 0 1 -1 1)[-2 1 1 0]
!* 3- (-1 0 1 1)[ 1 -2 1 0]
!* 4- ( 1 0 -1 1)[ 1 -2 1 0]
!* 5- (-1 1 0 1)[ 1 1 -2 0]
!* 6- ( 1 -1 0 1)[ 1 1 -2 0]
!* Sort?
data lattice_sd(:, 7,3)/-1, 0, 0/ ; data lattice_sn(:, 7,3)/ 0,-1, 1/
data lattice_sd(:, 8,3)/ 0,-1, 0/ ; data lattice_sn(:, 8,3)/ 0, 1, 1/
data lattice_sd(:, 9,3)/ 1, 1, 0/ ; data lattice_sn(:, 9,3)/-1, 0, 1/
data lattice_sd(:,10,3)/-1, 0, 0/ ; data lattice_sn(:,10,3)/ 1, 0, 1/
data lattice_sd(:,11,3)/ 0,-1, 0/ ; data lattice_sn(:,11,3)/-1, 1, 1/
data lattice_sd(:,12,3)/ 1, 1, 0/ ; data lattice_sn(:,12,3)/ 1,-1, 1/
!*** Twin systems for HCP structures (3) ***
!* System {1012}<1011>
!* Sort?
!* MISSING: not implemented yet
!*** Slip-Slip interactions for HCP structures (3) ***
data lattice_SlipIntType( 1,1:lattice_MaxNslipOfStructure(3),3)/1,2,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType( 2,1:lattice_MaxNslipOfStructure(3),3)/2,1,2,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType( 3,1:lattice_MaxNslipOfStructure(3),3)/2,2,1,2,2,2,2,2,2,2,2,2/
data lattice_SlipIntType( 4,1:lattice_MaxNslipOfStructure(3),3)/2,2,2,1,2,2,2,2,2,2,2,2/
data lattice_SlipIntType( 5,1:lattice_MaxNslipOfStructure(3),3)/2,2,2,2,1,2,2,2,2,2,2,2/
data lattice_SlipIntType( 6,1:lattice_MaxNslipOfStructure(3),3)/2,2,2,2,2,1,2,2,2,2,2,2/
data lattice_SlipIntType( 7,1:lattice_MaxNslipOfStructure(3),3)/2,2,2,2,2,2,1,2,2,2,2,2/
data lattice_SlipIntType( 8,1:lattice_MaxNslipOfStructure(3),3)/2,2,2,2,2,2,2,1,2,2,2,2/
data lattice_SlipIntType( 9,1:lattice_MaxNslipOfStructure(3),3)/2,2,2,2,2,2,2,2,1,2,2,2/
data lattice_SlipIntType(10,1:lattice_MaxNslipOfStructure(3),3)/2,2,2,2,2,2,2,2,2,1,2,2/
data lattice_SlipIntType(11,1:lattice_MaxNslipOfStructure(3),3)/2,2,2,2,2,2,2,2,2,2,1,2/
data lattice_SlipIntType(12,1:lattice_MaxNslipOfStructure(3),3)/2,2,2,2,2,2,2,2,2,2,2,1/
!*** Twin-twin interactions for HCP structures (3) ***
! MISSING: not implemented yet
CONTAINS
!****************************************
!* - lattice_Init
!* - lattice_SchmidMatrices
!****************************************
subroutine lattice_init()
!**************************************
!* Module initialization *
!**************************************
call lattice_SchmidMatrices()
end subroutine
subroutine lattice_SchmidMatrices()
!**************************************
!* Calculation of Schmid matrices *
!**************************************
use prec, only: pReal,pInt
use math, only: math_I3,nrmMandel,mapMandel
implicit none
!* Definition of variables
integer(pInt) i,j,k,l
real(pReal) norm_d,norm_t,norm_n
!* Iteration over the lattice structures
do l=1,lattice_MaxLatticeStructure
!* Iteration over the slip systems
do k=1,lattice_MaxNslipOfStructure(l)
!* Definition of transverse direction st for the frame (sd,st,sn)
lattice_st(1,k,l)=lattice_sn(2,k,l)*lattice_sd(3,k,l)-lattice_sn(3,k,l)*lattice_sd(2,k,l)
lattice_st(2,k,l)=lattice_sn(3,k,l)*lattice_sd(1,k,l)-lattice_sn(1,k,l)*lattice_sd(3,k,l)
lattice_st(3,k,l)=lattice_sn(1,k,l)*lattice_sd(2,k,l)-lattice_sn(2,k,l)*lattice_sd(1,k,l)
norm_d=dsqrt(lattice_sd(1,k,l)**2+lattice_sd(2,k,l)**2+lattice_sd(3,k,l)**2)
norm_t=dsqrt(lattice_st(1,k,l)**2+lattice_st(2,k,l)**2+lattice_st(3,k,l)**2)
norm_n=dsqrt(lattice_sn(1,k,l)**2+lattice_sn(2,k,l)**2+lattice_sn(3,k,l)**2)
lattice_sd(:,k,l)=lattice_sd(:,k,l)/norm_d
lattice_st(:,k,l)=lattice_st(:,k,l)/norm_t
lattice_sn(:,k,l)=lattice_sn(:,k,l)/norm_n
!* Defintion of Schmid matrix
forall (i=1:3,j=1:3) lattice_Sslip(i,j,k,l)=lattice_sd(i,k,l)*lattice_sn(j,k,l)
!* Vectorization of normalized Schmid matrix
forall (i=1:6) lattice_Sslip_v(i,k,l) = nrmMandel(i)/2.0_pReal * &
(lattice_Sslip(mapMandel(1,i),mapMandel(2,i),k,l)+lattice_Sslip(mapMandel(2,i),mapMandel(1,i),k,l))
enddo
!* Iteration over the twin systems
do k=1,lattice_MaxNtwinOfStructure(l)
!* Definition of transverse direction tt for the frame (td,tt,tn)
lattice_tt(1,k,l)=lattice_tn(2,k,l)*lattice_td(3,k,l)-lattice_tn(3,k,l)*lattice_td(2,k,l)
lattice_tt(2,k,l)=lattice_tn(3,k,l)*lattice_td(1,k,l)-lattice_tn(1,k,l)*lattice_td(3,k,l)
lattice_tt(3,k,l)=lattice_tn(1,k,l)*lattice_td(2,k,l)-lattice_tn(2,k,l)*lattice_td(1,k,l)
norm_d=dsqrt(lattice_td(1,k,l)**2+lattice_td(2,k,l)**2+lattice_td(3,k,l)**2)
norm_t=dsqrt(lattice_tt(1,k,l)**2+lattice_tt(2,k,l)**2+lattice_tt(3,k,l)**2)
norm_n=dsqrt(lattice_tn(1,k,l)**2+lattice_tn(2,k,l)**2+lattice_tn(3,k,l)**2)
lattice_td(:,k,l)=lattice_td(:,k,l)/norm_d
lattice_tt(:,k,l)=lattice_tt(:,k,l)/norm_t
lattice_tn(:,k,l)=lattice_tn(:,k,l)/norm_n
!* Defintion of Schmid matrix and transformation matrices
lattice_Qtwin(:,:,k,l)=-math_I3
forall (i=1:3,j=1:3)
lattice_Stwin(i,j,k,l)=lattice_td(i,k,l)*lattice_tn(j,k,l)
lattice_Qtwin(i,j,k,l)=lattice_Qtwin(i,j,k,l)+2*lattice_tn(i,k,l)*lattice_tn(j,k,l)
endforall
!* Vectorization of normalized Schmid matrix
lattice_Stwin_v(1,k,l)=lattice_Stwin(1,1,k,l)
lattice_Stwin_v(2,k,l)=lattice_Stwin(2,2,k,l)
lattice_Stwin_v(3,k,l)=lattice_Stwin(3,3,k,l)
!* be compatible with Mandel notation of Tstar
lattice_Stwin_v(4,k,l)=(lattice_Stwin(1,2,k,l)+lattice_Stwin(2,1,k,l))/dsqrt(2.0_pReal)
lattice_Stwin_v(5,k,l)=(lattice_Stwin(2,3,k,l)+lattice_Stwin(3,2,k,l))/dsqrt(2.0_pReal)
lattice_Stwin_v(6,k,l)=(lattice_Stwin(1,3,k,l)+lattice_Stwin(3,1,k,l))/dsqrt(2.0_pReal)
enddo
enddo
end subroutine
END MODULE

3
trunk/math.f90
View File

@ -79,6 +79,7 @@
call random_seed() call random_seed()
call get_seed(seed) call get_seed(seed)
seed = 1
call halton_seed_set(seed) call halton_seed_set(seed)
call halton_ndim_set(3) call halton_ndim_set(3)
@ -306,7 +307,6 @@
InvA = math_identity2nd(dimen) InvA = math_identity2nd(dimen)
B = A B = A
CALL Gauss(dimen,B,InvA,LogAbsDetA,AnzNegEW,error) CALL Gauss(dimen,B,InvA,LogAbsDetA,AnzNegEW,error)
RETURN RETURN
END SUBROUTINE math_invert END SUBROUTINE math_invert
@ -1878,6 +1878,5 @@ endif
END FUNCTION END FUNCTION
END MODULE math END MODULE math

69
trunk/mpie_cpfem_marc2005r3.f90
View File

@ -33,8 +33,9 @@
include "math.f90" include "math.f90"
include "IO.f90" include "IO.f90"
include "mesh.f90" include "mesh.f90"
include "crystal.f90" include "lattice.f90"
include "constitutive.f90" include "constitutive.f90"
include "crystallite.f90"
include "CPFEM.f90" include "CPFEM.f90"
! !
@ -74,7 +75,37 @@
! s stress - should be updated by user ! s stress - should be updated by user
! t state variables (comes in at t=n, must be updated ! t state variables (comes in at t=n, must be updated
! to have state variables at t=n+1) ! to have state variables at t=n+1)
! dt increment of state variables ! dt ! Marc common blocks are in fixed format so they have to be pasted in here
!
! Marc common blocks are in fixed format so they have to be pasted in here
! Beware of changes in newer Marc versions -- these are from 2005r3
! concom is needed for inc, subinc, ncycle, lovl
! include 'concom'
common/concom/ &
iacous, iasmbl, iautth, ibear, icompl, iconj, icreep, ideva(50), idyn, idynt,&
ielas, ielcma, ielect, iform, ifour, iharm, ihcps, iheat, iheatt, ihresp,&
ijoule, ilem, ilnmom, iloren, inc, incext, incsub, ipass, iplres, ipois,&
ipoist, irpflo, ismall, ismalt, isoil, ispect, ispnow, istore, iswep, ithcrp,&
itherm, iupblg, iupdat, jacflg, jel, jparks, largst, lfond, loadup, loaduq,&
lodcor, lovl, lsub, magnet, ncycle, newtnt, newton, noshr, linear, ivscpl,&
icrpim, iradrt, ipshft, itshr, iangin, iupmdr, iconjf, jincfl, jpermg, jhour,&
isolvr, jritz, jtable, jshell, jdoubl, jform, jcentr, imini, kautth, iautof,&
ibukty, iassum, icnstd, icnstt, kmakmas, imethvp,iradrte,iradrtp, iupdate,iupdatp,&
ncycnt, marmen ,idynme, ihavca, ispf, kmini, imixed, largtt, kdoela, iautofg,&
ipshftp,idntrc, ipore, jtablm, jtablc, isnecma,itrnspo,imsdif, jtrnspo,mcnear,&
imech, imecht, ielcmat, ielectt,magnett, imsdift,noplas, jtabls, jactch, jtablth,&
kgmsto ,jpzo, ifricsh, iremkin,iremfor, ishearp,jspf, machining, jlshell,icompsol,&
iupblgfo,jcondir,nstcrp, nactive,ipassref, nstspnt,ibeart,icheckmpc, noline, icuring,&
ishrink,ioffsflg,isetoff, iharmt, inc_incdat, iautspc,ibrake
! creeps is needed for timinc (time increment)
! include 'creeps'
common/marc_creeps/ &
cptim,timinc,timinc_p,timinc_s,timincm,timinc_a,timinc_b,creept(33),icptim,icfte,icfst,&
icfeq,icftm,icetem,mcreep,jcreep,icpa,icftmp,icfstr,icfqcp,icfcpm,icrppr,icrcha,icpb,iicpmt,iicpa
increment of state variables
! ngens size of stress - strain law ! ngens size of stress - strain law
! n element number ! n element number
! nn integration point number ! nn integration point number
@ -132,15 +163,18 @@
integer(pInt) computationMode integer(pInt) computationMode
dimension e(*),de(*),t(*),dt(*),g(*),d(ngens,*),s(*), n(2),coord(ncrd,*),disp(ndeg,*),matus(2),dispt(ndeg,*),ffn(itel,*),& dimension e(*),de(*),t(*),dt(*),g(*),d(ngens,*),s(*), n(2),coord(ncrd,*),disp(ndeg,*),matus(2),dispt(ndeg,*),ffn(itel,*),&
frotn(itel,*),strechn(itel),eigvn(itel,*),ffn1(itel,*),frotn1(itel,*),strechn1(itel),eigvn1(itel,*),kcus(2) frotn(itel,*),strechn(itel),eigvn(itel,*),ffn1(itel,*),frotn1(itel,*),strechn1(itel),eigvn1(itel,*),kcus(2)
! Marc common blocks are in fixed format so they have to be pasted in here ! Marc common blocks are in fixed format so they have to be pasted in here
! Beware of changes in newer Marc versions -- these are from 2005r3 ! Beware of changes in newer Marc versions -- these are from 2005r3
! concom is needed for inc, subinc, ncycle, lovl ! concom is needed for inc, subinc, ncycle, lovl
! include 'concom' ! include 'concom'
common/concom/ & common/marc_concom/ &
iacous, iasmbl, iautth, ibear, icompl, iconj, icreep, ideva(50), idyn, idynt,& iacous, iasmbl, iautth, ibear, icompl, iconj, icreep, ideva(50), idyn, idynt,&
ielas, ielcma, ielect, iform, ifour, iharm, ihcps, iheat, iheatt, ihresp,& ielas, ielcma, ielect, iform, ifour, iharm, ihcps, iheat, iheatt, ihresp,&
ijoule, ilem, ilnmom, iloren, inc, incext, incsub, ipass, iplres, ipois,& ijoule, ilem, ilnmom, iloren, inc, incext, incsub, ipass, iplres, ipois,&
@ -155,7 +189,8 @@
imech, imecht, ielcmat, ielectt,magnett, imsdift,noplas, jtabls, jactch, jtablth,& imech, imecht, ielcmat, ielectt,magnett, imsdift,noplas, jtabls, jactch, jtablth,&
kgmsto ,jpzo, ifricsh, iremkin,iremfor, ishearp,jspf, machining, jlshell,icompsol,& kgmsto ,jpzo, ifricsh, iremkin,iremfor, ishearp,jspf, machining, jlshell,icompsol,&
iupblgfo,jcondir,nstcrp, nactive,ipassref, nstspnt,ibeart,icheckmpc, noline, icuring,& iupblgfo,jcondir,nstcrp, nactive,ipassref, nstspnt,ibeart,icheckmpc, noline, icuring,&
ishrink,ioffsflg,isetoff, iharmt, inc_incdat, iautspc,ibrake ishrink,ioffsflg,isetoff, ioffsetm,iharmt, inc_incdat,iautspc,ibrake, icbush ,istream_input,&
iprsinp,ivlsinp,ifirst_time,ipin_m,jgnstr_glb, imarc_return,iqvcinp,nqvceid,istpnx,imicro1
! creeps is needed for timinc (time increment) ! creeps is needed for timinc (time increment)
! include 'creeps' ! include 'creeps'
@ -167,22 +202,32 @@
if (inc == 0) then if (inc == 0) then
cycleCounter = 0 cycleCounter = 0
else else
if (theInc /= inc .or. theCycle /= ncycle .or. theLovl /= lovl) cycleCounter = cycleCounter+1 if (theCycle > ncycle) cycleCounter = 0 ! reset counter for each cutback
if (theCycle /= ncycle .or. theLovl /= lovl) cycleCounter = cycleCounter+1 ! ping pong
endif
if (cptim > theTime .or. theInc /= inc) then ! reached convergence
lastIncConverged = .true.
outdatedByNewInc = .true.
endif endif
if (theInc /= inc) outdatedByNewInc = .true.
if (mod(cycleCounter,2) /= 0) computationMode = 4 ! recycle if (mod(cycleCounter,2) /= 0) computationMode = 4 ! recycle
if (mod(cycleCounter,4) == 2) computationMode = 3 ! collect if (mod(cycleCounter,4) == 2) computationMode = 3 ! collect
if (mod(cycleCounter,4) == 0) computationMode = 2 ! compute if (mod(cycleCounter,4) == 0) computationMode = 2 ! compute
if (computationMode == 4 .and. ncycle == 0 .and. .not. lastIncConverged) &
computationMode = 6 ! recycle but restore known good consistent tangent
if (computationMode == 4 .and. lastIncConverged) then
computationMode = 5 ! recycle and record former consistent tangent
lastIncConverged = .false.
endif
if (computationMode == 2 .and. outdatedByNewInc) then if (computationMode == 2 .and. outdatedByNewInc) then
outdatedByNewInc = .false.
computationMode = 1 ! compute and age former results computationMode = 1 ! compute and age former results
outdatedByNewInc = .false.
endif endif
theInc = inc theTime = cptim ! record current starting time
theCycle = ncycle theInc = inc ! record current increment number
theLovl = lovl theCycle = ncycle ! record current cycle count
theLovl = lovl ! record current lovl
call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,n(1),nn,s,mod(theCycle,2_pInt*ijaco)==0,d,ngens) call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,n(1),nn,s,mod(theCycle,2_pInt*ijaco)==0,d,ngens)
@ -196,7 +241,7 @@
END SUBROUTINE END SUBROUTINE
! !
!
SUBROUTINE plotv(v,s,sp,etot,eplas,ecreep,t,m,nn,layer,ndi,nshear,jpltcd) SUBROUTINE plotv(v,s,sp,etot,eplas,ecreep,t,m,nn,layer,ndi,nshear,jpltcd)
!******************************************************************** !********************************************************************
! This routine sets user defined output variables for Marc ! This routine sets user defined output variables for Marc

30
trunk/mpie_cpfem_marc2007r1.f90
View File

@ -33,8 +33,9 @@
include "math.f90" include "math.f90"
include "IO.f90" include "IO.f90"
include "mesh.f90" include "mesh.f90"
include "crystal.f90" include "lattice.f90"
include "constitutive.f90" include "constitutive.f90"
include "crystallite.f90"
include "CPFEM.f90" include "CPFEM.f90"
! !
@ -135,7 +136,6 @@
dimension e(*),de(*),t(*),dt(*),g(*),d(ngens,*),s(*), n(2),coord(ncrd,*),disp(ndeg,*),matus(2),dispt(ndeg,*),ffn(itel,*),& dimension e(*),de(*),t(*),dt(*),g(*),d(ngens,*),s(*), n(2),coord(ncrd,*),disp(ndeg,*),matus(2),dispt(ndeg,*),ffn(itel,*),&
frotn(itel,*),strechn(itel),eigvn(itel,*),ffn1(itel,*),frotn1(itel,*),strechn1(itel),eigvn1(itel,*),kcus(2) frotn(itel,*),strechn(itel),eigvn(itel,*),ffn1(itel,*),frotn1(itel,*),strechn1(itel),eigvn1(itel,*),kcus(2)
@ -172,22 +172,32 @@
if (inc == 0) then if (inc == 0) then
cycleCounter = 0 cycleCounter = 0
else else
if (theInc /= inc .or. theCycle /= ncycle .or. theLovl /= lovl) cycleCounter = cycleCounter+1 if (theCycle > ncycle) cycleCounter = 0 ! reset counter for each cutback
if (theCycle /= ncycle .or. theLovl /= lovl) cycleCounter = cycleCounter+1 ! ping pong
endif
if (cptim > theTime .or. theInc /= inc) then ! reached convergence
lastIncConverged = .true.
outdatedByNewInc = .true.
endif endif
if (theInc /= inc) outdatedByNewInc = .true.
if (mod(cycleCounter,2) /= 0) computationMode = 4 ! recycle if (mod(cycleCounter,2) /= 0) computationMode = 4 ! recycle
if (mod(cycleCounter,4) == 2) computationMode = 3 ! collect if (mod(cycleCounter,4) == 2) computationMode = 3 ! collect
if (mod(cycleCounter,4) == 0) computationMode = 2 ! compute if (mod(cycleCounter,4) == 0) computationMode = 2 ! compute
if (computationMode == 4 .and. ncycle == 0 .and. .not. lastIncConverged) &
computationMode = 6 ! recycle but restore known good consistent tangent
if (computationMode == 4 .and. lastIncConverged) then
computationMode = 5 ! recycle and record former consistent tangent
lastIncConverged = .false.
endif
if (computationMode == 2 .and. outdatedByNewInc) then if (computationMode == 2 .and. outdatedByNewInc) then
outdatedByNewInc = .false.
computationMode = 1 ! compute and age former results computationMode = 1 ! compute and age former results
outdatedByNewInc = .false.
endif endif
theInc = inc theTime = cptim ! record current starting time
theCycle = ncycle theInc = inc ! record current increment number
theLovl = lovl theCycle = ncycle ! record current cycle count
theLovl = lovl ! record current lovl
call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,n(1),nn,s,mod(theCycle,2_pInt*ijaco)==0,d,ngens) call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,n(1),nn,s,mod(theCycle,2_pInt*ijaco)==0,d,ngens)
@ -201,7 +211,7 @@
END SUBROUTINE END SUBROUTINE
! !
!
SUBROUTINE plotv(v,s,sp,etot,eplas,ecreep,t,m,nn,layer,ndi,nshear,jpltcd) SUBROUTINE plotv(v,s,sp,etot,eplas,ecreep,t,m,nn,layer,ndi,nshear,jpltcd)
!******************************************************************** !********************************************************************
! This routine sets user defined output variables for Marc ! This routine sets user defined output variables for Marc

10
trunk/prec.f90
View File

@ -16,11 +16,17 @@
integer(pInt), parameter :: ijaco = 1_pInt ! frequency of FEM Jacobi update integer(pInt), parameter :: ijaco = 1_pInt ! frequency of FEM Jacobi update
integer(pInt), parameter :: nCutback = 10_pInt ! cutbacks in time-step integration integer(pInt), parameter :: nCutback = 10_pInt ! cutbacks in time-step integration
integer(pInt), parameter :: nReg = 1_pInt ! regularization attempts for Jacobi inversion integer(pInt), parameter :: nReg = 1_pInt ! regularization attempts for Jacobi inversion
real(pReal), parameter :: pert_Fg = 1.0e-5_pReal ! strain perturbation for FEM Jacobi real(pReal), parameter :: pert_Fg = 1.0e-6_pReal ! strain perturbation for FEM Jacobi
integer(pInt), parameter :: nOuter = 10_pInt ! outer loop limit integer(pInt), parameter :: nOuter = 10_pInt ! outer loop limit
integer(pInt), parameter :: nInner = 1000_pInt ! inner loop limit integer(pInt), parameter :: nInner = 1000_pInt ! inner loop limit
real(pReal), parameter :: reltol_Outer = 1.0e-4_pReal ! relative tolerance in outer loop (state) real(pReal), parameter :: reltol_Outer = 1.0e-6_pReal ! relative tolerance in outer loop (state)
real(pReal), parameter :: reltol_Inner = 1.0e-6_pReal ! relative tolerance in inner loop (Lp) real(pReal), parameter :: reltol_Inner = 1.0e-6_pReal ! relative tolerance in inner loop (Lp)
real(pReal), parameter :: abstol_Inner = 1.0e-8_pReal ! absolute tolerance in inner loop (Lp) real(pReal), parameter :: abstol_Inner = 1.0e-8_pReal ! absolute tolerance in inner loop (Lp)
!
real(pReal), parameter :: resToler = 1.0e-6_pReal ! relative tolerance of residual in GIA iteration
real(pReal), parameter :: resAbsol = 1.0e+0_pReal ! absolute tolerance of residual in GIA iteration (corresponds to 1 Pa)
real(pReal), parameter :: resBound = 1.0e+2_pReal ! relative maximum value (upper bound) for GIA residual
integer(pInt), parameter :: NRiterMax = 24_pInt ! maximum number of GIA iteration
real(pReal), parameter :: subStepMin = 1.0e-3_pReal ! minimum (relative) size of sub-step allowed during cutback
END MODULE prec END MODULE prec