線彈性小變形空間板殼靜力有限元計算程序

1、元計算有限元自動生成系統(tǒng)所開發(fā)源代碼系列線彈性小變形空間板殼靜力有限元計算程序簡介元計算（www.ectec.asia）公司所開發(fā)的并行有限元程序自動生成系統(tǒng)（pFEPG）可根據(jù)用戶需要開發(fā)出各種有限元計算程序源代碼。該源代碼系列即為pFEPG所開發(fā)出來的求解各學科典型問題的有限元計算程序。該組程序為線彈性小變形空間板殼靜力有限元計算程序。starta.for，對位移場的數(shù)據(jù)進行初始化； implicit real*8 (a-h,o-z) character*12 fname,filename(20) common /aa/ ia(250000000) common /bb/ ib(12500

2、0000)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=kvar+1 kcoor=3 kelem=31250

3、000 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=knb2+1 knb3=kvar1*1 if (k

4、nb3/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 55555 endif knb0=1 knb1=knb1+

5、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(kna0),ia(kna1),ia(kna2),

6、 *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),COOR(KCOOR,KNODE),R(3), * U

7、0(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.O.FC . U0 U1 U2 INITIAL

8、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(*,*) NUMNOD =,NUMNOD, NODDOF =,NODD

9、OFc 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=ORDER.NOD,FORM=UNFORMATTED,

10、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 DO 30 JNOD=1,NUMNOD J=INO

11、DVAR(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 =,NUMNOD, NODDOF =,NODDO

12、Fc 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=UNFORMATTED,STATUS=unkno

13、wn) 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=UNFORMATTED,STATUS=unknown) WRIT

14、E(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 PREPROCESSOR inquire(file=disp1,e

15、xist=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(J,N),J=1,NODDOF),N=1,NUMN

16、OD) 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 = C WRITE(*,(6F13.3) (U1(J

17、,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 THE INITIAL VALUE FILE U

18、NOD 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(21, (1a) material read(21,

19、*) 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 WRITE(*,*) NODE =cc WRITE(*,

20、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(*,*) NUM =,NUM, NNODE =,N

21、NODEcc 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 (*,*) L,LM =,LC WRITE (*,(1X,

22、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= ,FORM=UNFORMATTED,STATUS=un

23、known) 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) goto 800 j=id(k,n) j0=0 if

24、 (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 100 I=1,ND110 IF(LM(I)-LS)

25、 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 RETURN ENDeshell3da.for，G

26、alerkin法求解位移場的主程序 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= ,FORM=UNFORMATTED,STATUS=OL

27、D) 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(*,(1X,6I7) KNODE,KDGOF,KVAR,K

28、COOR,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=100000*2 knc5=neq*2 knc4=kd

29、gof*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 eshell3da stop 55555 endif knb0=1 knb1=knb1+knb0 knb2=knb2+knb1 if (knb2-1.gt.250000000)

30、then write(*,*) exceed memory of array ia write(*,*) memory of ia = 250000000 write(*,*) memory needed = ,knb2, in prgram eshell3da stop 55555 endif knc0=1 knc1=knc1+knc0 knc2=knc2+knc1 knc3=knc3+knc2 knc4=knc4+knc3 knc5=knc5+knc4 knc6=knc6+knc5 knc7=knc7+knc6 knc8=knc8+knc7 if (knc8-1.gt.62500000)

31、then write(*,*) exceed memory of array ic write(*,*) memory of ic = 62500000 write(*,*) memory needed = ,knc8, in prgram eshell3da stop 55555 endif call eshell3da(knode,kdgof,kvar,kcoor, *numtyp,numel,neq,kelem,kemate,maxa, *maxt,neq1,ib(kna0),ib(kna1),ib(kna2), *ia(knb0),ia(knb1),ic(knc0),ic(knc1),

32、ic(knc2), *ic(knc3),ic(knc4),ic(knc5),ic(knc6),ic(knc7), *filename) end subroutine eshell3da(knode,kdgof,kvar,kcoor, *numtyp,numel,neq,kelem,kemate,maxa, *maxt,neq1,nodvar,jdiag,node,a, *b,u,coor,f,ubf,u1, *emate,eu,sml, *filename) implicit real*8 (a-h,o-z) character*12 filename(20) DIMENSION NODVAR

33、(KDGOF,KNODE),U(KDGOF,KNODE),COOR(KCOOR,KNODE), *eu(kdgof,knode), & F(NEQ),A(MAXA),B(MAXA),JDIAG(NEQ),EMATE(KEMATE), & NODE(KELEM),SML(100000),u1(neq),UBF(KDGOF,KNODE)6 FORMAT (1X,15I5)7 FORMAT (1X,5e15.5)1001 FORMAT(1X,9I7)C.OPEN TIME File OPEN(1,FILE= ,FORM=UNFORMATTED,STATUS=OLD) READ(1) TMAX,DT,

34、TIME,IT WRITE(*,*) TMAX,DT,TIME,IT =,TMAX,DT,TIME,IT CLOSE(1)C.OPEN NODVAR file OPEN (1,FILE= ,FORM=UNFORMATTED,STATUS=OLD) READ (1) (NODVAR(I,J),I=1,KDGOF),J=1,KNODE) CLOSE (1)cc WRITE(*,*) KDGOF =,KDGOF, KNODE =,KNODEcc WRITE (*,*) NODVAR =cc WRITE (*,6) (NODVAR(I,J),I=1,KDGOF),J=1,KNODE)C.OPEN CO

35、OR file OPEN (1,FILE= ,FORM=UNFORMATTED,STATUS=OLD) READ (1) NUMNOD,NCOOR,(COOR(I,J),I=1,NCOOR),J=1,NUMNOD) CLOSE(1)cc WRITE(*,*) NUMNOD,NCOOR=,NUMNOD,NCOORC.OPEN BFD file OPEN (1,FILE= ,FORM=UNFORMATTED,STATUS=OLD) READ (1) (UBF(J,I),J=1,KDGOF),I=1,KNODE) CLOSE (1)cc WRITE (*,*) BF =cc WRITE(*,7) (

36、UBF(J,I),J=1,KDGOF),I=1,KNODE) numtyp = 1C.OPEN DIAG file OPEN (2,FILE= ,FORM=UNFORMATTED,STATUS=OLD) READ(2) (JDIAG(I),I=1,NEQ) CLOSE(2)C.OPEN ELEM0 file OPEN (3,FILE= ,FORM=UNFORMATTED,STATUS=OLD) itime=01 continue itime=itime+1 if (itime.gt.1) then write(*,*) Nonlinear Iteration Times =,itime rew

37、ind(3) endif DO 111 I=1,KNODE DO 111 J=1,KDGOF U(J,I) = UBF(J,I)111 CONTINUEcc WRITE (*,*) BF =cc WRITE(*,7) (U(J,I),J=1,KDGOF),I=1,KNODE) DO 112 I=1,MAXA A(I) = 0.0 B(I) = 0.0112 CONTINUE DO 2300 I=1,NEQ2300 CONTINUE NUMEL=0C.OPEN EMATE+ENODE+ELOAD fileC OPEN (3,FILE= ,FORM=UNFORMATTED,STATUS=OLD)

38、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 WRITE(*,*) NODE =cc WRITE(*,6) (NODE(I-1)*NNODE+J),J=1,NNODE),I=1,NUM) NNE = NNODE nne = nne-1 K=0 DO 115 J=1,NNE JNOD = NODE(J) DO 115 L=1,KDGOF IF (NODVAR(L,JNOD

39、).NE.0) K=K+1115 CONTINUE WRITE(*,*) K =,K kk=k*k k0=1 k1=k0+k*k k2=k1+k k3=k2+k k4=k3+k*k k5=k4+k*k CALL ETSUB(KNODE,KDGOF,IT,KCOOR,KELEM,K,KK,NNODE,NNE, * ITYP,NCOOR,NUM,TIME,DT,neq,maxa,NODVAR,COOR,NODE,EMATE, & A,B,JDIAG, &sml(k0),sml(k1),sml(k2),sml(k3),sml(k4), &eu, *U)2000 CONTINUE DO 2050 IJ

