📄 beam0511.for
字号:
PROGRAM beam
!USE IMSL
REAL:: AGL(3,3)
INCLUDE 'COMTM.H'
OPEN(11,STATUS='SCRATCH',FORM='UNFORMATTED')
OPEN(5,FILE='TM0.DAT',STATUS='OLD')
OPEN(6,FILE='TM.OUT',STATUS='UNKNOWN')
OPEN(1,file='temp.log',STATUS='replace')
!AGL(1,1)=2.
!AGL(1,2)=0.5
!AGL(1,3)=-0.09375
!AGL(2,2)=1.8
!AGL(2,3)=-0.01875
!AGL(3,3)=0.4430
!AGL(2,1)=AGL(1,2)
!AGL(3,1)=AGL(1,3)
!AGL(3,2)=AGL(2,3)
!WRITE(*,*)AGL
!BGL=DET(AGL)
!WRITE(*,*)BGL,'=DET'
CALL DATA
CALL INITIAL
DO 30 IINCS=1,NINCS
write(1,*)'run sub inclod iincs=',iincs
CALL INCLOD
! WRITE(1,*)TLOAD,'???'
DO 10 IITER=1,NITER
write(1,*)'run sub nonal niter=',iiter
CALL NONAL
IF(KRESL.EQ.1) CALL STIFFB
write(1,*)'run sub assemb'
CALL ASSEMB
IF(KRESL.EQ.1) CALL GREDUC
IF(KRESL.EQ.2) CALL RESOLV
write(1,*)'run sub baksub'
CALL BAKSUB
write(1,*)'run sub reforb'
CALL REFORB
write(1,*)'run sub conund'
!CALL MONITR
CALL CONUND
IF(NCHEK.EQ.0) GO TO 20
IF(IITER.EQ.1.AND.NOUTP.EQ.1) CALL RESULT
IF(NOUTP.EQ.2) CALL RESULT
10 CONTINUE
WRITE(6,900)
900 FORMAT(1H,5X,'SOLUTION NOT CONVERGED')
STOP
20 CALL RESULT
30 CONTINUE
STOP
END
SUBROUTINE DATA
INCLUDE 'COMTM.H'
DIMENSION ICODE(2),VALUE(2),TITLE(18)
READ(5,965) TITLE
WRITE(6,965) TITLE
965 FORMAT(18A4)
READ(5,*) NPOIN,NELEM,NBOUN,NMATS, NPROP,NNODE,NINCS,NALGO,NDOFN
900 FORMAT(9I5)
WRITE(6,905)NPOIN,NELEM,NBOUN,NMATS,NPROP,NNODE,NINCS,NALGO,NDOFN
905 FORMAT(/1X,'NPOIN=',I5,3X,'NELEM=',I5,3X,'NBOUN='I5,3X,
1 'NMATS=',I5//1X,'NPROP=',I5,3X,'NNODE=',I5,3X,'NINCS=',I5,
1 3X,'NALGO=',I5//1X,'NDOFN=',I5)
NEVAB=NDOFN*NNODE
NSVAB=NDOFN*NPOIN
WRITE(6,910)
910 FORMAT(5X,'MATERIAL PROPERTIES')
DO IMATS=1,NMATS
READ(5,*) JMATS,(PROPS(JMATS,IPROP),IPROP=1,NPROP)
WRITE(6,915) JMATS,(PROPS(JMATS,IPROP),IPROP=1,NPROP)
enddo !IMATS
915 FORMAT(I4,4F20.5)
WRITE(6,920)
920 FORMAT(3X,'EL NODFS MAT.')
DO IELEM=1,NELEM
READ(5,*)JELEM,(LNODS(JELEM,INODE),INODE=1,NNODE),
1 MATNO(JELEM)
WRITE(6,925)JELEM,(LNODS(JELEM,INODE),INODE=1,NNODE),
1 MATNO(JELEM)
ENDDO !IELEM
925 FORMAT(4I5)
WRITE(6,930)
930 FORMAT(1H,5X,'NODE',5X,'COORD.')
DO 30 IPOIN=1,NPOIN
READ(5,*) JPOIN,COORD(JPOIN)
30 WRITE(6,935)JPOIN,COORD(JPOIN)
!ENDDO !IPOIN
935 FORMAT(I10,F15.5)
DO ISVAB=1,NSVAB
IFPRE(ISVAB)=0
PEFIX(ISVAB)=0
ENDDO !ISVAB
IF(NDOFN.EQ.1) WRITE(6,940)
940 FORMAT(1X,'RES.NODE',2X,'CODE',3X,'PRES.VALUES')
IF(NDOFN.EQ.2) WRITE(6,945)
945 FORMAT(1H,1X,'RES.NODE',2X,'CODE',3X,'PRES.VALUES',2X,
1 'CODE',3X,'PRES.VALUES')
DO IBOUN=1,NBOUN
READ(5,*) NODFX,(ICODE(IDOFN),VALUE(IDOFN),IDOFN=1,NDOFN)
WRITE(6,950) NODFX,(ICODE(IDOFN),VALUE(IDOFN),IDOFN=1,NDOFN)
950 FORMAT(I10,2(I5,F15.5))
NPOSN=(NODFX-1)*NDOFN
DO IDOFN=1,NDOFN
NPOSN=NPOSN+1
IFPRE(NPOSN)=ICODE(IDOFN)
PEFIX(NPOSN)=VALUE(IDOFN)
ENDDO !IDOFN
ENDDO !IBOUN
WRITE(6,955)
955 FORMAT(2X,'ELEMENT',10X,'NODAL LOADS')
DO 60 IELEM=1,NELEM
DO 60 IEVAB=1,NEVAB
60 RLOAD(IELEM,IEVAB)=0.0
!ENDDO !IEVAB
!ENDDO !IELEM
70 READ(5,*) JELEM,(RLOAD(JELEM,IEVAB),IEVAB=1,NEVAB)
IF(JELEM.NE.NELEM) GO TO 70
DO IELEM =1,NELEM
WRITE(6,960) IELEM,(RLOAD(IELEM,IEVAB),IEVAB=1,NEVAB)
ENDDO
960 FORMAT(I10,5F15.5)
RETURN
END
SUBROUTINE INITIAL
INCLUDE 'COMTM.H'
DO 20 IELEM=1,NELEM
PLAST(IELEM)=0.0
DO 10 IDOFN=1,NDOFN
10 STRES(IELEM,IDOFN)=0.0
DO 20 IEVAB=1,NEVAB
ELOAD(IELEM,IEVAB)=0.0
20 TLOAD (IELEM,IEVAB)=0.0
DO 30 IPOIN=1,NPOIN
DO 30 IDOFN=1,NDOFN
TDISP(IPOIN,IDOFN)=0.0
30 TREAC(IPOIN,IDOFN)=0.0
RETURN
END
SUBROUTINE INCLOD
INCLUDE 'COMTM.H'
READ (5,*) NITER,NOUTP,FACTO,TOLER
!900 FORMAT(2I5,2F15.5)
WRITE(6,905)IINCS,NITER,NOUTP,FACTO,TOLER
905 FORMAT(5X,'IINCS=',I5,3X,'NITER=',I5,3X,'NOUTP=',I5,
1 3X,'FACTO=',E14.6,3X,'TOLER=',E14.6)
DO IELEM=1,NELEM
DO IEVAB=1,NEVAB
ELOAD(IELEM,IEVAB)=ELOAD(IELEM,IEVAB)+
1 RLOAD(IELEM,IEVAB)*FACTO
TLOAD(IELEM,IEVAB)=TLOAD(IELEM,IEVAB)+
1 RLOAD(IELEM,IEVAB)*FACTO
!10 CONTINUE
ENDDO !IEVAB
ENDDO !IELEM
RETURN
END
SUBROUTINE NONAL
INCLUDE 'COMTM.H'
KRESL=2
IF(NALGO.EQ.1) KRESL=1
IF(NALGO.EQ.2) KRESL=1
IF(NALGO.EQ.3.AND.IINCS.EQ.1.AND.IITER.EQ.1) KRESL=1
IF(NALGO.EQ.4.AND.IITER.EQ.1) KRESL=1
IF(NALGO.EQ.5.AND.IINCS.EQ.1.AND.IITER.EQ.1) KRESL=1
IF(NALGO.EQ.5.AND.IITER.EQ.2) KRESL=1
IF(IITER.EQ.1.OR.NALGO.EQ.1) GO TO 20
DO 10 ISVAB=1,NSVAB
10 FIXED(ISVAB)=0.0
RETURN
20 DO 30 ISVAB=1,NSVAB
30 FIXED(ISVAB)=PEFIX(ISVAB)*FACTO
RETURN
END
SUBROUTINE STIFFB
INCLUDE 'COMTM.H'
REWIND 11
DO 20 IELEM=1,NELEM
LPROP=MATNO(IELEM)
EIVAL=PROPS(LPROP,1)
SVALU=PROPS(LPROP,2)
HARDS=PROPS(LPROP,4)
NODE1=LNODS(IELEM,1)
NODE2=LNODS(IELEM,2)
ELENG=ABS(COORD(NODE2)-COORD(NODE1))
IF(PLAST(IELEM).EQ.0.0)
1 EIVAL=EIVAL*(1.0-EIVAL/(EIVAL+HARDS))
VALU1=0.5*SVALU
VALU2=SVALU/ELENG
VALU3=EIVAL/ELENG
VALU4=0.25*SVALU*ELENG
! VALU4=3*EIVAL/ELENG
ESTIF(1,1)=VALU2
ESTIF(1,2)=VALU1
ESTIF(1,3)=-VALU2
ESTIF(1,4)=VALU1
ESTIF(2,2)=VALU3+VALU4
ESTIF(2,3)=-VALU1
ESTIF(2,4)=-VALU3+VALU4
ESTIF(3,3)=VALU2
ESTIF(3,4)=-VALU1
ESTIF(4,4)=VALU3+VALU4
DO 10 ISTIF=1,4
DO 10 JSTIF=ISTIF,4
10 ESTIF(JSTIF,ISTIF)=ESTIF(ISTIF,JSTIF)
WRITE(11) ESTIF
20 CONTINUE
RETURN
END
SUBROUTINE ASSEMB
INCLUDE 'COMTM.H'
REWIND 11
DO 10 ISVAB=1,NSVAB
10 ASLOD(ISVAB)=0.0
IF(KRESL.EQ.2) GO TO 30
DO 20 ISVAB=1,NSVAB
DO 20 JSVAB=1,NSVAB
20 ASTIF(ISVAB,JSVAB)=0.0
30 CONTINUE
DO 50 IELEM=1,NELEM
READ(11) ESTIF
DO 40 INODE=1,NNODE
NODEI=LNODS(IELEM,INODE)
DO 40 IDOFN=1,NDOFN
NROWS=(NODEI-1)*NDOFN+IDOFN
NROWE=(INODE-1)*NDOFN+IDOFN
ASLOD(NROWS)=ASLOD(NROWS)+ELOAD(IELEM,NROWE)
IF(KRESL.EQ.2) GO TO 40
DO JNODE=1,NNODE
NODEJ=LNODS(IELEM,JNODE)
DO JNOFN=1,NDOFN
NCOLS=(NODEJ-1)*NDOFN+JNOFN
NCOLE=(JNODE-1)*NDOFN+JNOFN
ASTIF(NROWS,NCOLS)=ASTIF(NROWS,NCOLS)+ESTIF(NROWE,NCOLE)
ENDDO
ENDDO
40 CONTINUE
50 CONTINUE
!WRITE(1,*)aslod
!WRITE(1,*)'????????????'
RETURN
END
SUBROUTINE GREDUC
INCLUDE 'COMTM.H'
KOUNT=0
NEQNS=NSVAB
!WRITE(1,*)ASTIF
!WRITE(1,*)'!!!!!!!!!!!!!!!!!!!'
!WRITE(1,*)(ASTIF(II,II),II=1,NEQNS)
DO 70 IEQNS=1,NEQNS
IF(IFPRE(IEQNS).EQ.1) GO TO 40
PIVOT=ASTIF(IEQNS,IEQNS) !+1.E-5
!WRITE(1,*)'?1HERE IS RUNNING PIVOT',PIVOT,IEQNS
IF(ABS(PIVOT).LT.1.0E-10) GO TO 60
IF(IEQNS.EQ.NEQNS) GO TO 70
IEQN1=IEQNS+1
DO 30 IROWS=IEQN1,NEQNS
KOUNT=KOUNT+1
FACTR=ASTIF(IROWS,IEQNS)/PIVOT
FRESV(KOUNT)=FACTR
IF(FACTR.EQ.0.0) GO TO 30
DO ICOLS=IEQNS,NEQNS
ASTIF(IROWS,ICOLS)=ASTIF(IROWS,ICOLS)-FACTR*ASTIF(IEQNS,ICOLS)
ENDDO
! WRITE(1,*)'1HERE IS RUNNING ASTIF(','=',IROWS,',',IROWS,')',
! &ASTIF(IROWS,IROWS)
!10 CONTINUE
ASLOD(IROWS)=ASLOD(IROWS)-FACTR*ASLOD(IEQNS)
!ENDDO
30 CONTINUE
GO TO 70
40 DO 50 IROWS=IEQNS,NEQNS
ASLOD(IROWS)=ASLOD(IROWS)-ASTIF(IROWS,IEQNS)*FIXED(IEQNS)
50 CONTINUE
GO TO 70
60 WRITE(6,900)
900 FORMAT(5X,15HINCORRECT PIVOT)
STOP
70 CONTINUE
RETURN
END
SUBROUTINE BAKSUB
INCLUDE 'COMTM.H'
DO IPOIN=1,NPOIN
DO IDOFN=1,NDOFN
TDISP(IPOIN,IDOFN)=0.
TREAC(IPOIN,IDOFN)=0.
ENDDO
ENDDO
NEQNS=NSVAB
DO 10 IEQNS=1,NEQNS
REACT(IEQNS)=0.0
10 CONTINUE
NEQN1=NEQNS+1
DO 40 IEQNS=1,NEQNS
NBACK=NEQN1-IEQNS
PIVOT=ASTIF(NBACK,NBACK)
RESID=ASLOD(NBACK)
IF(NBACK.EQ.NEQNS) GO TO 30
NBAC1=NBACK+1
DO 20 ICOLS=NBAC1,NEQNS
!WRITE(1,*)'1HERE IS RUNNING PIVOT=',XDISP(ICOLS)
RESID=RESID-ASTIF(NBACK,ICOLS)*XDISP(ICOLS)
20 CONTINUE
30 IF(IFPRE(NBACK).EQ.0) XDISP(NBACK)=RESID/PIVOT
IF(IFPRE(NBACK).EQ.1) XDISP(NBACK)=FIXED(NBACK)
IF(IFPRE(NBACK).EQ.1) REACT(NBACK)=-RESID
40 CONTINUE
KOUNT=0
! TDISP(1:NPOIN,1:NDOFN)=0.
! TREAC(1:NPOIN,1:NDOFN)=0.
DO 50 IPOIN=1,NPOIN
DO 50 IDOFN=1,NDOFN
KOUNT=KOUNT+1
!!!!!!!!!!!!!!!!!!REVISED 07-04-22!!!!!!!!!!!!!!
TDISP(IPOIN,IDOFN)=TDISP(IPOIN,IDOFN)+XDISP(KOUNT)
50 TREAC(IPOIN,IDOFN)=TREAC(IPOIN,IDOFN)+REACT(KOUNT)
! WRITE(1,*)'HERE IS RUNNING TLOAD=',TDISP
! TDISP(IPOIN,IDOFN)=XDISP(KOUNT)
!50 TREAC(IPOIN,IDOFN)=REACT(KOUNT)
!!!!!!!!!!!!ADD THE SECTION!!!!!!!!!!07-04-22!!!!!!!
DO 90 IPOIN=1,NPOIN
DO 60 IELEM=1,NELEM
DO 60 INODE=1,NNODE
NLOCA=LNODS(IELEM,INODE)
60 IF (IPOIN.EQ.NLOCA) GO TO 70
70 DO 80 IDOFN=1,NDOFN
NPOSN=(IPOIN-1)*NDOFN+IDOFN
IEVAB=(INODE-1)*NDOFN+IDOFN
!write(1,*)'2 here is running, IEVAB=',NPOSN,IDOFN
80 TLOAD(IELEM,IEVAB)=TLOAD(IELEM,IEVAB)+REACT(NPOSN)
90 CONTINUE
!!!!!!!!!!!!!!!!!ADD THE ABOVE SECTION!!!07-04-22!!!!
RETURN
END
SUBROUTINE RESOLV
INCLUDE 'COMTM.H'
KOUNT=0
NEQNS=NSVAB
DO 40 IEQNS=1,NEQNS
IF(IFPRE(IEQNS).EQ.1) GO TO 20
C
C* REDUCE RHS
C
IF(IEQNS.EQ.NEQNS) GO TO 40
IEQN1=IEQNS+1
DO 10 IROWS=IEQN1,NEQNS
KOUNT=KOUNT+1
FACTR=FRESV(KOUNT)
IF(FACTR.EQ.0) GO TO 10
ASLOD(IROWS)=ASLOD(IROWS)-FACTR*ASLOD(IEQNS)
10 CONTINUE
GO TO 40
20 DO 30 IROWS=IEQNS,NEQNS
ASLOD(IROWS)=ASLOD(IROWS)-ASTIF(IROWS,IEQNS)*FIXED(IEQNS)
30 CONTINUE
40 CONTINUE
RETURN
END
SUBROUTINE REFORB
INCLUDE 'COMTM.H'
DO 10 IELEM=1,NELEM
DO 10 IEVAB=1,NEVAB
10 ELOAD(IELEM,IEVAB)=0.0
DO 70 IELEM=1,NELEM
LPROP=MATNO(IELEM)
EIVAL=PROPS(LPROP,1)
SVALU=PROPS(LPROP,2)
YIELD=PROPS(LPROP,3)
HARDS=PROPS(LPROP,4)
NODE1=LNODS(IELEM,1)
NODE2=LNODS(IELEM,2)
ELENG=ABS(COORD(NODE2)-COORD(NODE1))
WNOD1=XDISP(NODE1*NDOFN-1)
WNOD2=XDISP(NODE2*NDOFN-1)
THTA1=XDISP(NODE1*NDOFN)
THTA2=XDISP(NODE2*NDOFN)
STRAN=(THTA1-THTA2)/ELENG
STLIN=STRAN*EIVAL
STCUR=STRES(IELEM,1)+STLIN
PREYS=YIELD+HARDS*ABS(PLAST(IELEM))
IF(ABS(STRES(IELEM,1)).GE.PREYS) GO TO 20
ESCUR=ABS(STCUR)-PREYS
IF(ESCUR.LE.0.0) GO TO 40
RFACT=ESCUR/ABS(STLIN)
GO TO 30
20 IF(STRES(IELEM,1).GT.0.0.AND.STLIN.LE.0.0) GO TO 40
IF(STRES(IELEM,1).LT.0.0.AND.STLIN.GT.0.0) GO TO 40
RFACT=1.0
30 REDUC=1.0-RFACT
STRES(IELEM,1)=STRES(IELEM,1)+REDUC*STLIN+RFACT*EIVAL*
1 (1.0-EIVAL/(EIVAL+HARDS))*STRAN
PLAST(IELEM)=PLAST(IELEM)+RFACT*STRAN*EIVAL/(EIVAL+HARDS)
GO TO 50
40 STRES(IELEM,1)=STRES(IELEM,1)+STLIN
50 STRES(IELEM,2)=STRES(IELEM,2)+(SVALU/ELENG)*(WNOD2-WNOD1)
1 -0.5*SVALU*(THTA1+THTA2)
ELOAD(IELEM,1)=ELOAD(IELEM,1)-STRES(IELEM,2)
ELOAD(IELEM,2)=ELOAD(IELEM,2)+STRES(IELEM,1)
1 -0.5*ELENG*STRES(IELEM,2)
ELOAD(IELEM,3)=ELOAD(IELEM,3)+STRES(IELEM,2)
ELOAD(IELEM,4)=ELOAD(IELEM,4)-STRES(IELEM,1)
1 -0.5*ELENG*STRES(IELEM,2)
70 CONTINUE
RETURN
END
SUBROUTINE CONUND
INCLUDE 'COMTM.H'
DIMENSION STFOR(52),TOFOR(52)
NCHEK=0
RESID=0.0
RETOT=0.0
DO 10 ISVAB=1,NSVAB
STFOR(ISVAB)=0.0
10 TOFOR(ISVAB)=0.0
DO 20 IELEM=1,NELEM
IEVAB=0
DO 20 INODE=1,NNODE
NODNO=LNODS(IELEM,INODE)
DO 20 IDOFN=1,NDOFN
IEVAB=IEVAB+1
NPOSN=(NODNO-1)*NDOFN+IDOFN
STFOR(NPOSN)=STFOR(NPOSN)+ELOAD(IELEM,IEVAB)
20 TOFOR(NPOSN)=TOFOR(NPOSN)+TLOAD(IELEM,IEVAB)
DO 30 ISVAB=1,NSVAB
REFOR=TOFOR(ISVAB)-STFOR(ISVAB)
RESID=RESID+REFOR*REFOR
30 RETOT=RETOT+TOFOR(ISVAB)*TOFOR(ISVAB)
DO 40 IELEM=1,NELEM
DO 40 IEVAB=1,NEVAB
40 ELOAD(IELEM,IEVAB)=TLOAD(IELEM,IEVAB)-ELOAD(IELEM,IEVAB)
WRITE(*,*)RESID,RETOT,'RESID,RETOT'
RATIO=100.0*SQRT(RESID/RETOT)
IF(RATIO.GT.TOLER) NCHEK=1
IF(IITER.EQ.1) GO TO 50
IF(RATIO.GT.PVALU) NCHEK=999
50 PVALU=RATIO
WRITE(6,900) IITER,NCHEK,RATIO
900 FORMAT(5X,'ITERATION NMBER=',I5,/
& 5X,'CONVERGENCE CODE=',
1 I4,3X,'NORM OF RESIDUAL SUM RATIO=',E14.6)
RETURN
END
SUBROUTINE MONITR
INCLUDE 'COMTM.H'
NCHEK=0
RCURR=0.0
DO 10 IPOIN=1,NPOIN
10 RCURR=RCURR+TDISP(IPOIN,1)*TDISP(IPOIN,1)
IF(IITER.EQ.1) RINTL=RCURR
IF(IITER.EQ.1) NCHEK=1
IF(IITER.EQ.1) GO TO 20
RATIO=100.0*SQRT(ABS(RCURR-PVALU))/SQRT(RINTL)
IF(RATIO.GT.TOLER) NCHEK=1
20 PVALU=RCURR
WRITE(6,900) NCHEK,RATIO
900 FORMAT(5X,'CNVERGENCE CODE=',I4,3X,'NORM OF
1 RESIDUAL SUM RATIO=',E14.6)
RETURN
END
SUBROUTINE RESULT
INCLUDE 'COMTM.H'
IF(NDOFN.EQ.1) WRITE(6,900)
900 FORMAT(5X,'NODE',4X,'DISPL.',12X,'REACTIONS')
IF(NDOFN.EQ.2) WRITE(6,910)
910 FORMAT(5X,'NODE',4X,'DISP.',12X,'REACTION',7X,'DISPL.',12X,
1 'REACTION')
DO 10 IPOIN=1,NPOIN
10 WRITE(6,920) IPOIN,(TDISP(IPOIN,IDOFN),TREAC(IPOIN,IDOFN),
1 IDOFN=1,NDOFN)
920 FORMAT(I10 ,2(E14.6,5X,E14.6))
IF(NDOFN.EQ.2) WRITE(6,930)
930 FORMAT(2X,'ELEMENT',12X,'STRESSES',12X,'PL.STRAIN')
IF(NDOFN.EQ.1) WRITE(6,940)
940 FORMAT(2X,'ELEMENT',5X,'STRESSES',5X,'PL.STRAIN')
DO 20 IELEM=1,NELEM
20 WRITE(6,950) IELEM,(STRES(IELEM,IDOFN),IDOFN=1,NDOFN),
1 PLAST(IELEM)
950 FORMAT(I10,3E14.6)
RETURN
END
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -