diff --git a/trunk/CPFEM.f90 b/trunk/CPFEM.f90 index 8f144fae3..cbbf4e02e 100644 --- a/trunk/CPFEM.f90 +++ b/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 diff --git a/trunk/mpie_cpfem_marc2005r3.f90 b/trunk/mpie_cpfem_marc2005r3.f90 index 86720a73d..c712ba784 100644 --- a/trunk/mpie_cpfem_marc2005r3.f90 +++ b/trunk/mpie_cpfem_marc2005r3.f90 @@ -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 ! diff --git a/trunk/mpie_cpfem_marc2007r1.f90 b/trunk/mpie_cpfem_marc2007r1.f90 index 6a32e250c..489693f94 100644 --- a/trunk/mpie_cpfem_marc2007r1.f90 +++ b/trunk/mpie_cpfem_marc2007r1.f90 @@ -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