From 6d721dc16c337f678934734e080e8f1098f272bc Mon Sep 17 00:00:00 2001 From: Luc Hantcherli Date: Fri, 14 Mar 2008 21:56:46 +0000 Subject: [PATCH] Modifications of cpfem_marc2005r3 according the ones in cpfem_marc2007r1 --- trunk/mpie_cpfem_marc2005r3.f90 | 290 +++++++++++++++++--------------- 1 file changed, 154 insertions(+), 136 deletions(-) diff --git a/trunk/mpie_cpfem_marc2005r3.f90 b/trunk/mpie_cpfem_marc2005r3.f90 index 7364ade17..dfdf7f643 100644 --- a/trunk/mpie_cpfem_marc2005r3.f90 +++ b/trunk/mpie_cpfem_marc2005r3.f90 @@ -1,15 +1,15 @@ !******************************************************************** ! Material subroutine for MSC.Marc Version 0.1 ! -! written by F. Roters, P. Eisenlohr, L. Hantcherli, W.A. Counts +! written by F. Roters, P. Eisenlohr, L. Hantcherli, W.A. Counts ! MPI fuer Eisenforschung, Duesseldorf ! ! last modified: 08.11.2007 !******************************************************************** ! Usage: ! - choose material as hypela2 -! - set statevariable 2 to index of material -! - set statevariable 3 to index of texture +! - set statevariable 2 to index of material +! - set statevariable 3 to index of texture ! - choose output of user variables if desired ! - make sure the file "mattex.mpie" exists in the working ! directory @@ -27,108 +27,119 @@ ! - creeps: timinc !******************************************************************** ! - include "prec.f90" - include "debug.f90" + include "prec.f90" + include "FEsolving.f90" + include "debug.f90" include "math.f90" - include "IO.f90" - include "mesh.f90" - include "crystal.f90" - include "constitutive.f90" - include "CPFEM.f90" - - 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,& + include "IO.f90" + include "mesh.f90" + include "crystal.f90" + include "constitutive.f90" + include "CPFEM.f90" +! + + 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) 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 -! +! ************* user subroutine for defining material behavior ************** ! - use prec, only: pReal,pInt +! +! 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, ijaco + use FEsolving use CPFEM, only: CPFEM_general + use math, only: invnrmMandel + implicit real(pReal) (a-h,o-z) -! + + integer(pInt) computationMode + + 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) + + ! 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' +! 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,& @@ -144,37 +155,44 @@ 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/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 -! - 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) -! -! 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 -! 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 + 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 + + + if (inc == 0) then + cycleCounter = 0 + else + if (theInc /= inc .or. theCycle /= ncycle .or. theLovl /= lovl) cycleCounter = cycleCounter+1 + endif + + if (theInc /= inc) outdatedByNewInc = .true. + + if (mod(cycleCounter,2) /= 0) computationMode = 4 ! recycle + if (mod(cycleCounter,4) == 2) computationMode = 3 ! collect + if (mod(cycleCounter,4) == 0) computationMode = 2 ! compute + if (computationMode == 2 .and. outdatedByNewInc) then + outdatedByNewInc = .false. + computationMode = 1 ! compute and age former results + endif + + theInc = inc + theCycle = ncycle + theLovl = lovl + + call CPFEM_general(computationMode,ffn,ffn1,t(1),timinc,n(1),nn,s,mod(theCycle,2_pInt*ijaco)==0,d,ngens) + + +! Mandel: 11, 22, 33, SQRT(2)*12, SQRT(2)*23, SQRT(2)*13 +! Marc: 11, 22, 33, 12, 23, 13 + forall(i=1:ngens) d(1:ngens,i) = invnrmMandel(i)*d(1:ngens,i)*invnrmMandel(1:ngens) + s(1:ngens) = s(1:ngens)*invnrmMandel(1:ngens) + + return END SUBROUTINE ! @@ -199,21 +217,21 @@ ! ndi (3) number of direct stress components ! nshear (3) number of shear stress components ! -!******************************************************************** - use prec, only: pReal,pInt - use CPFEM, only: CPFEM_results, CPFEM_Nresults - use constitutive, only: constitutive_maxNresults - use mesh, only: mesh_FEasCP +!******************************************************************** + 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(*) real(pReal) v, t(*) - integer(pInt) m, nn, layer, ndi, nshear, jpltcd -! -! assign result variable - v=CPFEM_results(mod(jpltcd-1_pInt, CPFEM_Nresults+constitutive_maxNresults)+1_pInt,& - (jpltcd-1_pInt)/(CPFEM_Nresults+constitutive_maxNresults)+1_pInt,& - nn, mesh_FEasCP('elem', m)) + integer(pInt) m, nn, layer, ndi, nshear, jpltcd +! +! assign result variable + v=CPFEM_results(mod(jpltcd-1_pInt, CPFEM_Nresults+constitutive_maxNresults)+1_pInt,& + (jpltcd-1_pInt)/(CPFEM_Nresults+constitutive_maxNresults)+1_pInt,& + nn, mesh_FEasCP('elem', m)) return END SUBROUTINE !