program beam2d!character*64 fname1,fname2dimension mea(100,4),jz(50,2),aai(10,2),aeu(10,2),qq(100) dimension r(6,6),rt(6,6),ra(6,6),pop(6),po(6)dimension akep(6,6),ake(6,6),p(300),x(100),y(100) dimension ld(300),is(6),a(50000),upe(6),ppe(6),ue(6)open(1,file='fname1_b.txt')open(2,file='fname2_b.txt')read(1,*) nn,ne,nm,na,ncwrite(*,*) NN,NE,NM,NA,NCnf=3nd=2nfd=nf*ndn=nn*nfdo 50 i=1,nn50 read(1,*) k,x(i),y(i),(p(nf*(i-1)+j),j=1,nf)read(1,*) ((jz(i,j),j=1,2),i=1,nc)do 100 i=1,ne100 read(1,*) ie,(mea(i,j),j=1,4),qq(i)read(1,*) ((aai(i,j),j=1,2),i=1,na)read(1,*) ((aeu(i,j),j=1,2),i=1,nm)close(1)write(2,460) nn,ne,nm,na,ncwrite(2,465) (i,x(i),y(i),p(3*i-2),p(3*i-1),p(3*i),i=1,nn)write(2,470) ((jz(i,j),j=1,2),i=1,nc)write(2,475) (i,(mea(i,j),j=1,4),i=1,ne)write(2,480) ((aai(i,j),j=1,2),i=1,na)write(2,485) ((aeu(i,j),j=1,2),i=1,nm)460 format(/5x,'the input of nn,ne,nm,na,nc'//(5x,5i5)) 465 format(/5x,'the input of x,y,p'//(5x,i5,5f10.2))470 format(/5x,'the input of jz'//(5x,i5,5x,i5))475 format(/5x,'the input of mea,qq'//(5x,i5,5x,4i5,f10.2)) 480 format(/5x,'the input of aai'//(5x,i5,5x,2e15.6))485 format(/5x,'the input of aeu'//(5x,i5,5x,2e15.6))!c structure stiffness statementcall fld(nn,ne,mea,nf,nd,n,nt,ld)do 500 i=1,nt500 a(i)=0do 600 ie=1,necall kep(nn,ie,mea,aeu,aai,x,y,akep)call cr(nn,ie,mea,x,y,r)call tran(6,6,r,rt)call dot(6,6,6,rt,akep,ra)call dot(6,6,6,ra,r,ake)call fis(ie,mea,nf,nd,nfd,is)do 560 i=1,nfddo 560 j=1,nfdif(is(i)-is(j)) 560,520,520520 ni=is(i)ij=ld(ni)-ni+is(j)a(ij)=a(ij)+ake(i,j)560 continue600 continue! c equivalent node forcedo 700 ie=1,necall fixf(nn,ne,ie,mea,x,y,qq,pop)call cr(nn,ie,mea,x,y,r)call tran(6,6,r,rt)call dot(6,6,1,rt,pop,po)call fis(ie,mea,nf,nd,nfd,is)do 650 i=1,nfdni=is(i)p(ni)=p(ni)-po(i)650 continue700 continuecall fcc(nc,n,nt,nf,jz,ld,a)call decom(n,nt,a,p,ld)write(2,850) (i,p(3*i-2),p(3*i-1),p(3*i),i=1,nn) write(2,860)do 800 ie=1,necall kep(nn,ie,mea,aeu,aai,x,y,akep)call cr(nn,ie,mea,x,y,r)call fis(ie,mea,nf,nd,nfd,is)do 750 i=1,nfdni=is(i)ue(i)=p(ni)750 continuecall dot(6,nfd,1,r,ue,upe)call dot(6,6,1,akep,upe,ppe)call fixf(nn,ne,ie,mea,x,y,qq,pop)do 760 i=1,nfdppe(i)=ppe(i)+pop(i)760 continuewrite(2,880) ie,ppe(1),ppe(2),ppe(3)write(2,900) ppe(4),ppe(5),ppe(6)800 continue850 format(//25x,'nodal displacement'//5x,'node no.',12x,'u',12x,'v',12x,'theta'/(7x,i5,3e15.6)) 860 format(//15x,'axial forces of element',/11x,'element no.','i-j',10x,'nx',12x,'ny',12x,'mz') 880 format(4x,i6,2x,'i',3e15.6)900 format(12x,'j',3e15.6)close(2)end!c is dimensionSUBROUTINE FIS(IE,MEA,NF,ND,NFD,IS)DIMENSION MEA(100,4),IS(NFD)DO 150 ID=1,NDDO 100 JF=1,NFI1=(ID-1)*NF+JFIS(I1)=(MEA(IE,ID)-1)*NF+JF100 CONTINUE150 CONTINUERETURNEND!c ld dimensionSUBROUTINE FLD(NN,NE,MEA,NF,ND,N,NT,LD)DIMENSION MEA(100,4),LD(N)LD(1)=1DO 300 K=1,NNNMIN=1000DO 200 I=1,NEDO 150 J=1,NDIF(MEA(I,J).NE.K) GOTO 150DO 100 L=1,NDIF (MEA(I,L).LT.NMIN) NMIN=MEA(I,L)100 CONTINUE150 CONTINUE200 CONTINUEDO 250 I=1,NFJ=(K-1)*NF+IIF(J.NE.1) LD(J)=LD(J-1)+(K-NMIN)*NF+I 250 CONTINUE300 CONTINUENT=LD(N)RETURNEND! c constrained conditionSUBROUTINE FCC(NC,N,NT,NF,JZ,LD,A)DIMENSION JZ(50,2),LD(N),A(NT)DO 350 K=1,NCI=JZ(K,1)J=JZ(K,2)NI=(I-1)*NF+JNJ=LD(NI)A(NJ)=1.0E25350 CONTINUERETURNEND! c linear equation by ldlt method SUBROUTINE DECOM(N,NT,A,B,LD)DIMENSION A(NT),B(N),LD(N)DO 20 I=1,NLDI=LD(I)IF (I.EQ.1) GOTO 10IO=LDI-IIM1=I-1MI=1-IO+LD(IM1)IF(MI.GT.IM1) GOTO 10DO 8 J=MI,IM1JO=LD(J)-JIJ=IO+JIF(J.EQ.1) GOTO 6JM1=J-1MJ=1-JO+LD(JM1)MIJ=MAX0(MI,MJ)IF(MIJ.GT.JM1) GOTO 6S=0.0DO 2 K=MIJ,JM1IK=IO+KJK=JO+KS=S+A(IK)*A(JK)2 CONTINUEA(IJ)=A(IJ)-S6 B(I)=B(I)-A(IJ)*B(J)8 CONTINUE10 ALDI=A(LDI)IF(I.EQ.1.OR.MI.GT.IM1) GOTO 13DO 12 J=MI,IM1IJ=IO+JLDJ=LD(J)S=A(IJ)A(IJ)=S/A(LDJ)ALDI=ALDI-S*A(IJ)12 CONTINUE13 A(LDI)=ALDIB(I)=B(I)/ALDI20 CONTINUEDO 40 LDI=2,NI=N-LDI+2IO=LD(I)-IMI=1-IO+LD(I-1)IM1=I-1IF(MI.GT.IM1) GOTO 40DO 30 J=MI,IM1IJ=IO+JB(J)=B(J)-A(IJ)*B(I)30 CONTINUE40 CONTINUERETURNENDSUBROUTINE zero(M,N,A)DIMENSION A(M,N)DO 10 I=1,MDO 10 J=1,NA(I,J)=0.010 CONTINUERETURNENDsubroutine tran(m,n,a,b)dimension a(m,n),b(n,m)do 20 i=1,mdo 20 j=1,nb(j,i)=a(i,j)20 continuereturnendsubroutine dot(m,n,l,a,b,c)! dimension a(m,n),b(n,l),c(m,l)do 50 i=1,mdo 50 j=1,lc(i,j)=0.0do 40 k=1,nc(i,j)=c(i,j)+a(i,k)*b(k,j)40 continue50 continuereturnend! c beam element(2-d) coordinate transformation subroutine cr(nn,ie,mea,x,y,r)dimension x(nn),y(nn),mea(100,4),r(6,6)call zero(6,6,r)i=mea(ie,1)j=mea(ie,2)sl=sqrt((x(j)-x(i))**2+(y(j)-y(i))**2)r(1,1)=(x(j)-x(i))/slr(1,2)=(y(j)-y(i))/slr(2,1)=-r(1,2)r(2,2)=r(1,1)r(3,3)=1.0r(4,4)=r(1,1)r(4,5)=r(1,2)r(5,4)=r(2,1)r(5,5)=r(2,2)r(6,6)=1.0returnend! c beam element(2-d)subroutine kep(nn,ie,mea,aeu,aai,x,y,akep)dimension x(nn),y(nn),mea(100,4),aeu(10,2),aai(10,2),akep(6,6) real izcall zero(6,6,akep)nm1=mea(ie,4)eo=aeu(nm1,1)na1=mea(ie,3)ao=aai(na1,1)iz=aai(na1,2)i=mea(ie,1)j=mea(ie,2)sl=sqrt((x(j)-x(i))**2+(y(j)-y(i))**2)akep(1,1)=eo*ao/slakep(4,4)=akep(1,1)akep(1,4)=-akep(1,1)akep(4,1)=akep(1,4)akep(2,2)=12.0*eo*iz/sl**3akep(2,5)=-akep(2,2)akep(5,2)=akep(2,5)akep(5,5)=akep(2,2)akep(2,3)=6.0*eo*iz/sl**2akep(3,2)=akep(2,3)akep(2,6)=akep(2,3)akep(6,2)=akep(2,6)akep(3,5)=-akep(2,3)akep(5,3)=akep(3,5)akep(5,6)=akep(5,3)akep(6,5)=akep(5,6)akep(3,3)=4.0*eo*iz/slakep(6,6)=akep(3,3)akep(3,6)=2.0*eo*iz/slakep(6,3)=akep(3,6)returnend! c fixed and forcesubroutine fixf(nn,ne,ie,mea,x,y,qq,pop)dimension x(nn),y(nn),qq(ne),mea(100,4),pop(6)i=mea(ie,1)j=mea(ie,2)sl=sqrt((x(j)-x(i))**2+(y(j)-y(i))**2)pop(1)=0.0pop(4)=0.0pop(2)=-qq(ie)*sl/2.0 pop(5)=pop(2)pop(3)=-qq(ie)*sl*sl/12.0 pop(6)=-pop(3)returnend4 3 1 1 51 0 0 0 0 02 3 0 0 -200 03 6 0 0 0 -504 12 0 0 0 01 11 21 33 24 21 12 1 1 02 23 1 1 03 34 1 1 -25.08.611E-3 2.17e-42.0E8 0.3。