diff --git a/trunk/CPFEM.f90 b/trunk/CPFEM.f90 index 06b945e1e..b9107d266 100644 --- a/trunk/CPFEM.f90 +++ b/trunk/CPFEM.f90 @@ -105,18 +105,18 @@ CPFEM_jacobi_all = 0.0_pReal ! ! *** User defined results !!! MISSING incorporate consti_Nresults *** - allocate(CPFEM_results(CPFEM_Nresults+constitutive_maxNresults,texture_maxNgrains,mesh_maxNips,mesh_NcpElems)) + allocate(CPFEM_results(CPFEM_Nresults+constitutive_maxNresults,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) CPFEM_results = 0.0_pReal ! ! *** Second Piola-Kirchoff stress tensor at (t=t0) and (t=t1) *** - allocate(CPFEM_sigma_old(6,texture_maxNgrains,mesh_maxNips,mesh_NcpElems)) - allocate(CPFEM_sigma_new(6,texture_maxNgrains,mesh_maxNips,mesh_NcpElems)) + allocate(CPFEM_sigma_old(6,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) + allocate(CPFEM_sigma_new(6,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) CPFEM_sigma_old = 0.0_pReal CPFEM_sigma_new = 0.0_pReal ! ! *** Plastic deformation gradient at (t=t0) and (t=t1) *** - allocate(CPFEM_Fp_old(3,3,texture_maxNgrains,mesh_maxNips,mesh_NcpElems)) - allocate(CPFEM_Fp_new(3,3,texture_maxNgrains,mesh_maxNips,mesh_NcpElems)) + allocate(CPFEM_Fp_old(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) + allocate(CPFEM_Fp_new(3,3,constitutive_maxNgrains,mesh_maxNips,mesh_NcpElems)) CPFEM_Fp_old = 0.0_pReal CPFEM_Fp_new = 0.0_pReal ! @@ -187,7 +187,7 @@ ! ! *** Initialization of the matrices for t=t0 *** ! data from constitutive? - vf = constitutive_volfrac(iori,CPFEM_in,cp_en) !ÄÄÄ + vf = constitutive_matVolFrac(iori,CPFEM_in,cp_en)*constitutive_texVolFrac(iori,CPFEM_in,cp_en) !ÄÄÄ ! *** Calculation of the solution at t=t1 *** ! QUESTION use the mod() as flag parameter in the call ?? @@ -282,10 +282,10 @@ real(pReal) cs(6), cd(6,6), CPFEM_dt integer(pInt) cp_en ,CPFEM_in, iori, ising, icut, iconv, isjaco ! *** Local variables *** - real(pReal) Fp_old(3,3), Fp_new(3,3), state_old(constitutive_Nstatevars) - real(pReal) state_new(constitutive_Nstatevars), Tstar_v(6), CPFEM_ffn(3,3), CPFEM_ffn1(3,3) - real(pReal) Tstar_v_h(6), state_new_h(constitutive_Nstatevars), phi1, PHI, phi2, dt_i, delta_Fg(3,3), Fg_i(3,3) - real(pReal) state_new_i(constitutive_Nstatevars), time + real(pReal) Fp_old(3,3), Fp_new(3,3), state_old(constitutive_Nstatevars(iori, CPFEM_in, cp_en)) + real(pReal) state_new(constitutive_Nstatevars(iori, CPFEM_in, cp_en)), Tstar_v(6), CPFEM_ffn(3,3), CPFEM_ffn1(3,3) + real(pReal) Tstar_v_h(6), state_new_h(constitutive_Nstatevars(iori, CPFEM_in, cp_en)), phi1, PHI, phi2, dt_i + real(pReal) delta_Fg(3,3), Fg_i(3,3), state_new_i(constitutive_Nstatevars(iori, CPFEM_in, cp_en)), time integer(pInt) jcut ! *** Numerical parameters *** integer(pInt), parameter :: ncut=7_pInt @@ -404,14 +404,14 @@ ! ! *** Definition of variables *** ! *** Subroutine parameters *** + integer(pInt) cp_en, CPFEM_in, iori, ising, icut, iconv, isjaco real(pReal) cs(6), dcs_de(6,6), dt, phi1, PHI, phi2, Fg_old(3,3), Fg_new(3,3) - real(pReal) Fp_old(3,3), Fp_new(3,3), state_old(constitutive_Nstatevars) - real(pReal) state_new(constitutive_Nstatevars), Tstar_v(6) - integer(pInt) cp_en, CPFEM_in, iori, ising, icut, iconv, isjaco + real(pReal) Fp_old(3,3), Fp_new(3,3), state_old(constitutive_Nstatevars(iori, CPFEM_in, cp_en)) + real(pReal) state_new(constitutive_Nstatevars(iori, CPFEM_in, cp_en)), Tstar_v(6) ! *** Local variables *** integer(pInt) ic real(pReal) Fe(3,3), R(3,3), U(3,3), dev(6), dF(3,3), Fg2(3,3), sgm2(6) - real(pReal) state2(constitutive_Nstatevars), Fp2(3,3), cs1(6) + real(pReal) state2(constitutive_Nstatevars(iori, CPFEM_in, cp_en)), Fp2(3,3), cs1(6) ! *** Numerical parameters *** real(pReal), parameter :: pert_ct=1.0e-5_pReal ! *** Error treatment *** @@ -483,19 +483,19 @@ !*** NEWTON-RAPHSON Calculation *** !*********************************************************************** use prec, only: pReal,pInt - use constitutive, only: constitutive_Nstatevars + use constitutive, only: constitutive_Nstatevars, constitutive_HomogenizedC, constitutive_dotState use math implicit none ! *** Definition of variables *** ! *** Subroutine parameters *** - real(pReal) dt,Fg_old(3,3),Fg_new(3,3),Fp_old(3,3),Fp_new(3,3), Fe(3,3) - real(pReal) state_old(constitutive_Nstatevars), state_new(constitutive_Nstatevars) - real(pReal) Tstar_v(6), cs(6) integer(pInt) cp_en, CPFEM_in, iori, iconv, ising + real(pReal) dt,Fg_old(3,3),Fg_new(3,3),Fp_old(3,3),Fp_new(3,3), Fe(3,3) + real(pReal) state_old(constitutive_Nstatevars(iori, CPFEM_in, cp_en)), state_new(constitutive_Nstatevars(iori, CPFEM_in, cp_en)) + real(pReal) Tstar_v(6), cs(6) ! *** Local variables *** real(pReal) crite, tol_in, tol_out, invFp_old(3,3), det, A(3,3), C_66(6,6), Lp(3,3), dLp(3,3,6) real(pReal) I3tLp(3,3), help(3,3), help1(6), Tstar0_v(6), R1(6), norm1, tdLp(3,3) - real(pReal) dstate(constitutive_Nstatevars), R2(6), norm2, invFp_new(3,3), Estar(3,3) + real(pReal) dstate(constitutive_Nstatevars(iori, CPFEM_in, cp_en)), R2(6), norm2, invFp_new(3,3), Estar(3,3) real(pReal) Estar_v(6), Jacobi(6,6), invJacobi(6,6), dTstar_v(6) integer(pInt) iouter, iinner , dummy, err, i, j, k ! *** Numerical parameters *** @@ -505,11 +505,11 @@ real(pReal), parameter :: tol_inner = 1.0e-3_pReal real(pReal), parameter :: eta = 13.7_pReal - crite=eta*constitutive_s0_slip/constitutive_n_slip !ÄÄÄ +! crite=eta*constitutive_s0_slip/constitutive_n_slip !ÄÄÄ ! ! *** Tolerances *** - tol_in = tol_inner*s0_slip !ÄÄÄ - tol_out = tol_outer*s0_slip !ÄÄÄ +! tol_in = tol_inner*s0_slip !ÄÄÄ +! tol_out = tol_outer*s0_slip !ÄÄÄ ! ! *** Error treatment *** iconv = 0 @@ -527,7 +527,7 @@ ! *** Calculation of A and T*0 (see Kalidindi) *** A = matmul(Fg_new,invFp_old) A = matmul(transpose(A), A) - C_66=constitutive_homogenizedC(iori, CPFEM_in, cp_en) !ÄÄÄ + C_66=constitutive_HomogenizedC(iori, CPFEM_in, cp_en) !ÄÄÄ ! ! *** Second level of iterative procedure: Resistences *** do iouter=1,nouter @@ -554,7 +554,7 @@ enddo enddo enddo - Jacobi= matmul(C_66, help1) + math_identity(6) +! Jacobi= matmul(C_66, help1) + math_identity(6) call math_invert6x6(Jacobi, invJacobi, dummy, err) !ÄÄÄ if (err==1_pInt) then forall (i=1:6) Jacobi(i,i)=1.05d0*maxval(Jacobi(i,:)) ! regularization diff --git a/trunk/mpie_cpfem_marc.f90 b/trunk/mpie_cpfem_marc.f90 index 60ff60080..b3b827ba0 100644 --- a/trunk/mpie_cpfem_marc.f90 +++ b/trunk/mpie_cpfem_marc.f90 @@ -4,7 +4,7 @@ ! written by F. Roters, P. Eisenlohr, L. Hantcherli, W.A. Counts ! MPI fuer Eisenforschung, Duesseldorf ! -! last modified: 26.03.2007 +! last modified: 28.03.2007 !******************************************************************** ! Usage: ! - choose material as hypela2 @@ -32,102 +32,91 @@ include "mesh.f90" include "constitutive.f90" include "CPFEM.f90" -! - subroutine hypela2(d,g,e,de,s,t,dt,ngens,n,nn,kc,mats,ndi,nshear,& - disp,dispt,coord,ffn,frotn,strechn,eigvn,ffn1,& - frotn1,strechn1,eigvn1,ncrd1,itel,ndeg1,ndm,& - nnode,jtype,lclass,ifr,ifu) +! + subroutine hypela2(d,g,e,de,s,t,dt,ngens,n,nn,kcus,matus,ndi,& + nshear,disp,dispt,coord,ffn,frotn,strechn,eigvn,ffn1,& + frotn1,strechn1,eigvn1,ncrd,itel,ndeg,ndm,& + nnode,jtype,lclass,ifr,ifu) !******************************************************************** ! This is the Marc material routine !******************************************************************** ! -! ************* user subroutine for defining material behavior ************** -! -! -!CAUTION: Due to calculation of the Deformation gradients, Stretch Tensors and -! Rotation tensors at previous and current states, the analysis can be -! computationally expensive. Please use the user subroutine -> hypela -! if these kinematic quantities are not needed in the constitutive model -! -! -! IMPORTANT NOTES : -! -! -! (1) You must include the -> process,1,1,1, card in the parameter section -! of MARC input deck. -! -! (2) For total Lagrangian formulation use the -> 'large disp' card in the -! parameter section of MARC input deck. -! For updated Lagrangian formulation use the -> 'large disp' and 'update' -! cards in the parameter section of MARC input deck. However for any -! large strain calculation (whether elasticity or inelasticity) must entail -! the use of 'finite' parameter card also. -! -! (3) For Plasticity, the 2nd or 3rd cards in 'geometry' option in the model -! definition sections must be flagged for correct behavior in incompressible -! deformation. -! -! (4) The kinematic quantities are calculated for the following continuum -! elements (both lower and higher order) : -! plane stress, plane strain, generalized plane strain, axisymmetric, -! axisymmetric with twist and brick elements. -! -! -! -! -! d stress strain law to be formed -! g change in stress due to temperature effects -! e total strain -! de increment of strain -! s stress - should be updated by user -! t state variables (comes in at t=n, must be updated -! to have state variables at t=n+1) -! dt increment of state variables -! ngens size of stress - strain law -! n element number -! nn integration point number -! kc layer number -! mats material identification number -! ndi number of direct components -! nshear number of shear components -! disp incremental displacements -! dispt displacements at t=n (at assembly, lovl=4) and -! displacements at t=n+1 (at stress recovery, lovl=6) -! coord coordinates -! ncrd number of coordinates -! ndeg number of degrees of freedom -! itel dimension of F and R, either 2 or 3 -! nnode number of nodes per element -! jtype element type -! lclass element class -! ifr set to 1 if R has been calculated -! ifu set to 1 if strech has been calculated -! -! at t=n : -! -! ffn deformation gradient -! frotn rotation tensor -! strechn square of principal stretch ratios, lambda(i) -! eigvn(i,j) i principal direction components for j eigenvalues -! -! at t=n+1 : -! -! ffn1 deformation gradient -! frotn1 rotation tensor -! strechn1 square of principal stretch ratios, lambda(i) -! eigvn1(i,j) i principal direction components for j eigenvalues -! -! The following operation obtains U (stretch tensor) at t=n+1 : -! -! call scla(un1,0.d0,itel,itel,1) -! do 3 k=1,3 -! do 2 i=1,3 -! do 1 j=1,3 -! un1(i,j)=un1(i,j)+dsqrt(strechn1(k))*eigvn1(i,k)*eigvn1(j,k) -!1 continue -!2 continue -!3 continue -! +! ************* user subroutine for defining material behavior ************** +! +! +! CAUTION : Due to calculation of the Deformation gradients, Stretch Tensors and +! Rotation tensors at previous and current states, the analysis can be +! computationally expensive. Please use the user subroutine -> hypela +! if these kinematic quantities are not needed in the constitutive model +! +! +! IMPORTANT NOTES : +! +! (1) F,R,U are only available for continuum and membrane elements (not for +! shells and beams). +! +! (2) For total Lagrangian formulation use the -> 'Elasticity,1' card(= +! total Lagrange with large disp) in the parameter section of input deck. +! For updated Lagrangian formulation use the -> 'Plasticity,3' card(= +! update+finite+large disp+constant d) in the parameter section of +! input deck. +! +! +! d stress strain law to be formed +! g change in stress due to temperature effects +! e total elastic strain +! de increment of strain +! s stress - should be updated by user +! t state variables (comes in at t=n, must be updated +! to have state variables at t=n+1) +! dt increment of state variables +! ngens size of stress - strain law +! n element number +! nn integration point number +! kcus(1) layer number +! kcus(2) internal layer number +! matus(1) user material identification number +! matus(2) internal material identification number +! ndi number of direct components +! nshear number of shear components +! disp incremental displacements +! dispt displacements at t=n (at assembly, lovl=4) and +! displacements at t=n+1 (at stress recovery, lovl=6) +! coord coordinates +! ncrd number of coordinates +! ndeg number of degrees of freedom +! itel dimension of F and R, either 2 or 3 +! nnode number of nodes per element +! jtype element type +! lclass element class +! ifr set to 1 if R has been calculated +! ifu set to 1 if strech has been calculated +! +! at t=n : +! +! ffn deformation gradient +! frotn rotation tensor +! strechn square of principal stretch ratios, lambda(i) +! eigvn(i,j) i principal direction components for j eigenvalues +! +! at t=n+1 : +! +! ffn1 deformation gradient +! frotn1 rotation tensor +! strechn1 square of principal stretch ratios, lambda(i) +! eigvn1(i,j) i principal direction components for j eigenvalues +! +! The following operation obtains U (stretch tensor) at t=n+1 : +! +! call scla(un1,0.d0,itel,itel,1) +! do 3 k=1,3 +! do 2 i=1,3 +! do 1 j=1,3 +! un1(i,j)=un1(i,j)+dsqrt(strechn1(k))*eigvn1(i,k)*eigvn1(j,k) +!1 continue +!2 continue +!3 continue +! ! use prec, only: pReal,pInt use CPFEM, only : CPFEM_stress_all, CPFEM_jaco_old @@ -162,8 +151,8 @@ real(pReal) mpie_timefactor, mpie_stress(ngens) real(pReal) mpie_jacobi(ngens,ngens) ! - dimension e(*),de(*),t(*),dt(*),g(*),d(ngens,ngens),s(ngens),n(2),coord(ncrd,*),disp(ndeg,*),dispt(ndeg,*),ffn(itel,*),& - frotn(itel,*),strechn(itel),eigvn(itel,*),ffn1(itel,*),frotn1(itel,*),strechn1(itel),eigvn1(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) ! ! call general material routine only in increment 0 and for lovl==6 (stress recovery) @@ -182,13 +171,13 @@ ! mpie_in intergration point number ! mpie_dimension dimension of stress/strain vector !******************************************************************** - cp_en=mesh_mapFEtoCPelement(n(1)) + cp_en=mesh_FEasCP('elem', n(1)) if ((lovl==6).or.(inc==0)) then call cpfem_general(ffn, ffn1, inc, incsub, ncycle, timinc, cp_en, nn) endif ! return stress and jacobi - s=CPFEM_stress_all(1:ngens, nn, cp_en) - d=CPFEM_jaco_old(1:ngens,1:ngens, nn, cp_en) + s(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) return end @@ -217,6 +206,8 @@ !******************************************************************** use prec, only: pReal,pInt use CPFEM, only: CPFEM_results, CPFEM_Nresults + use constitutive, only: constitutive_maxNresults + use mesh, only: mesh_FEasCP implicit none ! real(pReal) s(*),etot(*),eplas(*),ecreep(*),sp(*) @@ -224,23 +215,23 @@ integer(pInt) m, nn, layer, ndi, nshear, jpltcd ! ! assign result variable - v=CPFEM_result(mod(jpltcd, CPFEM_Nresults+constitutive_Nresults),& - int(jpltcd/(CPFEM_Nresults+constitutive_Nresults)),& - nn, mesh_mapFEtoCPelement(m)) + v=CPFEM_results(mod(jpltcd, CPFEM_Nresults+constitutive_maxNresults),& + int(jpltcd/(CPFEM_Nresults+constitutive_maxNresults)),& + nn, mesh_FEasCP('elem', m)) return end ! ! - subroutine utimestep(timestep,timestepold,icall,time,timeloadcase) +! subroutine utimestep(timestep,timestepold,icall,time,timeloadcase) !******************************************************************** ! This routine modifies the addaptive time step of Marc !******************************************************************** - use prec, only: pReal,pInt - use CPFEM, only : CPFEM_timefactor_max - implicit none +! use prec, only: pReal,pInt +! use CPFEM, only : CPFEM_timefactor_max +! implicit none ! - real(pReal) timestep, timestepold, time,timeloadcase - integer(pInt) icall +! real(pReal) timestep, timestepold, time,timeloadcase +! integer(pInt) icall ! ! user subroutine for modifying the time step in auto step ! @@ -261,22 +252,22 @@ ! ! user coding ! reduce timestep during increment in case mpie_timefactor is too large - if(icall==2_pInt) then - if(mpie_timefactor_max>1.25_pReal) then - timestep=min(timestep,timestepold*0.8_pReal) - end if - return +! if(icall==2_pInt) then +! if(mpie_timefactor_max>1.25_pReal) then +! timestep=min(timestep,timestepold*0.8_pReal) +! end if +! return ! modify timestep at beginning of new increment - else if(icall==3_pInt) then - if(mpie_timefactor_max<=0.8_pReal) then - timestep=min(timestep,timestepold*1.25_pReal) - else if (mpie_timefactor_max<=1.0_pReal) then - timestep=min(timestep,timestepold/mpie_timefactor_max) - else if (mpie_timefactor_max<=1.25_pReal) then - timestep=min(timestep,timestepold*1.01_pReal) - else - timestep=min(timestep,timestepold*0.8_pReal) - end if - end if - return - end +! else if(icall==3_pInt) then +! if(mpie_timefactor_max<=0.8_pReal) then +! timestep=min(timestep,timestepold*1.25_pReal) +! else if (mpie_timefactor_max<=1.0_pReal) then +! timestep=min(timestep,timestepold/mpie_timefactor_max) +! else if (mpie_timefactor_max<=1.25_pReal) then +! timestep=min(timestep,timestepold*1.01_pReal) +! else +! timestep=min(timestep,timestepold*0.8_pReal) +! end if +! end if +! return +! end