c...Plane stress/strain advanced elasticity "main engine" c...for use with FELIPE version 3 - sample datafiles are: c... eladbm1.dat (3-noded beam in elastic block) c... eladft3.dat (footing on semi-infinite soil mass) c... eladcv2.dat (cavern excavated in infinite stressed rock mass) c...These illustrate use of different element and loading types c program eladv cb..."big" version of elast.for 2D elasticity program. cb...Includes coding for full range of element types in PREFEL, cb...and for loading by spec. displts, body forces and cb...excavation. Lines in CAPITALS are added in eladv, as are cb...comments starting cb... c parameter(melem=800,mpoin=2000,mnode=8, 1msets=4,mcomp=6,mtotv=4000,mfixv=150, 2mfron=150,mstif=11325,mbufa=100) 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 (e.g. Young's mod, Poisson ratio) c...mtotv = max no. of degrees of freedom (= 2*mpoin for elasticity) cb...mfixv = max no. of fixed degrees of freedom cb...mfron = max frontwidth in solver FRONT cb...mstif = (mfron+1)*mfron/2 cb...mbufa = size of buffer used in FRONT c 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),lmtyp(melem) 4 ,nvar(mpoin),nfvar(mfixv),spdis(mfixv),TLOAD(MTOTV) 5 ,strsg(4,4,melem) 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(l) = 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...ndof(np) = no. of dof.s at node np cp..nvar(np) is a pointer array for global dof of node np cb...nfvar(i) is the global dof of the i'th fixed deg of freedom cb...spdis(i) is the specified value of the i'th specified dof cb...tload(i) is total load so far (used with excavations) cb...strsg(i,ij,l) is stress component i at gausspt ij of element l 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... 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 c real*8 pival,eqrhs(mbufa) 1 ,vecrv(mfron),gstif(mstif),gload(mfron),equat(mfron,mbufa) DIMENSION NPIVO(MBUFA),NAMEV(MBUFA),NACVA(MFRON) cb...arrays above are used in the FRONT frontal solver c write(6,*)'Program ELADV: advanced elasticity.' write(6,*)'Handles all 2D finite and infinite elements,' write(6,*)'as well as 2-noded and 3-noded 1D beam elements.' write(6,*)'Dimensioned for',melem,' elements, and',mpoin,' nodes.' write(6,*) 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 write(6,*) 1'Is this a plane stress, plane strain or axisymmetric analysis?' 3 write(6,*) 1'Type 1 for plane stress, 2 for plane strain, or 3 for axisymm...' read(5,'(i1)',err=3) mstyp write(16,*) if (mstyp.eq.1) then write(6,*) 'Plane stress analysis' write(16,*)'Plane stress analysis' else if (mstyp.eq.2) then write(6,*) 'Plane strain analysis' write(16,*)'Plane strain analysis' else if (mstyp.eq.3) then write(6,*) 'Axisymmetric analysis' write(16,*)'Axisymmetric analysis' else write(6,*)'Invalid input!' goto 3 endif write(16,*) c cp...create the pointer array nvar. Needed because a node may cp...have 2 or 3 d.o.f.s; this info stored in ndof cp...also store info on fixed d.o.f.s in nfvar and spdis, cp...rather than in nkode ntotv = 0 nfixv = 0 do 5 np=1,npoin ndofnp = ndof(np) nvar(np) = ntotv kode = nkode(np) C...SORT OUT KODE AND STORE FIXED VARIABLE NO. IN NFVAR() KOD = KODE DO 4 M=1,NDOFNP KD = KOD/10 K = KOD - KD*10 KOD = KD IF(K.ne.0) then NFIXV=NFIXV+1 IF (NFIXV.GT.MFIXV) THEN WRITE(6,*) 1 'No. of fixed dof.s exceeds parameter MFIXV' write(6,*)'Increase MFIXV and recompile.' call prompt stop endif M1 = NDOFNP + 1 - M NFVAR(NFIXV) = ntotv + M1 spdis(nfixv) = 0.0 endif 4 continue ntotv = ntotv + ndofnp 5 continue nsdis = nfixv cp...s/r LOAd will add spec displts, total no will be held in nsdis cp...Now local dof i of node np is global dof nvar(np)+i cp...and nfvar(i) is global dof of i'th fixed dof cp...and ntotv is total no. of degrees of freedom c c c...open a scratch file to write the element stiffness matrices on open(unit=18,status='scratch',form='unformatted') write(6,*)'Scratch file opened ...' write(6,*) c cb...open further scratch files for use in FRONT OPEN(19,STATUS='SCRATCH',FORM='UNFORMATTED') OPEN(20,STATUS='SCRATCH',FORM='UNFORMATTED') OPEN(21,STATUS='SCRATCH',FORM='UNFORMATTED') PIVAL = 1.d-8 c call prompt write(6,*) write(6,*)'Starting subroutine STIFF...' call stiff(melem,mpoin,mnode,msets,mcomp,mstyp, 1 nelem,npoin,lnods,lnodn,lsidn,coord,lpros,props 2 ,lmtyp,ndof,nvar) c call prompt write(6,*) write(6,*)'Starting subroutine LOAD...' call load(melem,mpoin,mnode,msets,mcomp,mtotv,mstyp, 1 nelem,npoin,lnods,lnodn,lsidn,coord, 2 surft,aslod,lmtyp,ndof,nvar,tload,strsg, 3 lpros,props,mfixv,nfixv,nsdis,nfvar,spdis) c cb...the next subroutine replaces ASSMB and SOLVE in ELAST.FOR call prompt write(6,*) write(6,*)'Starting subroutine FRONT...' irsol = 0 cb...irsol=0 for first solution, irsol=1 for re-solution CALL FRONT(MTOTV,MPOIN,MELEM,MFIXV,MNODE,MFRON,MSTIF, 1NPOIN,NELEM,NSDIS,LNODS,LNODN,NDOF,NVAR,NFVAR,ASLOD,ASDIS,SPDIS, 2VECRV,GSTIF,GLOAD,EQUAT,NACVA,lsidn, 3PIVAL,IRSOL,EQRHS,NPIVO,NAMEV,MBUFA) c call prompt write(6,*) write(6,*)'Starting subroutine RESULT...' call result(melem,mpoin,mnode,msets,mcomp,mtotv,mstyp, 1 nelem,npoin,lnods,lnodn,lsidn,coord,nkode,lpros,props, 2 aslod,asdis,lmtyp,ndof,nvar,strsg,tload) c write(6,*) write(6,*)'ELADV 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