changed structure of CPFEM_general to allow for non local algorithm and parrallelisation
This commit is contained in:
Franz Roters
parent
8619a92009
commit
896c37ede2
115
trunk/CPFEM.f90
115
trunk/CPFEM.f90
|
|
@ -22,11 +22,14 @@
|
|||
real(pReal), dimension (:,:,:,:,:), allocatable :: CPFEM_Fp_old
|
||||
real(pReal), dimension (:,:,:,:,:), allocatable :: CPFEM_Fp_new
|
||||
real(pReal), dimension (:,:,:,:), allocatable :: CPFEM_jaco_old
|
||||
real(pReal), dimension(6,6) :: CPFEM_dummy_jacobian
|
||||
integer(pInt) :: CPFEM_inc_old = 0_pInt
|
||||
integer(pInt) :: CPFEM_subinc_old = 1_pInt
|
||||
integer(pInt) :: CPFEM_cycle_old = -1_pInt
|
||||
integer(pInt) :: CPFEM_Nresults = 4_pInt ! three Euler angles plus volume fraction
|
||||
logical :: CPFEM_first_call = .true.
|
||||
|
||||
|
||||
CONTAINS
|
||||
|
||||
!*********************************************************
|
||||
|
|
@ -36,7 +39,7 @@
|
|||
SUBROUTINE CPFEM_init()
|
||||
!
|
||||
use prec, only: pReal,pInt
|
||||
use math, only: math_EulertoR
|
||||
use math, only: math_EulertoR, math_I3, math_identity2nd
|
||||
use mesh
|
||||
use constitutive
|
||||
!
|
||||
|
|
@ -46,8 +49,9 @@
|
|||
!
|
||||
! *** mpie.marc parameters ***
|
||||
allocate(CPFEM_Temperature (mesh_maxNips,mesh_NcpElems)) ; CPFEM_Temperature = 0.0_pReal
|
||||
allocate(CPFEM_ffn_all (3,3,mesh_maxNips,mesh_NcpElems)) ; CPFEM_ffn_all = 0.0_pReal
|
||||
allocate(CPFEM_ffn1_all (3,3,mesh_maxNips,mesh_NcpElems)) ; CPFEM_ffn1_all = 0.0_pReal
|
||||
allocate(CPFEM_ffn_all (3,3,mesh_maxNips,mesh_NcpElems))
|
||||
forall(e=1:mesh_NcpElems,i=1:mesh_maxNips) CPFEM_ffn_all(:,:,i,e) = math_I3
|
||||
allocate(CPFEM_ffn1_all (3,3,mesh_maxNips,mesh_NcpElems)) ; CPFEM_ffn1_all = CPFEM_ffn_all
|
||||
allocate(CPFEM_stress_all( 6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_stress_all = 0.0_pReal
|
||||
allocate(CPFEM_jacobi_all(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_jacobi_all = 0.0_pReal
|
||||
!
|
||||
|
|
@ -68,6 +72,9 @@
|
|||
! *** Old jacobian (consistent tangent) ***
|
||||
allocate(CPFEM_jaco_old(6,6,mesh_maxNips,mesh_NcpElems)) ; CPFEM_jaco_old = 0.0_pReal
|
||||
!
|
||||
! *** dummy Jacobian returned in odd cycles
|
||||
CPFEM_dummy_jacobian=1.0e50_pReal*math_identity2nd(6)
|
||||
|
||||
! *** Output to MARC output file ***
|
||||
write(6,*)
|
||||
write(6,*) 'Arrays allocated:'
|
||||
|
|
@ -93,49 +100,75 @@
|
|||
!*** perform initialization at first call, update variables and ***
|
||||
!*** call the actual material model ***
|
||||
!***********************************************************************
|
||||
SUBROUTINE CPFEM_general(ffn, ffn1, Temperature, CPFEM_inc, CPFEM_subinc, CPFEM_cn, CPFEM_dt, CPFEM_en, CPFEM_in)
|
||||
SUBROUTINE CPFEM_general(ffn, ffn1, Temperature, CPFEM_inc, CPFEM_subinc, CPFEM_cn, CPFEM_stress_recovery, CPFEM_dt,&
|
||||
CPFEM_en, CPFEM_in, CPFEM_stress, CPFEM_jaco, CPFEM_ngens)
|
||||
!
|
||||
use prec, only: pReal,pInt
|
||||
use math, only: math_init
|
||||
use mesh, only: mesh_init,mesh_FEasCP
|
||||
use math, only: math_init, invnrmMandel
|
||||
use mesh, only: mesh_init,mesh_FEasCP, mesh_NcpElems, FE_Nips, FE_mapElemtype, mesh_element
|
||||
use crystal, only: crystal_Init
|
||||
use constitutive, only: constitutive_init,constitutive_state_old,constitutive_state_new
|
||||
implicit none
|
||||
!
|
||||
real(pReal) ffn(3,3), ffn1(3,3), Temperature, CPFEM_dt
|
||||
integer(pInt) CPFEM_inc, CPFEM_subinc, CPFEM_cn, CPFEM_en, CPFEM_in, cp_en
|
||||
integer(pInt) CPFEM_inc, CPFEM_subinc, CPFEM_cn, CPFEM_en, CPFEM_in, cp_en, CPFEM_ngens, i, e
|
||||
real(pReal) ffn(3,3), ffn1(3,3), Temperature, CPFEM_dt, CPFEM_stress(CPFEM_ngens), CPFEM_jaco(CPFEM_ngens,CPFEM_ngens)
|
||||
logical CPFEM_stress_recovery
|
||||
!
|
||||
! initialization step
|
||||
if (CPFEM_first_call) then
|
||||
! three dimensional stress state ?
|
||||
call math_init()
|
||||
call mesh_init()
|
||||
call crystal_Init()
|
||||
call constitutive_init()
|
||||
call CPFEM_init()
|
||||
CPFEM_first_call = .false.
|
||||
endif
|
||||
if (CPFEM_inc==CPFEM_inc_old) then ! not a new increment
|
||||
! case of a new subincrement:update starting with subinc 2
|
||||
if (CPFEM_subinc > CPFEM_subinc_old) then
|
||||
CPFEM_sigma_old = CPFEM_sigma_new
|
||||
CPFEM_Fp_old = CPFEM_Fp_new
|
||||
constitutive_state_old = constitutive_state_new
|
||||
CPFEM_subinc_old = CPFEM_subinc
|
||||
endif
|
||||
else ! new increment
|
||||
CPFEM_sigma_old = CPFEM_sigma_new
|
||||
CPFEM_Fp_old = CPFEM_Fp_new
|
||||
constitutive_state_old = constitutive_state_new
|
||||
CPFEM_inc_old = CPFEM_inc
|
||||
CPFEM_subinc_old = 1_pInt
|
||||
endif
|
||||
|
||||
cp_en = mesh_FEasCP('elem',CPFEM_en)
|
||||
CPFEM_Temperature(CPFEM_in, cp_en) = Temperature
|
||||
CPFEM_ffn_all(:,:,CPFEM_in, cp_en) = ffn
|
||||
CPFEM_ffn1_all(:,:,CPFEM_in, cp_en) = ffn1
|
||||
call CPFEM_stressIP(CPFEM_cn, CPFEM_dt, cp_en, CPFEM_in)
|
||||
! calculate only every second cycle
|
||||
if(mod(CPFEM_cn,2)==0) then
|
||||
! really calculate only in first call of new cycle and when in stress recovery
|
||||
if(CPFEM_cn/=CPFEM_cycle_old .and. (CPFEM_stress_recovery .or. CPFEM_cn==0)) then
|
||||
! initialization step
|
||||
if (CPFEM_first_call) then
|
||||
! three dimensional stress state ?
|
||||
call math_init()
|
||||
call mesh_init()
|
||||
call crystal_Init()
|
||||
call constitutive_init()
|
||||
call CPFEM_init()
|
||||
CPFEM_Temperature = Temperature
|
||||
CPFEM_first_call = .false.
|
||||
endif
|
||||
if (CPFEM_inc==CPFEM_inc_old) then ! not a new increment
|
||||
! case of a new subincrement:update starting with subinc 2
|
||||
if (CPFEM_subinc > CPFEM_subinc_old) then
|
||||
CPFEM_sigma_old = CPFEM_sigma_new
|
||||
CPFEM_Fp_old = CPFEM_Fp_new
|
||||
constitutive_state_old = constitutive_state_new
|
||||
CPFEM_subinc_old = CPFEM_subinc
|
||||
endif
|
||||
else ! new increment
|
||||
CPFEM_sigma_old = CPFEM_sigma_new
|
||||
CPFEM_Fp_old = CPFEM_Fp_new
|
||||
constitutive_state_old = constitutive_state_new
|
||||
CPFEM_inc_old = CPFEM_inc
|
||||
CPFEM_subinc_old = 1_pInt
|
||||
endif
|
||||
CPFEM_cycle_old=CPFEM_cn
|
||||
! this shall be done in a parallel loop in the future
|
||||
do e=1,mesh_NcpElems
|
||||
do i=1,FE_Nips(FE_mapElemtype(mesh_element(2,e)))
|
||||
call CPFEM_stressIP(CPFEM_cn, CPFEM_dt, i, e)
|
||||
enddo
|
||||
enddo
|
||||
end if
|
||||
! return stress and jacobi
|
||||
! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13
|
||||
! Marc: 11, 22, 33, 12, 23, 13
|
||||
cp_en = mesh_FEasCP('elem', CPFEM_en)
|
||||
CPFEM_stress(1:CPFEM_ngens)=invnrmMandel(1:CPFEM_ngens)*CPFEM_stress_all(1:CPFEM_ngens, CPFEM_in, cp_en)
|
||||
CPFEM_jaco(1:CPFEM_ngens,1:CPFEM_ngens)=CPFEM_jaco_old(1:CPFEM_ngens,1:CPFEM_ngens, CPFEM_in, cp_en)
|
||||
forall(i=1:CPFEM_ngens) CPFEM_jaco(1:CPFEM_ngens,i)=CPFEM_jaco(1:CPFEM_ngens,i)*invnrmMandel(1:CPFEM_ngens)
|
||||
else
|
||||
! record data for use in second cycle and return fixed result
|
||||
cp_en = mesh_FEasCP('elem',CPFEM_en)
|
||||
CPFEM_Temperature(CPFEM_in, cp_en) = Temperature
|
||||
CPFEM_ffn_all(:,:,CPFEM_in, cp_en) = ffn
|
||||
CPFEM_ffn1_all(:,:,CPFEM_in, cp_en) = ffn1
|
||||
CPFEM_stress(1:CPFEM_ngens)=1.0e5_pReal
|
||||
CPFEM_jaco(1:CPFEM_ngens,1:CPFEM_ngens)=CPFEM_dummy_jacobian(1:CPFEM_ngens,1:CPFEM_ngens)
|
||||
end if
|
||||
return
|
||||
|
||||
END SUBROUTINE
|
||||
|
|
@ -147,8 +180,8 @@
|
|||
SUBROUTINE CPFEM_stressIP(&
|
||||
CPFEM_cn,& ! Cycle number
|
||||
CPFEM_dt,& ! Time increment (dt)
|
||||
cp_en,& ! Element number
|
||||
CPFEM_in) ! Integration point number
|
||||
CPFEM_in,& ! Integration point number
|
||||
cp_en) ! Element number
|
||||
|
||||
use prec, only: pReal,pInt,ijaco,nCutback
|
||||
use math, only: math_pDecomposition,math_RtoEuler, inDeg
|
||||
|
|
@ -169,7 +202,7 @@
|
|||
real(pReal), dimension(3,3,2) :: Fg,Fp
|
||||
real(pReal), dimension(constitutive_maxNstatevars,2) :: state
|
||||
|
||||
updateJaco = (mod(CPFEM_cn,ijaco)==0) ! update consistent tangent every ijaco'th iteration
|
||||
updateJaco = (mod(CPFEM_cn,2_pInt*ijaco)==0) ! update consistent tangent every ijaco'th iteration
|
||||
|
||||
CPFEM_stress_all(:,CPFEM_in,cp_en) = 0.0_pReal ! average Cauchy stress
|
||||
if (updateJaco) CPFEM_jaco_old(:,:,CPFEM_in,cp_en) = 0.0_pReal ! average consistent tangent
|
||||
|
|
|
|||
|
|
@ -4,7 +4,7 @@
|
|||
! written by F. Roters, P. Eisenlohr, L. Hantcherli, W.A. Counts
|
||||
! MPI fuer Eisenforschung, Duesseldorf
|
||||
!
|
||||
! last modified: 16.10.2007
|
||||
! last modified: 08.11.2007
|
||||
!********************************************************************
|
||||
! Usage:
|
||||
! - choose material as hypela2
|
||||
|
|
@ -121,9 +121,7 @@
|
|||
!
|
||||
!
|
||||
use prec, only: pReal,pInt
|
||||
use math, only: invnrmMandel, nrmMandel
|
||||
use mesh, only: mesh_FEasCP
|
||||
use CPFEM, only: CPFEM_general,CPFEM_stress_all, CPFEM_jaco_old
|
||||
use CPFEM, only: CPFEM_general
|
||||
implicit real(pReal) (a-h,o-z)
|
||||
!
|
||||
! Marc common blocks are in fixed format so they have to be pasted in here
|
||||
|
|
@ -152,35 +150,30 @@
|
|||
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
|
||||
!
|
||||
integer(pInt) cp_en, i
|
||||
!
|
||||
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)
|
||||
!
|
||||
! call general material routine only in cycle 0 and for lovl==6 (stress recovery)
|
||||
|
||||
! subroutine cpfem_general(mpie_ffn, mpie_ffn1, mpie_cn, mpie_tinc, mpie_enp, mpie_in)
|
||||
!********************************************************************
|
||||
! This routine calculates the material behaviour
|
||||
!********************************************************************
|
||||
! mpie_ffn deformation gradient for t=t0
|
||||
! mpie_ffn1 deformation gradient for t=t1
|
||||
! mpie_cn number of cycle
|
||||
! mpie_tinc time increment
|
||||
! mpie_en element number
|
||||
! mpie_in intergration point number
|
||||
! subroutine cpfem_general(mpie_ffn, mpie_ffn1, temperature, mpie_inc, mpie_subinc, mpie_cn,
|
||||
! mpie_stress_recovery, mpie_tinc, mpie_en, mpie_in, mpie_s, mpie_d, mpie_ngens)
|
||||
!********************************************************************
|
||||
if ((lovl==6).or.(ncycle==0)) then
|
||||
call CPFEM_general(ffn, ffn1, t(1), inc, incsub, ncycle, timinc, n(1), nn)
|
||||
endif
|
||||
! return stress and jacobi
|
||||
! Mandel: 11, 22, 33, 12, 23, 13
|
||||
! Marc: 11, 22, 33, 12, 23, 13
|
||||
cp_en = mesh_FEasCP('elem', n(1))
|
||||
s(1:ngens)=invnrmMandel(1:ngens)*CPFEM_stress_all(1:ngens, nn, cp_en)
|
||||
d(1:ngens,1:ngens)=CPFEM_jaco_old(1:ngens,1:ngens, nn, cp_en)
|
||||
forall(i=1:ngens) d(1:ngens,i)=d(1:ngens,i)*invnrmMandel(1:ngens)
|
||||
return
|
||||
! This routine calculates the material behaviour
|
||||
!********************************************************************
|
||||
! mpie_ffn deformation gradient for t=t0
|
||||
! mpie_ffn1 deformation gradient for t=t1
|
||||
! temperature temperature
|
||||
! mpie_inc increment number
|
||||
! mpie_subinc subincrement number
|
||||
! mpie_cn number of cycle
|
||||
! mpie_stress_recovery indicates wether we are in stiffness assemly(lovl==4) or stress recovery(lovl==6)
|
||||
! mpie_tinc time increment
|
||||
! mpie_en element number
|
||||
! mpie_in intergration point number
|
||||
! mpie_s stress vector in Marc notation, i.e. 11 22 33 12, 23, 13
|
||||
! mpie_d jacoby in Marc notation
|
||||
! mpie_ngens size of stress strain law
|
||||
!********************************************************************
|
||||
call CPFEM_general(ffn, ffn1, t(1), inc, incsub, ncycle, stress_recovery, timinc, n(1), nn, s, d, ngens)
|
||||
return
|
||||
|
||||
END SUBROUTINE
|
||||
!
|
||||
|
|
|
|||
|
|
@ -4,7 +4,7 @@
|
|||
! written by F. Roters, P. Eisenlohr, L. Hantcherli, W.A. Counts
|
||||
! MPI fuer Eisenforschung, Duesseldorf
|
||||
!
|
||||
! last modified: 28.03.2007
|
||||
! last modified: 08.11.2007
|
||||
!********************************************************************
|
||||
! Usage:
|
||||
! - choose material as hypela2
|
||||
|
|
@ -120,10 +120,8 @@
|
|||
!3 continue
|
||||
!
|
||||
!
|
||||
use prec, only: pReal,pInt
|
||||
use math, only: invnrmMandel, nrmMandel
|
||||
use mesh, only: mesh_FEasCP
|
||||
use CPFEM, only: CPFEM_general,CPFEM_stress_all, CPFEM_jaco_old
|
||||
use prec, only: pReal,pInt
|
||||
use CPFEM, only: CPFEM_general
|
||||
implicit real(pReal) (a-h,o-z)
|
||||
!
|
||||
! Marc common blocks are in fixed format so they have to be pasted in here
|
||||
|
|
@ -153,35 +151,38 @@
|
|||
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
|
||||
!
|
||||
integer(pInt) cp_en, i
|
||||
!
|
||||
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)
|
||||
|
||||
logical stress_recovery
|
||||
|
||||
if(lovl==6) then
|
||||
stress_recovery = .true.
|
||||
else
|
||||
stress_recovery = .false.
|
||||
endif
|
||||
!
|
||||
! call general material routine only in cycle 0 and for lovl==6 (stress recovery)
|
||||
|
||||
! subroutine cpfem_general(mpie_ffn, mpie_ffn1, mpie_cn, mpie_tinc, mpie_enp, mpie_in)
|
||||
! subroutine cpfem_general(mpie_ffn, mpie_ffn1, temperature, mpie_inc, mpie_subinc, mpie_cn,
|
||||
! mpie_stress_recovery, mpie_tinc, mpie_en, mpie_in, mpie_s, mpie_d, mpie_ngens)
|
||||
!********************************************************************
|
||||
! This routine calculates the material behaviour
|
||||
!********************************************************************
|
||||
! mpie_ffn deformation gradient for t=t0
|
||||
! mpie_ffn1 deformation gradient for t=t1
|
||||
! mpie_cn number of cycle
|
||||
! mpie_tinc time increment
|
||||
! mpie_en element number
|
||||
! mpie_in intergration point number
|
||||
! mpie_ffn deformation gradient for t=t0
|
||||
! mpie_ffn1 deformation gradient for t=t1
|
||||
! temperature temperature
|
||||
! mpie_inc increment number
|
||||
! mpie_subinc subincrement number
|
||||
! mpie_cn number of cycle
|
||||
! mpie_stress_recovery indicates wether we are in stiffness assemly(lovl==4) or stress recovery(lovl==6)
|
||||
! mpie_tinc time increment
|
||||
! mpie_en element number
|
||||
! mpie_in intergration point number
|
||||
! mpie_s stress vector in Marc notation, i.e. 11 22 33 12, 23, 13
|
||||
! mpie_d jacoby in Marc notation
|
||||
! mpie_ngens size of stress strain law
|
||||
!********************************************************************
|
||||
if ((lovl==6).or.(ncycle==0)) then
|
||||
call CPFEM_general(ffn, ffn1, t(1), inc, incsub, ncycle, timinc, n(1), nn)
|
||||
endif
|
||||
! return stress and jacobi
|
||||
! Mandel: 11, 22, 33, 12, 23, 13
|
||||
! Marc: 11, 22, 33, 12, 23, 13
|
||||
cp_en = mesh_FEasCP('elem', n(1))
|
||||
s(1:ngens)=invnrmMandel(1:ngens)*CPFEM_stress_all(1:ngens, nn, cp_en)
|
||||
d(1:ngens,1:ngens)=CPFEM_jaco_old(1:ngens,1:ngens, nn, cp_en)
|
||||
forall(i=1:ngens) d(1:ngens,i)=d(1:ngens,i)*invnrmMandel(1:ngens)
|
||||
call CPFEM_general(ffn, ffn1, t(1), inc, incsub, ncycle, stress_recovery, timinc, n(1), nn, s, d, ngens)
|
||||
return
|
||||
|
||||
END SUBROUTINE
|
||||
|
|
|
|||
Loading…
Reference in New Issue