yuandaima.rar.rar_4 3 2 1_txt_有限元 fortran_混凝土

C FEM.FOR C 平面四节点等参单元程序 IMPLICIT REAL*4(A-H,O-Z),INTEGER*2(I-N) CHARACTER CH (80)*1,DFALE*12,OFALE*12,YN*1 INTEGER*2 LI(20000) REAL*4 A(10000) EQUIVALENCE (A,LI) DFALE='FEM.DAT' OFALE='FEM.OUT' CALL OPE('平面四节点等参单元程序','(FEM 1.0)',36,18, $DFALE,OFALE,YN) IF (YN.EQ.'N'.OR.YN.EQ.'n') THEN ITYPE=0 IMESH=1 ISG=1 INW=1 WRITE(*,'(2X,4(A,12),A\)') $ '问题类型，平面形状，应力输出，刚度输出[ $ ',ITYPE,',',IMESH,',',ISG,',',INW,']:' CALL ICDA(CH,L1,80,80,NIS,ITYPE,IMESH,ISG,INW,K5,K6) NX=3 XY=2 NFP=3 ND=0 NTP=0 NPP=2 WRITE(*,'(2X,2A/2X,6(A,I3),A\)') $ 'X方向分点数，Y方向分点数，支撑约束数,', $ '附力结点数，附加单元数，节点荷载数','[' $ NX,',',NY,',',NFP,',',ND,',',NTP,',',NPP,']:' CALL ICDA(CH,LI,80,80,NIS,NX,NY,NEP,ND,NTP,NPP) OPEN(7,FILE=DFALE,STATUS='NEW') WRITE(7,*)ITYPE,IMESH,ISG,INW,NX,NY,NFP,ND,NTP,NPP ELSE OPEN(7,FILE=DFALE,STATUS='OLD') READ(7,*)ITYPE,IMESH,ISG,INW,NX,NY,NFP,ND,NTP,NPP END IF WRITE(11,'(/4X,//4X,A)')'输入数据:' WRITE(11,'(/4X,4A12/4I12)') $'问题类型','平面形状','应力输出','刚度输出', $ITYPE,IMESH,ISG,INW WRITE(11,'(/4X,6A12/6I12)') $'X方向分点数','Y方向分点数','支撑约束数', $'附力节点数','附加单元数','节点荷载数',NX,NY,NFP,ND,NTP,NPP C NNOP=NC*NY+ND NNEL=(NX-1)*(NY-1)+NTP NUVD=2*NNOP ixt=1 iyt=ixt+nnop ibkkt=iyt+nnop irnod＝ibkkt+nuvd*nuvd ixx1=irnod+2*nd iutixx1+nuvd ivt=iut+nnop ixa=iroad+npp iya＝ixa+4 mal=iya+4 jntoe=2*mal+2 jnod=jntoe+4*nnel jntoes=jnod+nd jload=jntoes+5*ntp jnsp=jload+2*npp nal=jnsp+2*nfp mal=nal/2+2 IF(YN.EQ.'N'.OR.YN.EQ.'n')THEN IF(IMESH.EQ.1)THEN A(IXA)=20.0 A(IXA+1)=15.0 A(IXA+2)=0.0 A(IXA+3)=0.5 WRITE(*,'(2X,4(A,F6.2).A\)') $ '高,宽,始点坐标X,始点坐标Y $ [',(A(IXA+J),',',J=0.2),A(IXA+3),']:' CALL RCDA(CH,A(IXA+J),J=0,3) ELSE DO 54 J=IXA,IYA+3 A(J)=0.0 DO 55 J=0,1 A(IXA+J+J+1)=10.0 A(IYA+J+2)=10. 55 CONTINUE WRITE(*,'(2X,2A/2X,A,4(F6.2,A,F6.2).A\)') $ '输入四个角点坐标(左下,右下,左上,右上)', $ 'X1,Y1 X2,Y2 X3,Y3 X4,Y4', $ '[',(A(IXA+J),',',A(IXA+J),J=0,3),']:' CALL RCDA(CH,A(MAL+1),80,80,NIS,X1,X2,X3,X4,X5,X6) DO 60 J=1,4 J2=J+J J1=J2-1 IF(J1.GT.NIS)GO TO 70 A(IXA+J-1)=A(MAL+J1) IF (J2.GT.NIS)GO TO 70 A(IYA+J-1)=A(MAL+J2) 60 CONTINUE 70 WRITE(7,*)(A(IXA+J),A(IYA+J),J=0,3) END IF ELSE IF(IMESH.EQ.1)THEN READ(7,*)(A(IXA+J),J=0,3) ELSE READ(7,*)(A(IXA+J),A(IYA+J),J=0,3) END IF END IF IF(IMESH.EQ.1)THEN WRITE(11,'(/2X,4A12/4F12.3)') $ '高','宽','始点坐标X','始点坐标Y',(A(IXA+J),J=0,3) ELSE WRITE(11,'(4X,A/1X,8A8/2X,8F8.2)') $ '四个角点坐标(左下,右下,左上,右上)', $ 'X1,Y1 X2,Y2 X3,Y3 X4,Y4',(A(IXA+J),A(IYA+J),J=0,3) END IF C IF(ND.GT.0)THEN IF(YN.EQ.'N'.OR.YN.EQ.'n')THEN LI(JNOD)=NNOP+1-ND A(IRNOD)=1.0 A(IRNOD+1)=1.0 LH=0 20 J1=LH J2=J1+ND IF(LH.GT.0) THEN LI(JNOD+LH)=LI(GNOD+LH-1)+1 A(IRNOD+J1)=A(IRNOD+J1-1) A(IRNOD+J2)=A(IRNOD+J2-1) END IF WRITE(*,'(2X,A,I2,A,I3,A,F8.3,A\)') $ '第',LH+1,'个附加点：结点号,X方向坐标,Y方向坐标[,', $ LI(JNOD+LH),',',A(IRNOD+J1),',',=A(IRNOD+J2),']:' CALL RCDA(CH,A(IRNOD+J1),80,80,NIS,X1,X2,X3,X4,X5,X6) JI=0 DO 23 J=1,NIS JI=(J-1)/3 JL=J-3*JI X1=A(MAL+J) IF(JK.EQ.1)LI(JNOD+LH+JI)=X1 IF(JK.EQ.2)A(IRNOD+J1+JI)=X1 IF(JK.EQ.3)A(IRNOD+J2+JI)=X1 23 CONTINUE LH=LH+JI+1 IF(LH.LT.ND)GO TO 20 WRITE(7,*)(LI(JNOD+J),A(IRNOD+J),A(IRNOD+J+ND),J=0,ND-1) ELSE READ(7,*)(LI(JNOD+J),A(IRNOD+J),A(IRNOD+J+ND),J=0,ND-1) END IF END IF IF(NTP.GT.0)THEN IF(YN.EQ.'N'.OR.YN.EQ.'n')THEN J3=4 L1(JNTOES)=NNEL+1-NTP DO 28 J=1,2 28 LI(JNTOES+J*NTP)=NNOP-ND-NTP-1+J LI(JNTOES+J*NTP+2*NTP)=NNOP-NTP+2-J LH=0 30 IF(LH.GT.0)THEN DO 32 J=0,J3*NTP,NTP 32 LI(JNTOES+LH+J)=LI(JNTOES+LH+J-1)+1 ENDIF WRITE(*,'(2X,A,I2,A\)') $ '第',LH+1,'个附加单元:单元号,点1,点2,点3', WRITE(*,'(A\)')'点4' WRITE(*,'(4(A,I3)\)')'[', $ (LI(JNTOES+LH+J*NTP),',',J=0,2),LI(JNTOES+LH+3*NTP) WRITE(*,'(A,I3\)')',',LI(JNTOES+LH+4*NTP) WRITE(*,'(A\)')']:' CALL ICDA(CH,LI(NAL+1),80,80,NIS,I1,I2,I3,I4,I5,I6) JI=0 DO 34 J=1,NIS JI=(J-1)/(J3+1) JK=J-JI*(J3+1) DO 33 K=1,J3+1 IF(JK.EQ.K) THEN LI(JNTOES+LH+K*NTP-NTP)=LI(NAL+J) GO TO 34 END IF 33 CONTINUE 34 CONTINUE LH=LH+JI+1 IF(LH.LT.NTP)GO TO 30 WRITE(7,*)(LI(JNTOES+J),J=0,J3*NTP+NTP-1) ELSE J4+5 READ(7,*)(LI(JNTOES+J),J=0,J4*NTP+NTO-1) END IF END IF C IF(ND.GT.0)THEN WRITE(11,'(/4X,A)')'附加节点信息:' WRITE(11,'(3A12)')'结点号,X方向坐标,Y方向坐标' WRITE(11,'(I10,F12.2,F12.2)') $ (LI(JNOD+J),A(IRNOD+J),A(IRNOD+J+ND),J=0,ND-1) END IF IF(NTP.GT.0)THEN WRITE(11,'(/4X,A)')'附加单元信息:' WRITE(11,'(4X,A6,3A10\)')'单元号','点1','点2','点3' WRITE(11,'(A\)')'点4' WRITE(11,'(A)')' WRITE(11,'(I8,I11,3I10)') $ ((LI(JNTOES+J*NTP+I),J=0,4)I=0,NTP-1) END IF C IF(YN.EQ.'N'.OR.YN.EQ.'n')THEN EP=0.0001 EMOD=2.1E6 GNU=0.22 THKU=1.0 WRITE(*,'(2X,A,F7.5,A,F10.0,A,2(F6.2,A)\)') $ '迭代精度,弹性模量,泊桑比,单元厚度[', $ EP,',',EMOD,',',GNU,',',THKU,']:' CALL RCDA(CH,A(MAL+1),80,80,NIS,EP,EMOD,GNU,THKU,X5,X6) WRITE(7,*)EP,EMOD,GNU,THKU ELSE READ(7,*)EP,EMOD,GNU,THKU END IF WRITE(11,'(/2X,6A12)')'迭代精度','弹性模量','泊桑比','单元厚度', WRITE(11,'(F12.4,F14.0,F10.3,F12.2)')EP,EMOD,GNU,THKU C IF(YN.EQ.'N'.OR.YN.EQ.'n')THEN WRITE(*,'(4X,A)')'输入荷载信息:' LH=0 40 LI(JLODA+LH)=LH+1 LI(JLODA+LH+NPP)=2 A(ILODA+LH)=1.0 WRITE(*,'(4X,AI2,A,2(I3,A)\)') $ '第',LH+1,'行:节点号,方向(1--X,2--Y),荷载值[', $ LI(JLOAD+LH),',',LI(JLOAD+LH+NPP),',',A(IROAD+LH),']:' CALL RCDA(CH,A(MAL+1),80,80,NIS,X1,X2,X3,X4,X5,X6) JI=0 DO 45 J=1,NIS JI=(J-1)/3 JL=J-3*JI X1=A(MAL+J) IF(JK.EQ.1)LI(JLOAD+LH+JI)=X1 IF(JK.EQ.2)LI(JLOAD+LH+JI+NPP)=X1 IF(JK.EQ.3)A(IROAD+LH+JI)=X1 45 CONTINUE LH=LH+JI+1 IF(LH.LT.NPP)GO TO 40 C LH=0 WRITE(*,'(4X,8A)')'输入支座约束信息:' 10 LI(JNSP+LH)=LH+1 LI(JNSP+LH+NFP)=2 WIRTE(*,'(2X,3(A,I2),A\)') $ '第',LH+1,'个约束:节点号,方向(1--X,2--Y),荷载值[', $ LI(JNSP+LH),',',LI(JNSP+LH+NFP),']:' CALL ICDA(CH,LI(NAL+1),80,80,NIS,I1,I2,I3,I4,I5,I6) JI=0 DO 25 J=1,NIS JI=(J-1)/2 JK=J-2*JI IF(JK.EQ.1)LI(JNSP+LH+JI)=LI(NAL+J) IF(JK.EQ.2)LI(JNSP+LH+JI+NFP)=LI(NAL+J) 25 CONTINUE LH=LH+JI+1 IF(LH.LT.NTP)GO TO 10 WRITE(7,*)(LI(JLOAD+I),LI(JLOAD+I+NPP),A(AROAD+I),I=0,NPP-1) WRITE(7,*)(LI(JNSP+I),LI(JNSP+I+NFP),I=0,NFP-1) ELSE IF WRITE(11,'(4X,A/4,4A10)')'荷载信息:', $'行号','节点号','方向','荷载值' WRITE(11,'(4X,A,I3,A,I8,/10,F12.3)') $('(',I+1')',(LI(JLOAD+I+J),J=0,NPP,NPP),A(AROAD+I),I=0,NPP-1) WRITE(11,'(/4X,A)')'约束信息(1--X方向约束,2--Y方向约束)' WRITE(11,'(4X,2A12/(2I12))')'约束节点号,约束方向', $(LI(JNSP+I),LI(JNSP+I+NFP),I=0,NFP-1) WRITE(11,'(/4X,A)')'结算结果部分:' CALL FEM(A(IXT),A(IYT),A(IBKKT),A(IRNOD),A(IXX1), $A(IUT),A(IVT),A(IROAD),A(IXA),A(IYA),EP,EMOD,GNU,THKU, $LI(JNTOE),LI(JNOD),LI(JNTOES),LI(JLOAD),LI(JNSP), $ITYPE,IMESH,ISG,INW,NX,NY,NFP,ND,NTP,NPP,NNOP,NNEL,