c...Plane strain elasto-plasticity "main engine" c...for use with FELIPE version 3 - sample datafile is plcyl4.dat c...pladv is same as plast, but includes a range of solution methods c program pladv c...pladv uses row-based skyline storage of the global stiffness matrix c...and allows solution by: c... 0) T Choleski decomposition c... 1) L.D.L decomposition c... 2) conjugate gradients with diagonal preconditioning c... 3) conjugate gradients with Incomplete Choleski preconditioning c cc...This program is for advanced elasticity analyses (ndofn=-2) using cc...the theory of plasticity. It uses Mohr-Coulomb yield function, and cc...is written for associated flow rule. cc...Plane strain analysis, with material thickness = 1.0 c CC...THE PROGRAM (AND FELSUB.FOR) SHOULD BE COMPILED IN DOUBLE PRECISION: CC... ftn77 /dreal plast.for in Salford Fortran77 c parameter(melem=100,mpoin=400,mnode=8, 1msets=4,mcomp=6,mtotv=800,mstif=80000) c...these parameters fix the dimensions of the arrays: c...melem = max. no. of elements c...mpoin = max. no. of nodes c...mnode = max. no. of nodes in a single element c...msets = max. no. of mat. prop. sets (also used for loading sets) c...mcomp = max. no. of mat. props = 6 for (visco)plasticity c...mtotv = max no. of degrees of freedom (= 2*mpoin for elasticity) c...mstif = max. length of vector gstif storing global matrix using c skyline storage c character char*1 character*80 title,infile,outfil,prtfil c dimension lnods(melem,mnode),lnodn(melem),lsidn(melem), 1 coord(mpoin,2),nkode(mpoin),aslod(mtotv),ndof(mpoin), 2 lpros(melem),props(msets,mcomp),asdis(mtotv), 3 surft(msets,2),gstif(mstif),lmtyp(melem), 4 nminvc(mtotv),lpoint(mtotv),astmp(mtotv),diagk(mtotv) 5 ,z(mtotv),r(mtotv),p(mtotv),y(mtotv) c...lnods(l,n) = global no. of n'th node of element l c...lnodn(l) = no. of nodes in element l ( = 8 for 8-noded quad.) c...lsidn = no. of sides in element l ( = 4 for 8-noded quad.) c...coord(np,1), coord(np,2) = x,y coords of node np c...nkode(np) = fixity code of node np c...aslod is global nodal load vector c...astmp is a scratch vector used in the equation solution substitutions c...ndof(np) = no. of dof.s at node np (not used in this program) c...lpros(l) = mat. prop. set no. of element l c...props(m,i) = value of i'th mat. property in mat. prop. set m c...asdis is global nodal displacement vector c...surft(m,1),surft(m,2) = normal,tangential surface tractions in set m c...gstif holds assembled global stiffness matrix in compact form c...lmtyp(l) = material type of elt l (not used in this program) c...nminvc(i)=j means that G(i,j) is the first nonzero entry in row i of gstif c...lpoint(i)=ii means that G(i,i) is stored at gstif(ii) in skyline storage c dimension strsg(3,4,melem),strsg1(3,4,melem),totld(mtotv), 1 accld(mtotv),accds(mtotv),aslod1(mtotv),resld(mtotv), 2 fincs(9),mohrg(4,melem) cc..strsg holds stresses at Gauss points (Cartesian) cc..totld is the total load to be applied cc..accld is the accumulated load applied so far cc..aslod is the load being applied in the current increment cc..accds is the accumulated displacements cc..resld is the residual force vector cc..fincs is the set of load increments cc..mohrg is the yield code at Gauss pts: 0=elastic, 1=plastic c cc..tolerances for convergence tests: gtoler = 1.e-4 ftoler = 1.e-6 stoler = 1.e-8 cc..gtoler is tolerance for the global (displacement) iterations cc..ftoler is tolerance for the yield function F on yield surface cc..stoler is tolerance for stresses in local (stress) iterations cc..ftoler and stoler are used in subroutine stress c c write(6,*)'Program PLADV: advanced plasticity (Mohr-Coulomb).' write(6,*)'Handles 8-noded quadrilateral elements.' write(6,*)'Solution by Choleski, LDLt or conjugate gradients' write(6,*)'with diagonal or Incomplete Cholesk preconditioning.' write(6,*)'Dimensioned for',melem,' elements, and',mpoin,' nodes.' write(6,*) c c c...read an existing .dat file c 1 continue c c...get the input file name in infile call askfil(infile,iflag) if (iflag.eq.1) then call yesno(char,'Do you want to exit?') if (char.eq.'y') stop goto 1 endif c write(6,*)'Opening file ',infile,'...' lenin = lenstr(infile) open(15,file=infile(1:lenin),status='old',err=999) write(6,*) lm4 = lenin - 4 prtfil = infile(1:lm4)//'.prt' write(6,*)'Creating file ',prtfil(1:lenin),' for results...' write(6,*) open(16,file=prtfil(1:lenin),status='unknown',err=999) outfil = infile(1:lm4)//'.out' write(6,*)'Creating file ',outfil(1:lenin),' for output...' write(6,*) open(17,file=outfil(1:lenin),status='unknown',err=999) c c...read the input mesh data, and echo it to the output files call input(15,melem,mpoin,mnode,msets,mcomp,nelem, 1 npoin,lnods,lnodn,lsidn,lmtyp,coord,nkode,lpros,props 2 ,ndofn,npros,title,ndof) call print(16,melem,mpoin,mnode,msets,mcomp,nelem,npoin, 1 lnods,lnodn,lsidn,lmtyp,coord,nkode,lpros,props,ndofn, 2 npros,title,ndof) call output(17,melem,mpoin,mnode,msets,mcomp,nelem,npoin, 1 lnods,lnodn,lsidn,lmtyp,coord,nkode,lpros,props,ndofn, 2 npros,title,ndof) call yesno(char,'Do you want the mesh details listed?') if (char.eq.'y') then call print(6,melem,mpoin,mnode,msets,mcomp,nelem,npoin, 1 lnods,lnodn,lsidn,lmtyp,coord,nkode,lpros,props,ndofn, 2 npros,title,ndof) endif c iflag = 0 do 3 l=1,nelem nnode = lnodn(l) nside = lsidn(l) nodsid = 10*nnode + nside if (nodsid.ne.84) then iflag = 1 write(6,*)'ERROR: Cannot handle element ',l write(6,*)'which has ',nside,' sides and ',nnode,' nodes' endif 3 continue if (iflag.eq.1) then call prompt stop endif c write(6,*)'You can use a direct or an iterative solution method.' write(6,*)'Choose the parameter ISOLA from the following:' write(6,*) write(6,*)'ISOLA = 0 for direct solution using Choleski (LLt)' write(6,*)' decomposition and skyline storage;' write(6,*)'ISOLA = 1 for direct solution using L.D.Lt' write(6,*)' decomposition and skyline storage;' write(6,*)'ISOLA = 2 for iterative solution using conjugate' write(6,*)' gradients with diagonal preconditioning;' write(6,*)'ISOLA = 3 for Incomplete Choleski preconditioning' write(6,*)' with conjugate gradient solution.' write(6,*) write(6,*)'Type choice of ISOLA now:' read(5,*) isola if (isola.gt.1) then write(6,*)'Enter tolerance for convergence of c.g. iterations:' read(5,*) cgtol if (isola.eq.3) then write(6,*)'Enter factor for determining if fill-in will occur:' read(5,*) factor endif endif c cc..zero the accumulating variables, stresses, etc. ntotv = 2*npoin c do 10 i=1,mtotv accds(i) = 0.0 accld(i) = 0.0 resld(i) = 0.0 10 continue do 11 i=1,9 fincs(i) = 0.0 11 continue do 12 l=1,melem do 12 ij=1,4 do 12 i=1,3 strsg(i,ij,l) = 0.0 strsg1(i,ij,l) = 0.0 12 continue do 13 l=1,melem do 13 ij=1,4 mohrg(ij,l) = 0 13 continue cc..read the equivalent load vector, into totld, and then cc..read the load increment sizes call prompt write(6,*) write(6,*)'Starting subroutine LOAD...' call load(melem,mpoin,mnode,msets,mcomp,mtotv, 1 nelem,npoin,lnods,lnodn,coord, 2 surft,totld) c c...read boundary data header (no boundary set data expected) read(15,*) read(15,*)nbrad,nbdir if (nbrad.gt.0.or.nbdir.gt.0) then write(6,*) 1 'ERROR: boundary set data present - cannot cope!' call prompt stop endif c...read data defining load increments read(15,*) read(15,*) nincs,(fincs(i),i=1,nincs) c c c...open a scratch file to write the element stiffness matrices on open(unit=18,status='scratch',form='unformatted') c c c...analyse the "mask" of the global stiffness matrix, i.e. the range c...of non-zero values in each row. This information is stored in c...nminvc: row i starts in column nminvc(i). This is used to construct c...lpoint: g(i,i) will be located at gstif(lpoint(i)) call mask(melem,mpoin,mnode,mtotv,mstif,ntotv, 1 nelem,npoin,lnods,lnodn,nminvc,lpoint) c...doing this now, avoids repeating it within assmb at each iteration! c c cc..start of load increment loop cc..============================ c do 200 inc=1,nincs finc = fincs(inc) write(6,*)' Load increment ',inc write(6,*) write(6,*)' Increment size: ',finc write(6,*) cc..increment load, adding residual from previous increment do 20 i=1,ntotv asdis(i) = 0.0 aslod(i) = finc*totld(i) accld(i) = accld(i) + aslod(i) aslod(i) = aslod(i) + resld(i) 20 continue do 21 l=1,nelem do 21 ij=1,4 do 21 i=1,3 strsg1(i,ij,l) = strsg(i,ij,l) 21 continue cc..strsg will stay fixed as the start-of-increment stress cc..strsg1 will hold the iterates for the end-of-increment stress c kount = 0 cc..start of global iteration for displacements cc..=========================================== 100 continue kount = kount + 1 cc..note: stiffness depends on current stress state call stiff(melem,mpoin,mnode,msets,mcomp, 1 nelem,npoin,lnods,lnodn,coord,lpros,props 2 ,strsg1,mohrg) c c...don't need to assemble gstif if using PCG with diagonal precond. if (isola.ne.2) then call assmb(melem,mpoin,mnode,mtotv,mstif,ntotv, 1 nelem,npoin,lnods,lnodn, 2 gstif,nminvc,lpoint) endif c c...store the diagonal of K: used as preconditioner if isola=2, c...and used to calculate reactions call assdk(melem,mpoin,mnode,mtotv, 1 nelem,npoin,lnods,lnodn,nkode,diagk) c cc..test residual, to see if iteration has converged call check(mtotv,melem,mpoin,mnode, 1 ntotv,nelem,npoin,nkode,lnods,lnodn, 1 asdis,aslod,resld,rnorm,nminvc,lpoint) write(6,*)' Residual norm = ',rnorm write(6,*) if (kount.ge.2.and.rnorm.lt.gtoler) goto 150 c c do 110 i=1,ntotv aslod1(i) = aslod(i) 110 continue c if (isola.le.1) then c...direct solution using Choleski decomposition call solve(mpoin,mtotv,mstif,npoin,nkode, 1 aslod1,asdis,astmp,gstif,nminvc,lpoint,isola) else c...iterative solution using preconditioned conjugate gradients call pcgsol(mtotv,melem,mpoin,mnode,ntotv,nelem, 1 npoin,lnods,lnodn,diagk,nkode,mstif,gstif,nminvc,lpoint, 2 aslod1,asdis, z , r , p , y , astmp, 3 cgtol,itscg,rnorm0,isola,factor) endif c cc..now evaluate new stresses, from this displt increment, in strsg1 cc..and update mohrg if new points have yielded call stress(melem,mpoin,mnode,msets,mcomp,mtotv, 1 nelem,npoin,lnods,lnodn,coord,lpros,props, 2 asdis,strsg,strsg1,mohrg,ftoler,stoler) c goto 100 cc..end of global iteration cc..======================= 150 continue write(6,*) kount,' global iterations' write(6,*) c c...update total displacements and stresses do 160 i=1,ntotv accds(i) = accds(i) + asdis(i) 160 continue do 161 l=1,nelem do 161 ij=1,4 do 161 i=1,3 strsg(i,ij,l) = strsg1(i,ij,l) 161 continue c cc..compute residual force vector call residu(melem,mpoin,mnode,mtotv,nelem,npoin,ntotv, 1 lnods,lnodn,coord,nkode,accld,strsg,resld) c 200 continue cc..end of load increment loop cc..========================== c call result(melem,mpoin,mnode,msets,mcomp,mtotv, 1 nelem,npoin,lnods,lnodn,coord,nkode,lpros,props, 2 accld,accds,strsg,mohrg) c write(6,*) write(6,*)'PLADV completed. Results are in file ',prtfil(1:lenin) write(6,*) write(6,*) call prompt stop 999 write(6,*) 1'Cannot open file with this name. Terminating program.' call prompt stop end