40、=1,NEQ if (itime.le.1) u1(IJ) = 0.0 F(IJ)=0.0D02050 CONTINUE DO 2200 I=1,KNODE DO 2100 J=1,KDGOF IJ=NODVAR(J,I) IF (IJ.LE.0) GOTO 2100 F(IJ)=F(IJ)+U(J,I) U1(IJ)=F(IJ)2100 CONTINUE2200 CONTINUECC IF (IT.GT.0) THENcc WRITE (*,*) U =cc WRITE (*,7) (U(J,I),J=1,KDGOF),I=1,KNODE)cc WRITE (*,*) NEQ =,NEQ,

41、F =cc WRITE(*,7) (F(I),I=1,NEQ) if (itime.le.1) thenC.OPEN LMATRIX FILE OPEN (2,FILE= ,FORM=UNFORMATTED,STATUS=unknown) CLOSE (2) endif WRITE(*,*) NIN_SOLVER MEMORY REQUIRED . ,MAXA IF (MAXA.GT.MAXT) THEN WRITE(*,*) WARNING MATRIX A EXCEED CORE MEMORY . ,MAXTc STOP 0000 ENDIF CALL REDU(A,B,U1,JDIAG,

42、NEQ,MAXA,1)C WRITE(*,*) U1 = C WRITE(*,7) (A(I),I,MAXA)C WRITE(*,7) (F(I),I=1,NEQ) NOUT = 20 OPEN(NOUT,FILE= ,FORM=FORMATTED,STATUS=unknown) DO 3200 INOD=1,KNODE DO 3100 IDFG=1,KDGOF N=NODVAR(IDFG,INOD)C WRITE (*,*) N =,N if(n.le.0) then eu(IDFG,INOD)=u(IDFG,INOD) else eu(IDFG,INOD)=u1(N) endif3100

43、CONTINUE3200 CONTINUE DO 3400 N=1,KNODE WRITE (NOUT,3600) N,(eu(I,N),I=1,KDGOF)3400 CONTINUE3600 FORMAT (1X,I5,1X,6E11.4,9(/6X,6E11.4) CLOSE (NOUT) open(10,file=unod,form=unformatted,status=unknown) write(10) (eu(j,i),i=1,knode),j=1,kdgof) close(10) close (3) RETURN END SUBROUTINE ETSUB(KNODE,KDGOF,

44、IT,KCOOR,KELEM,K,KK,NNODE,NNE, *ITYP,NCOOR,NUM,TIME,DT,neq,maxa,NODVAR,COOR,NODE,EMATE, &A,B,JDIAG, *es,em,ef,Estifn,Estifv,eu, *U) implicit real*8 (a-h,o-z) DIMENSION NODVAR(KDGOF,KNODE),COOR(KCOOR,KNODE),NODE(KELEM), *U(KDGOF,KNODE),EMATE(300), &A(MAXa),B(MAXa),JDIAG(neq), *es(k,k),em(k),ef(k),eu(

45、kdgof,knode), *Estifn(k,k),Estifv(kk), *R(500),PRMT(500),COEF(500),LM(500)17 FORMAT (1X,15I5)18 FORMAT (1X,8e9.2) READ (3) MMATE,NMATE,(EMATE(I-1)*NMATE+J),J=1,NMATE), * I=1,MMATE) WRITE(*,*) MMATE =,MMATE, NMATE =,NMATE WRITE (*,*) EMATE = WRITE (*,18) (EMATE(I-1)*NMATE+J),J=1,NMATE), * I=1,MMATE)

46、DO 1000 NE=1,NUM NR=0 DO 130 J=1,NNE JNOD = NODE(NE-1)*NNODE+J) IF (JNOD.LT.0) JNOD = -JNOD PRMT(NMATE+7+J) = JNOD DO 120 I=1,NCOOR NR=NR+1120 R(NR) = COOR(I,JNOD)130 CONTINUE IMATE = NODE(NNODE*NE) DO 140 J=1,NMATE140 PRMT(J) = EMATE(IMATE-1)*NMATE+J) PRMT(NMATE+1)=TIME PRMT(NMATE+2)=DT PRMT(NMATE+

47、3)=IMATE prmt(NMATE+4)=NE prmt(NMATE+5)=NUM prmt(NMATE+6)=IT prmt(NMATE+7)=NMATE prmt(NMATE+8)=ITIME prmt(NMATE+9)=ITYP goto 11 call csugt3m(r,coef,prmt,es,em,ec,ef,ne) goto 22 continueC WRITE(*,*) ES EM EF =C DO 555 I=1,KC555 WRITE(*,18) (ES(I,J),J=1,K)C WRITE(*,18) (EM(I),I=1,K)C WRITE(*,18) (EF(I

48、),I=1,K)CC IF (IT.GT.0) THEN do 201 i=1,k do 201 j=1,k Estifn(i,j)=0.0201 continue do 202 i=1,k Estifn(i,i)=Estifn(i,i) do 202 j=1,k Estifn(i,j)=Estifn(i,j)+es(i,j)202 continue L=0 M=0 I=0 DO 700 INOD=1,NNE NODI=NODE(NE-1)*NNODE+INOD) DO 600 IDGF=1,KDGOF INV=NODVAR(IDGF,NODI) IF (INV.EQ.0) GOTO 600

49、I=I+1 IF (INV.LT.0) GOTO 305 L=L+1 LM(L)=INV U(IDGF,NODI)=U(IDGF,NODI) *+ef(i)305 J=0 DO 500 JNOD=1,NNE NODJ=NODE(NE-1)*NNODE+JNOD) DO 400 JDGF=1,KDGOF JNV=NODVAR(JDGF,NODJ) IF (JNV.EQ.0) GOTO 400 J=J+1 IF (JNV.LT.0) GOTO 400 IF (INV.LT.0) GOTO 310 M=M+1 Estifv(m)=Estifn(i,j)310 CONTINUE IF (INV.LT.

50、0) * U(JDGF,NODJ)=U(JDGF,NODJ)-ESTIFN(I,J)*U(IDGF,NODI)400 CONTINUE500 CONTINUE600 CONTINUE700 CONTINUEC WRITE (*,*) U =C WRITE (*,18) (U(J,I),J=1,KDGOF),I=1,KNODE) LRD=M NER=NUMEL+NEC WRITE(*,*) *C WRITE(*,*) (ESTIFV(I),I=1,LRD)C WRITE (*,*) Einform .C WRITE (*,(1X,15I5) L,LRD,(LM(I),I=1,L) DO 800

51、I=1,L J=LM(I)800 CONTINUE call ADDA(A,B,JDIAG,L,LM,ESTIFV,NEQ,MAXA)1000 CONTINUE RETURN END SUBROUTINE ADDA(A,B,JDIAG,ND,LM,ESTIF,NEQ,MAXA) implicit real*8 (a-h,o-z) DIMENSION A(MAXA),B(MAXA),JDIAG(NEQ),LM(ND),ESTIF(ND,ND)C WRITE (*,*) ND, (LM(I),I=1,ND)C WRITE (*,*) (ESTIF(I,J),J=1,ND),I=1,ND) DO 3

