c...Plane beam analysis "main engine" c...for use with FELIPE version 3 c... - sample datafile is framexa1.dat or framexb1.dat c program frame c parameter(melem=100,mpoin=200,mnode=2, 1msets=4,mcomp=3,mtotv=600) c 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) c...mtotv = max no. of degrees of freedom (= 3*mpoin for beams) 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), 4 surft(msets,2),lmtyp(melem),diagk(mtotv), 6 nvar(mpoin),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 ( = 2 for beam element) c...lsidn = no. of sides in element l ( = 1 for beam element) 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 cb...nvar(np) indicates the global d.o.f. at which node np starts 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...astmp is a scratch vector used in the equation solution substitutions c...surft(m,1),surft(m,2) = normal,tangential surface tractions in set m c...gstif holds assembled global stiffness matrix as a vector in skyline storage c... or as a matrix using band storage 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...lmtyp(l) = material type of elt l (not used in this program) c...diagk holds the diagonal of the global stiffness matrix, to use as c... preconditioner in the conjugate gradient solution c...z,r,p,y are scratch vectors used during the c.g. solution c write(6,*)'Program FRAME: solves beam and frame problems.' write(6,*)'Handles 2-noded cubic 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 iflag = 0 do 2 l=1,nelem nnode = lnodn(l) nside = lsidn(l) nodsid = 10*nnode + nside if (nodsid.ne.21) then iflag = 1 write(6,*)'ERROR: Cannot handle element ',l write(6,*)'which has ',nside,' sides and ',nnode,' nodes' endif 2 continue if (iflag.eq.1) then call prompt stop endif c cgtol = 1.e-5 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 nvar(1) = 0 ntotv = 0 do 5 np=1,npoin ndofnp = ndof(np) nvar(np) = ntotv ntotv = ntotv + ndofnp 5 continue cp...Now local dof i of node np is global dof nvar(np)+i 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 call prompt write(6,*) write(6,*)'Starting subroutine STIFF...' call stiff(melem,mpoin,mnode,msets,mcomp,mstyp, 1 nelem,npoin,lnods,lnodn,coord,lpros,props) c call prompt write(6,*) write(6,*)'Starting subroutine LOAD...' call load(melem,mpoin,mnode,msets,mcomp,mtotv,mstyp, 1 nelem,npoin,lnods,lnodn,coord, 2 surft,aslod,lpros,props) c call prompt c...store the diagonal of K, to use as the preconditioner write(6,*) write(6,*)'Starting subroutine ASSDK...' call assdk(melem,mpoin,mnode,mtotv, 1 nelem,npoin,lnods,lnodn,ndof,nvar,nkode,diagk) c c...solve using conjugate gradients! call prompt write(6,*) write(6,*)'Starting subroutine PCGSOL...' call pcgsol(mtotv,melem,mpoin,mnode,ntotv,nelem, 1 npoin,lnods,lnodn,diagk,nkode,ndof,nvar, 2 aslod,asdis, z , r , p , y ,cgtol,itscg,rnorm0) c c...calculate nodal reactions by forming gstif*asdis in diagk, c...then add this to negative of equivalent nodal loads in aslod call prompt write(6,*) write(6,*)'Calculating reactions at nodes...' call ematgv(18,mtotv,melem,mpoin,mnode,ntotv,nelem, 1 npoin,lnods,lnodn,ndof,nvar,asdis,diagk) do 10 i=1,ntotv aslod(i) = -aslod(i) + diagk(i) 10 continue c call prompt write(6,*) write(6,*)'Starting subroutine RESULT...' call result(melem,mpoin,mnode,msets,mcomp,mtotv,mstyp, 1 nelem,npoin,lnods,lnodn,coord,nkode,lpros,props, 2 aslod,asdis) c write(6,*) write(6,*)'FRAME completed. Results are in file ',prtfil(1:lenin) write(6,*) write(6,*) close(15) close(16) close(17) call prompt stop 999 write(6,*) 1'Cannot open file with this name. Terminating program.' call prompt stop end c c c subroutine stiff(melem,mpoin,mnode,msets,mcomp,mstyp, 1 nelem,npoin,lnods,lnodn,coord,lpros,props) c c...calculates the element stiffness matrices, and writes them c...to the scratch file. c dimension props(msets,mcomp),lnods(melem,mnode),encod(2,8), 1 estif(6,6), 2 lpros(melem),lnodn(melem),coord(mpoin,2) c...encod holds the nodal coords for the nodes around the current element c...estif will be the stiffness matrix for the current element c rewind 18 c c...start loop over all elements do 80 l=1,nelem c c...extract Young's modulus be and moment of inertia bi from props lset = lpros(l) be = props(lset,1) bi = props(lset,2) c...also extract cross-sectional area ba of beam ba = props(lset,3) if (ba.eq.0.0) ba = 1.0 c c...put nodal coords in encod nnode = lnodn(l) do 15 n=1,nnode nn = lnods(l,n) encod(1,n) = coord(nn,1) encod(2,n) = coord(nn,2) 15 continue c...zero stiffness matrix do 20 j=1,6 do 20 i=1,6 estif(i,j) = 0.0 20 continue c c...find length, and angle of inclination of beam dx = encod(1,2) - encod(1,1) dy = encod(2,2) - encod(2,1) bl = sqrt(dx*dx + dy*dy) alpha = atan2(dy,dx) c = cos(alpha) s = sin(alpha) c2 = c*c s2 = s*s cs = c*s bl2 = bl*bl bil = bi/bl bil2 = bi/bl2 c c...now assemble stiffness matrix estif(1,1) = ba*c2 + 12.*bil2*s2 estif(1,2) = (ba - 12.*bil2)*cs estif(1,3) = - 6.*bil*s estif(1,4) = -(ba*c2 + 12.*bil2*s2) estif(1,5) = -(ba - 12.*bil2)*cs estif(1,6) = - 6.*bil*s estif(2,2) = ba*s2 + 12.*bil2*c2 estif(2,3) = 6.*bil*c estif(2,4) = -(ba - 12.*bil2)*cs estif(2,5) = -(ba*s2 + 12.*bil2*c2) estif(2,6) = 6.*bil*c estif(3,3) = 4.*bi estif(3,4) = 6.*bil*s estif(3,5) = - 6.*bil*c estif(3,6) = 2.*bi estif(4,4) = ba*c2 + 12.*bil2*s2 estif(4,5) = (ba - 12.*bil2)*cs estif(4,6) = 6.*bil*s estif(5,5) = ba*s2 + 12.*bil2*c2 estif(5,6) = - 6.*bil*c estif(6,6) = 4.*bi c c...all terms are multiplied by E/l bel = be/bl do 25 i=1,6 do 25 j=i,6 estif(i,j) = bel*estif(i,j) 25 continue c c...now fill in lower triangle of estif do 30 i=2,6 do 30 j=1,i-1 estif(i,j) = estif(j,i) 30 continue c c...write estif to file write (18) estif write(6,*)'Stiffness matrix for elt ',l,' written to file.' 80 continue c...end of element loop c write(6,*)'Elt stiffness matrices written to file' endfile 18 c return end