各向同性线弹性小变形平面应力静力有限元计算程序

元计算有限元自动生成系统所开发源代码系列 各向同性线弹性小变形平面应力静力有限元计算程序 1. 简介 元计算（www.ectec.asia）公司所开发的并行有限元程序自动生成系统（pFEPG）可根据用户需要开发出各 种有限元计算程序源代码。该源代码系列即为 pFEPG 所开发出来的求解各学科典型问题的有限元计算程序。 该组程序为各向同性线弹性小变形平面应力静力有限元计算程序。 2. starta.for，对位移场的数据进行初始化； implicit real*8 (a-h,o-z) character*12 fname,filename(20) common /aa/ ia(250000000) common /bb/ ib(125000000) c.... open disp0 file to get the numbers of nodes and degree of freedom c.... 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=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=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 55555 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(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),COOR(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 filflg C ................................................................. C ..... KDGOF NUMBER OF D.O.F C ..... KNODE NUMBER OF NODES C ..... INODVAR ID DATA C ..... NODVAR DENOTE THE EQUATION NUMBER CORRESPONDING THE D.O.F C ..... U0 U1 U2 INITIAL VALUE C ..... COOR COORDINATES C ..... NODE ELEMENT NODAL CONNECTION C ................................................................. 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 =',NODDOF c 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)=N 12 CONTINUE C..... 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',STATUS='OLD') READ (1) (INODVAR(I),I=1,NUMNOD) CLOSE(1) WRITE(*,*) 'NODORDER =' WRITE(*,6) (INODVAR(I),I=1,NUMNOD) endif C..... 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) = NEQ 18 CONTINUE 20 CONTINUE 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 CONTINUE 30 CONTINUE C..... 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 =',NODDOF c WRITE(*,6) ((NODVAR(I,J),I=1,NODDOF),J=1,NUMNOD) C.... WRITE THE BOUNDAY CONDITION FILE BFD ACCORDING TO THE DISP0 FILE C....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='unknown') WRITE(1) ((U0(I,J),I=1,NODDOF),J=1,NUMNOD) CLOSE(1) C...... GET THE INITIAL TIME FROM TIME0 FILE C.......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') 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 PREPROCESSOR 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(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) endif c WRITE(*,*) ' U0 = ' c WRITE(*,'(6F13.3)') ((U0(J,N),J=1,NODDOF),N=1,NUMNOD) C WRITE(*,*) ' U1 = ' C WRITE(*,'(6F13.3)') ((U1(J,N),J=1,NODDOF),N=1,NUMNOD) C...... COMPUTE THE INITIAL VALUE BY BOUND.FOR zo = 0.0d0 c DO 321 N=1,NUMNOD c DO 100 J=1,NCOOR c100 R(J) = COOR(J,N) c DO 200 J=1,NODDOF c U0(J,N) = BOUND(R,zo,J) c U1(J,N) = BOUND1(R,zo,J) c U2(J,N) = BOUND2(R,zo,J) c200 CONTINUE c321 CONTINUE C.......OPEN AND WRITE 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(21, '(1a)') material read(21,*) numtyp close(21) DO I=1,NEQ NUMCOL(i)=1 ENDDO C.......OPEN ELEM0 file OPEN (3,FILE=' ',FORM='UNFORMATTED',STATUS='OLD') NUMEL=0 KELEM=0 KEMATE=0 DO 2000 ITYP=1,NUMTYP C.......INPUT ENODE READ (3) NUM,NNODE, * ((NODE((I-1)*NNODE+J),J=1,NNODE),I=1,NUM) cc WRITE(*,*) 'NUM =',NUM,' NNODE =',NNODE cc WRITE(*,*) '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 =',NMATE cc WRITE(*,*) 'NUM =',NUM,' NNODE =',NNODE cc 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)=INV 600 CONTINUE 700 CONTINUE NUMEL=NUMEL+1 C WRITE (*,*) 'L,LM =',L C WRITE (*,'(1X,15I5)') (LM(I),I=1,L) if (l.gt.0) call ACLH(NEQ,NUMCOL,l,lm) 1000 continue 2000 CONTINUE c 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='unknown') write(2) (NUMCOL(I),I=1,NEQ) CLOSE(2) c write(*,*) 'NEQ,NUMCOL=',NEQ c 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 (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 endif 800 continue 1000 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,ND 110 IF(LM(I)-LS) 120,100,100 120 LS=LM(I) 100 CONTINUE DO 200 I=1,ND II=LM(I) ME=II-LS IF(ME.GT.NUMCOL(II)) NUMCOL(II)=ME 200 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,NEQ 490 NUMCOL(I) = NUMCOL(I) + NUMCOL(I-1) + 1 RETURN END 3. esolid2dda.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 C.......OPEN SYS File OPEN (2,FILE=' ',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=3 C KELEM=31250000 WRITE(*,*) 'KNODE,KDGOF,KVAR,KCOOR,KELEM =' WRITE(*,'(1X,6I7)') KNODE,KDGOF,KVAR,KCOOR,KELEM kna1=kdgof*knode*1 if (kna1/2*2 .lt. kna1) kna1=kna1+1 knc3=kdgof*knode*2 knc1=kcoor*knode*2 knc7=kdgof*knode*2 knc2=neq*2 knb1=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=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 esolid2dda' stop 55555 endif knb0=1 knb1=knb1+knb0 if (knb1-1.gt.250000000) then write(*,*) 'exceed memory of array ia' write(*,*) 'memory of ia = 250000000' write(*,*) 'memory needed = ',knb1,' in prgram esolid2dda' 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) then write(*,*) 'exceed memory of array ic' write(*,*) 'memory of ic = 62500000' write(*,*) 'memory needed = ',knc8,' in prgram esolid2dda' stop 55555 endif call esolid2dda(knode,kdgof,kvar,kcoor, *numtyp,numel,neq,kelem,kemate,maxa, *maxt,neq1,ib(kna0),ib(kna1),ib(kna2), *ia(knb0),ic(knc0),ic(knc1),ic(knc2),ic(knc3), *ic(knc4),ic(knc5),ic(knc6),ic(knc7), *filename) end subroutine esolid2dda(knode,kdgof,kvar,kcoor, *numtyp,numel,neq,kelem,kemate,maxa, *maxt,neq1,nodvar,jdiag,node,a, *coor,f,u,ubf,u1,emate, *eu,sml, *filename) implicit real*8 (a-h,o-z) character*12 filename(20) DIMENSION NODVAR(KDGOF,KNODE),U(KDGOF,KNODE),COOR(KCOOR,KNODE), *eu(kdgof,knode), & F(NEQ),A(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,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 =',KNODE cc WRITE (*,*) 'NODVAR =' cc WRITE (*,6) ((NODVAR(I,J),I=1,KDGOF),J=1,KNODE) 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) cc WRITE(*,*) 'NUMNOD,NCOOR=',NUMNOD,NCOOR C.......OPEN BF 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) ((U(J,I),J=1,KDGOF),I=1,KNODE) numtyp = 2 C.......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=0 1 continue itime=itime+1 if (itime.gt.1) then write(*,*) 'Nonlinear Iteration Times ========',itime rewind(3) endif DO 111 I=1,KNODE DO 111 J=1,KDGOF U(J,I) = UBF(J,I) 111 CONTINUE cc WRITE (*,*) 'BF =' cc WRITE(*,7) ((U(J,I),J=1,KDGOF),I=1,KNODE) DO 112 I=1,MAXA A(I) = 0.0 112 CONTINUE DO 2300 I=1,NEQ 2300 CONTINUE NUMEL=0 C.......OPEN EMATE+ENODE+ELOAD file C OPEN (3,FILE=' ',FORM='UNFORMATTED',STATUS='OLD') DO 2000 ITYP=1,NUMTYP C.......INPUT ENODE READ (3) NUM,NNODE, * ((NODE((I-1)*NNODE+J),J=1,NNODE),I=1,NUM) cc WRITE(*,*) 'NUM =',NUM,' NNODE =',NNODE cc 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).NE.0) K=K+1 115 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,JDIAG, &sml(k0),sml(k1),sml(k2),sml(k3),sml(k4), &eu, *U) 2000 CONTINUE DO 2050 IJ=1,NEQ if (itime.le.1) u1(IJ) = 0.0 F(IJ)=0.0D0 2050 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 CONTINUE 2200 CONTINUE cc WRITE (*,*) 'U =' cc WRITE (*,7) ((U(J,I),J=1,KDGOF),I=1,KNODE) cc WRITE (*,*) 'NEQ =',NEQ,' F =' cc WRITE(*,7) (F(I),I=1,NEQ) if (itime.le.1) then C.......OPEN LMATRIX FILE OPEN (2,FILE=' ',FORM='UNFORMATTED',STATUS='unknown') CLOSE (2) endif WRITE(*,*) 'SIN_SOLVER MEMORY REQUIRED .... ',MAXA IF (MAXA.GT.MAXT) THEN WRITE(*,*) 'WARNING MATRIX A EXCEED CORE MEMORY .... ',MAXT c STOP 0000 ENDIF CALL REDU(A,U1,JDIAG,NEQ,MAXA,1) C SUBROUTINE REDU(A,B,U,JDIAG,NEQ,MAXA,KKK) 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) endif 3100 CONTINUE 3200 CONTINUE DO 3400 N=1,KNODE WRITE (NOUT,3600) N,(eu(I,N),I=1,KDGOF) 3400 CONTINUE 3600 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,IT,KCOOR,KELEM,K,KK,NNODE,NNE, *ITYP,NCOOR,NUM,TIME,DT,neq,maxa,NODVAR,COOR,NODE,EMATE, &A,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),JDIAG(neq), *es(k,k),em(k),ef(k),eu(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) 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+1 120 R(NR) = COOR(I,JNOD) 130 CONTINUE IMATE = NODE(NNODE*NE) DO 140 J=1,NMATE 140 PRMT(J) = EMATE((IMATE-1)*NMATE+J) PRMT(NMATE+1)=TIME PRMT(NMATE+2)=DT PRMT(NMATE+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 (1,2), ityp 1 call aeq4g2(r,coef,prmt,es,em,ec,ef,ne) goto 3 2 call agl2g2(r,coef,prmt,es,em,ec,ef,ne) goto 3 3 continue C WRITE(*,*) 'ES EM EF =' C DO 555 I=1,K C555 WRITE(*,18) (ES(I,J),J=1,K) C WRITE(*,18) (EM(I),I=1,K) C WRITE(*,18) (EF(I),I=1,K) CC IF (IT.GT.0) THEN do 201 i=1,k do 201 j=1,k Estifn(i,j)=0.0 201 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)