线弹性小变形平面刚架静力有限元计算程序.docx

1、线弹性小变形平面刚架静力有限元计算程序元计算有限元自动生成系统所开发源代码系列线弹性小变形平面刚架静力有限元计算程序1. 简介元计算（www.ectec.asia）公司所开发的并行有限元程序自动生成系统（pFEPG）可根据用户需要开发出各种有限元计算程序源代码。该源代码系列即为pFEPG所开发出来的求解各学科典型问题的有限元计算程序。该组程序为线弹性小变形平面刚架静力有限元计算程序。2. starta.for，对位移场的数据进行初始化； implicit real*8 (a-h,o-z) character*12 fname,filename(20) common /aa/ ia(250000

2、000) common /bb/ ib(125000000)c. open disp0 file to get the numbers of nodes and degree of freedomc. knode . number of nodes, kdgof . number of d.o.f. open(1,file= ,form=unformatted) read(1) knode,kdgof close(1) kvar=knode*kdgof write(*,*) knode,kdgof,kvar = write(*,(1x,4i7) knode,kdgof,kvar kvar1=k

3、var+1 kcoor=3 kelem=31250000 knb1=kdgof*knode*1 if (knb1/2*2 .lt. knb1) knb1=knb1+1 kna4=kcoor*knode*2 kna1=kdgof*knode*2 kna2=kdgof*knode*2 kna3=kdgof*knode*2 kna5=knode*1 if (kna5/2*2 .lt. kna5) kna5=kna5+1 knb4=kelem*1 if (knb4/2*2 .lt. knb4) knb4=knb4+1 knb2=kvar1*1 if (knb2/2*2 .lt. knb2) knb2=

4、knb2+1 knb3=kvar1*1 if (knb3/2*2 .lt. knb3) knb3=knb3+1 kna0=1 kna1=kna1+kna0 kna2=kna2+kna1 kna3=kna3+kna2 kna4=kna4+kna3 kna5=kna5+kna4 if (kna5-1.gt.250000000) then write(*,*) exceed memory of array ia write(*,*) memory of ia = 250000000 write(*,*) memory needed = ,kna5, in prgram start stop 5555

5、5 endif knb0=1 knb1=knb1+knb0 knb2=knb2+knb1 knb3=knb3+knb2 knb4=knb4+knb3 if (knb4-1.gt.125000000) then write(*,*) exceed memory of array ib write(*,*) memory of ib = 125000000 write(*,*) memory needed = ,knb4, in prgram start stop 55555 endif call start(knode,kdgof,kcoor,kvar, *kelem,maxt,kvar1,ia

6、(kna0),ia(kna1),ia(kna2), *ia(kna3),ia(kna4),ib(knb0),ib(knb1),ib(knb2), *ib(knb3), *filename) end subroutine start(knode,kdgof,kcoor,kvar, *kelem,maxt,kvar1,u0,u1,u2, *coor,inodvar,nodvar,numcol,lm,node, *filename) implicit real*8 (a-h,o-z) character*12 filename(20) DIMENSION NODVAR(KDGOF,KNODE),CO

7、OR(KCOOR,KNODE),R(3), * U0(KDGOF,KNODE),U1(KDGOF,KNODE),U2(KDGOF,KNODE), * INODVAR(KNODE),node(kelem) DIMENSION NUMCOL(KVAR1),LM(KVAR1) CHARACTER*1 MATERIAL logical filflgC .C . KDGOF NUMBER OF D.O.FC . KNODE NUMBER OF NODESC . INODVAR ID DATAC . NODVAR DENOTE THE EQUATION NUMBER CORRESPONDING THE D

8、.O.FC . U0 U1 U2 INITIAL VALUEC . COOR COORDINATESC . NODE ELEMENT NODAL CONNECTIONC .6 FORMAT (1X, 15I4)7 FORMAT (1X,8F9.3)C.OPEN ID file OPEN (1,FILE= ,FORM=UNFORMATTED,STATUS=OLD) READ (1) NUMNOD,NODDOF,(NODVAR(I,J),I=1,NODDOF),J=1,NUMNOD) CLOSE (1) call chms(kdgof,knode,NODVAR)c WRITE(*,*) NUMNO