52、00 I=1,ND II = LM(I) DO 280 J=1,I JJ = LM(J) IF (II.LT.JJ) GOTO 240 K = JDIAG(II) - II + JJ A(K) = A(K) + ESTIF(I,J) B(K) = B(K) + ESTIF(J,I) GOTO 280240 K = JDIAG(JJ) - JJ + II B(K) = B(K) + ESTIF(I,J) A(K) = A(K) + ESTIF(J,I)280 CONTINUE300 CONTINUE RETURN END SUBROUTINE REDU(A,B,U,JDIAG,NEQ,MAXA,

53、KKK) implicit real*8 (a-h,o-z) DIMENSION A(MAXA),B(MAXA),JDIAG(NEQ),U(NEQ) DOUBLE PRECISION CC.COMPUTE L U & u, L*U(T) = A, A*u = f DO 500 I=2,NEQ I1 = I-1 NI = JDIAG(I) LI = I-NI+JDIAG(I-1)+1 IF (KKK.GT.1) GOTO 333 DO 200 J=LI,I NJ = JDIAG(J) LJ = J-NJ+1 IF (J.GT.1) LJ = LJ+JDIAG(J-1) LIJ = MAX0(LI

54、,LJ) J1 = J-1 IF (J.EQ.I) GOTO 130 C = 0.0D0 DO 100 L=LIJ,J1100 C = C + A(NI-I+L)*B(NJ-J+L) A(NI-I+J) = (A(NI-I+J) - C)/B(NJ)130 C = 0.0D0 DO 150 L=LIJ,J1150 C = C + B(NI-I+L)*A(NJ-J+L)200 B(NI-I+J) = B(NI-I+J) - C333 C = 0.0D0 DO 400 L=LI,I1400 C = C + A(NI-I+L)*U(L) U(I) = U(I)-C500 CONTINUE J = M

55、AXA + 1 N = NEQ + 1700 N = N - 1 J = J - 1 U(N) = U(N)/B(J) IF (N.EQ.1) GOTO 999 M = J-JDIAG(N-1)-1 DO 800 I=1,M J = J - 1800 U(N-I) = U(N-I) - B(J)*U(N) GOTO 700999 RETURN ENDcsult3m.for，計算板殼單元剛度矩陣和荷載向量的子程序 subroutine csult3m(coorr,coefr, & prmt,estif,emass,edamp,eload,num)c . coorr - nodal coordin

56、ate valuec . coefr - nodal coef value implicit real*8 (a-h,o-z) common /csult3mcom1/ det,d1,d2,d3,b1,b2,b3 dimension estif(18,18),elump(18),emass(18), & eload(18) dimension prmt(*), & eekx(18),eeky(18),eekxy(18),egmx(18), & egmy(18),eex(18),eey(18),eexy(18), & coorr(2,3),coor(2) common /rcsult3m/ru(

57、3,18),rv(3,18),rw(9,18), & rs(9,18),ro(9,18),rc(3,18), & cu(3,3),cv(3,3),cw(9,3),cs(9,3), & co(9,3),cc(3,3)c . store shape functions and their partial derivativesc . for all integral points common /vcsult3m/rctr(2,2),crtr(2,2) common /dcsult3m/ refc(2,6),gaus(6), & nnode,ngaus,ndisp,nrefc,ncoor,nvar

58、, & nvard(6),kdord(6),kvord(18,6)c . nnode - the number of nodesc . nrefc - the number of numerical integral pointsc . ndisp - the number of unknown functionsc . nrefc - the number of reference coordinatesc . nvar - the number of unknown varibles varc . refc - reference coordinates at integral point

59、sc . gaus - weight number at integral pointsc . nvard - the number of var for each unknownc . kdord - the highest differential order for each unknownc . kvord - var number at integral points for each unknown pe=prmt(1) pv=prmt(2) thick=prmt(3) fu=prmt(4) fv=prmt(5) fw=prmt(6) rou=prmt(7) alpha=prmt(

60、8) time=prmt(9) dt=prmt(10) imate=prmt(11)+0.5 ielem=prmt(12)+0.5 nelem=prmt(13)+0.5 it=prmt(14)+0.5 nmate=prmt(15)+0.5 itime=prmt(16)+0.5 ityp=prmt(17)+0.5 fact = pe*thick*3/12./(1.-pv*pv) fact2 = pe/(1.+pv)/(1.-pv)*thick b1=coorr(2,2)-coorr(2,3) b2=coorr(2,3)-coorr(2,1) b3=coorr(2,1)-coorr(2,2) d1