9、D =,NUMNOD, NODDOF =,NODDOFc WRITE (*,*) ID =c WRITE (*,6) (NODVAR(I,J),I=1,NODDOF),J=1,NUMNOD)C. GET THE NATURAL NODAL ORDER DO 12 N=1,KNODE INODVAR(N)=N12 CONTINUEC. OPEN ORDER.NOD FILE AND READ THE NODAL ORDER IF THE FILE EXIST inquire(file=ORDER.NOD,exist=filflg) if (filflg) then OPEN (1,FILE=OR

10、DER.NOD,FORM=UNFORMATTED,STATUS=OLD) READ (1) (INODVAR(I),I=1,NUMNOD) CLOSE(1) WRITE(*,*) NODORDER = WRITE(*,6) (INODVAR(I),I=1,NUMNOD) endifC. GET NV BY ID NEQ=0 DO 20 JNOD=1,NUMNOD J=INODVAR(JNOD) DO 18 I=1,NODDOF IF (NODVAR(I,J).NE.1) GOTO 18 NEQ = NEQ + 1 NODVAR(I,J) = NEQ18 CONTINUE20 CONTINUE

11、DO 30 JNOD=1,NUMNOD J=INODVAR(JNOD) DO 28 I=1,NODDOF IF (NODVAR(I,J).GE.-1) GOTO 28 N = -NODVAR(I,J)-1 NODVAR(I,J) = NODVAR(I,N)28 CONTINUE30 CONTINUEC. OPEN AND WRITE THE NV FILE OPEN(8,STATUS=unknown,FILE= ,FORM=UNFORMATTED) WRITE(8) (NODVAR(I,J),I=1,NODDOF),J=1,NUMNOD) CLOSE(8)c WRITE(*,*) NUMNOD

12、 =,NUMNOD, NODDOF =,NODDOFc WRITE(*,6) (NODVAR(I,J),I=1,NODDOF),J=1,NUMNOD)C. WRITE THE BOUNDAY CONDITION FILE BFD ACCORDING TO THE DISP0 FILEC.OPEN DISP0 FILE OPEN(1,FILE= ,FORM=UNFORMATTED,STATUS=OLD) READ(1) NUMNOD,NODDOF,(U0(I,J),I=1,NODDOF),J=1,NUMNOD) CLOSE(1)C.OPEN BFD FILE OPEN(1,FILE= ,FORM

13、=UNFORMATTED,STATUS=unknown) WRITE(1) (U0(I,J),I=1,NODDOF),J=1,NUMNOD) CLOSE(1)C. GET THE INITIAL TIME FROM TIME0 FILEC.OPEN TIME0 File OPEN(1,FILE= ,FORM=FORMATTED) READ(1,*) T0,TMAX,DT TIME = T0 IT = 0 WRITE(*,*) TMAX,DT,TIME,IT =,TMAX,DT,TIME,IT CLOSE(1)C.OPEN TIME File OPEN(1,FILE= ,FORM=UNFORMA

14、TTED,STATUS=unknown) WRITE(1) TMAX,DT,TIME,IT CLOSE(1)C.OPEN COOR file OPEN (1,FILE= ,FORM=UNFORMATTED,STATUS=OLD) READ (1) NUMNOD,NCOOR,(COOR(I,J),I=1,NCOOR),J=1,NUMNOD) CLOSE(1)c WRITE(*,*) COOR =c WRITE(*,7) (COOR(I,J),I=1,NCOOR),J=1,NUMNOD)C. GET THE INITIAL VALUE FROM THE DATA FILES BY PREPROCE

15、SSOR inquire(file=disp1,exist=filflg) if (filflg) then open(16,file=disp1,form=unformatted,status=old) read(16) numnod,noddof,(U0(J,N),J=1,NODDOF),N=1,NUMNOD) close(16) endif inquire(file=disp2,exist=filflg) if (filflg) then open(16,file=disp2,form=unformatted,status=old) read(16) numnod,noddof,(U1(

16、J,N),J=1,NODDOF),N=1,NUMNOD) close(16) endif inquire(file=disp3,exist=filflg) if (filflg) then open(16,file=disp3,form=unformatted,status=old) read(16) numnod,noddof,(U2(J,N),J=1,NODDOF),N=1,NUMNOD) close(16) endifc WRITE(*,*) U0 = c WRITE(*,(6F13.3) (U0(J,N),J=1,NODDOF),N=1,NUMNOD)C WRITE(*,*) U1 =

17、 C WRITE(*,(6F13.3) (U1(J,N),J=1,NODDOF),N=1,NUMNOD)C. COMPUTE THE INITIAL VALUE BY BOUND.FOR zo = 0.0d0c DO 321 N=1,NUMNODc DO 100 J=1,NCOORc100 R(J) = COOR(J,N)c DO 200 J=1,NODDOFc U0(J,N) = BOUND(R,zo,J)c U1(J,N) = BOUND1(R,zo,J)c U2(J,N) = BOUND2(R,zo,J)c200 CONTINUEc321 CONTINUEC.OPEN AND WRITE

18、 THE INITIAL VALUE FILE UNOD OPEN (1,FILE= ,FORM=UNFORMATTED,STATUS=unknown) WRITE(1) (U0(I,J),J=1,NUMNOD),I=1,NODDOF), * (U1(I,J),J=1,NUMNOD),I=1,NODDOF), * (U2(I,J),J=1,NUMNOD),I=1,NODDOF), * (U0(I,J),J=1,NUMNOD),I=1,NODDOF) CLOSE (1)c. open IO file open(21,file= ,form=formatted,status=old) read(2

19、1, (1a) material read(21,*) numtyp close(21) DO I=1,NEQ NUMCOL(i)=1 ENDDOC.OPEN ELEM0 file OPEN (3,FILE= ,FORM=UNFORMATTED,STATUS=OLD) NUMEL=0 KELEM=0 KEMATE=0 DO 2000 ITYP=1,NUMTYPC.INPUT ENODE READ (3) NUM,NNODE, * (NODE(I-1)*NNODE+J),J=1,NNODE),I=1,NUM)cc WRITE(*,*) NUM =,NUM, NNODE =,NNODEcc WRI

20、TE(*,*) NODE =cc WRITE(*,6) (NODE(I-1)*NNODE+J),J=1,NNODE),I=1,NUM) IF (KELEM.LT.NUM*NNODE) KELEM = NUM*NNODE NNE = NNODE IF (MATERIAL.EQ.Y .OR. MATERIAL.EQ.y) THEN READ (3) MMATE,NMATE IF (KEMATE.LT.MMATE*NMATE) KEMATE = MMATE*NMATE NNE = NNE-1 ENDIF WRITE(*,*) MMATE =,MMATE, NMATE =,NMATEcc WRITE(

21、*,*) NUM =,NUM, NNODE =,NNODEcc WRITE(*,*) NODE =cc WRITE(*,6) (NODE(I-1)*NNODE+J),J=1,NNODE),I=1,NUM) DO 1000 NE=1,NUM L=0 DO 700 INOD=1,NNE NODI=NODE(NE-1)*NNODE+INOD) DO 600 IDGF=1,KDGOF INV=NODVAR(IDGF,NODI) IF (INV.LE.0) GOTO 600 L=L+1 LM(L)=INV600 CONTINUE700 CONTINUE NUMEL=NUMEL+1C WRITE (*,*

22、) L,LM =,LC WRITE (*,(1X,15I5) (LM(I),I=1,L) if (l.gt.0) call ACLH(NEQ,NUMCOL,l,lm)1000 continue2000 CONTINUEc CLOSE(1) CLOSE(3) call BCLH(NEQ,NUMCOL) MAXA=NUMCOL(NEQ)C.OPEN SYS File OPEN (2,FILE= ,FORM=UNFORMATTED,STATUS=unknown) WRITE(2) NUMEL,NEQ,NUMTYP,MAXA,KELEM,KEMATE CLOSE (2) OPEN(2,FILE= ,F

23、ORM=UNFORMATTED,STATUS=unknown) write(2) (NUMCOL(I),I=1,NEQ) CLOSE(2)c write(*,*) NEQ,NUMCOL=,NEQc write(*,6) (NUMCOL(i),i=1,NEQ) END subroutine chms(kdgof,knode,id) dimension id(kdgof,knode),ms(1000),is(1000) do 1000 k=1,kdgof m = 0 do 800 n=1,knode if (id(k,n).le.-1) id(k,n)=-1 if (id(k,n).le.1) g

24、oto 800 j=id(k,n) j0=0 if (m.gt.0) then do i=1,m if (j.eq.ms(i) j0=is(i) enddo endif if (j0.eq.0) then m=m+1 ms(m)=j is(m)=n id(k,n)=1 else id(k,n)=-j0-1 endif800 continue1000 continue return end SUBROUTINE ACLH(NEQ,NUMCOL,ND,LM) implicit real*8 (a-h,o-z) DIMENSION LM(ND),NUMCOL(NEQ) LS=LM(1)+1 DO 1

25、00 I=1,ND110 IF(LM(I)-LS) 120,100,100120 LS=LM(I)100 CONTINUE DO 200 I=1,ND II=LM(I) ME=II-LS IF(ME.GT.NUMCOL(II) NUMCOL(II)=ME200 CONTINUE RETURN END SUBROUTINE BCLH(NEQ,NUMCOL) implicit real*8 (a-h,o-z) DIMENSION NUMCOL(NEQ)C NUMCOL(1) = 1 DO 490 I=2,NEQ490 NUMCOL(I) = NUMCOL(I) + NUMCOL(I-1) + 1

26、RETURN END3. ebeam2da.for，Galerkin法求解位移场的主程序 implicit real*8 (a-h,o-z) character*12 fname,filename(20) common /aa/ ia(250000000) common /bb/ ib(125000000) common /cc/ ic(62500000) open(1,file= ,form=unformatted,status=old) read(1) knode,kdgof close(1) MAXT=250000000/2/2C.OPEN SYS File OPEN (2,FILE=

27、,FORM=UNFORMATTED,STATUS=OLD) read(2) NUMEL,NEQ,NUMTYP,MAXA,KELEM,KEMATE CLOSE (2) IF (MAXA.GT.MAXT) THEN WRITE(*,*) MATRIX A EXCEED CORE MEMERY . ,MAXA WRITE(*,*) REQUIRED CORE MEMERY . ,MAXT STOP 0000 ENDIF KVAR=KNODE*KDGOF KCOOR=3C KELEM=31250000 WRITE(*,*) KNODE,KDGOF,KVAR,KCOOR,KELEM = WRITE(*,

28、(1X,6I7) KNODE,KDGOF,KVAR,KCOOR,KELEM kna1=kdgof*knode*1 if (kna1/2*2 .lt. kna1) kna1=kna1+1 knc1=kdgof*knode*2 knc2=kcoor*knode*2 knc7=kdgof*knode*2 knc3=neq*2 knb1=maxa*2 knb2=maxa*2 kna2=neq*1 if (kna2/2*2 .lt. kna2) kna2=kna2+1 knc6=kemate*2 kna3=kelem*1 if (kna3/2*2 .lt. kna3) kna3=kna3+1 knc8=

29、100000*2 knc5=neq*2 knc4=kdgof*knode*2 kna0=1 kna1=kna1+kna0 kna2=kna2+kna1 kna3=kna3+kna2 if (kna3-1.gt.125000000) then write(*,*) exceed memory of array ib write(*,*) memory of ib = 125000000 write(*,*) memory needed = ,kna3, in prgram ebeam2da stop 55555 endif knb0=1 knb1=knb1+knb0 knb2=knb2+knb1 if (knb2-1.gt.250000000) then write(*,*) excee