结构力学电算程序

DIMENSION JE(2,100),JN(3,100),JC(6),EA(100),EI(100),X(100),Y(100),PJ(2,50),PF(4,100)
REAL*8 KE(6,6),KD(6,6),T(6,6),P(300),KB(200,20),F(6),F0(6),D(6),BL,SI,CO,S,C
OPEN(5,FILE='FRAM2.DAT')


READ(5,*) NE,NJ,N,NW,NPJ,NPF
READ(5,*) (X(J),Y(J),(JN(I,J),I=1,3),J=1,NJ)
READ(5,*) ((JE(I,J),I=1,2),EA(J),EI(J),J=1,NE)
IF(NPJ.NE.0) READ(5,*) ((PJ(I,J),I=1,2),J=1,NPJ)
IF(NPF.NE.0) READ(5,*) ((PF(I,J),I=1,4),J=1,NPF)
WRITE(*,10) NE,NJ,N,NW,NPJ,NPF
WRITE(*,20) (J,X(J),Y(J),(JN(I,J),I=1,3),J=1,NJ)
WRITE(*,30) (J,(JE(I,J),I=1,2),EA(J),EI(J),J=1,NE)
IF(NPJ.NE.0) WRITE(*,40) ((PJ(I,J),I=1,2),J=1,NPJ)
IF(NPF.NE.0) WRITE(*,50) ((PF(I,J),I=1,4),J=1,NPF)
10 FORMAT(/6X,'NE=',I5,2X,'NJ=',I5,2X,'N=',I5,2X,'NW=',I5,2X,'NPJ=',I5,2X,'NPF=',I5)
20 FORMAT(/7X,'NODE',7X,'X',11X,'Y',12X,'XX',8X,'YY',8X,'ZZ'/(1X,I10,2F12.4,3I10))
30 FORMAT(/4X,'ELEMENT',4X,'NODE-I',4X,'NODE-J',11X,'EA',13X,'EI'/(1X,3I10,2E15.6))
40 FORMAT(/7X,'CODE',7X,'PX-PY-PM'/(1X,F10.0,F15.4))
50 FORMAT(/4X,'ELEMENT',7X,'IND',10X,'A',14X,'Q',/(1X,2F10.0,2E15.4))
DO 55 I=1,N
55 P(I)=0.D0
IF(NPJ.EQ.0) GO TO 65
DO 60 I=1,NPJ
L=PJ(1,I)
60 P(L)=PJ(2,I)
65 IF(NPF.EQ.0) GO TO 90
DO 70 I=1,NPF
M=PF(1,I)
CALL SCL(M,NE,NJ,BL,SI,CO,JE,X,Y)
CALL EFX(I,NPF,BL,PF,F0)
CALL CTM(SI,CO,T)
CALL EJC(M,NE,NJ,JE,JN,JC)
DO 75 L=1,6
S=0.D0
DO 80 K=1,6
80 S=S-T(K,L)*F0(K)
F(L)=S
75 CONTINUE
DO 85 J=1,6
L=JC(J)
IF(L.EQ.0) GO TO 85
P(L)=P(L)+F(J)
85 CONTINUE
70 CONTINUE
90 DO 95 I=1,N
DO 100 J=1,NW
100 KB(I,J)=0.D0
95 CONTINUE
DO 105 M=1,NE
CALL SCL(M,NE,NJ,BL,SI,CO,JE,X,Y)
CALL CTM(SI,CO,T)
CALL ESM(M,NE,BL,EA,EI,KD)
CALL EJC(M,NE,NJ,JE,JN,JC)
DO 110 I=1,6
DO 115 J=1,6
S=0.D0
DO 120 L=1,6
DO 125 K=1,6
125 S=S+T(L,I)*KD(L,K)*T(K,J)
120 CONTINUE
KE(I,J)=S
115 CONTINUE
110 CONTINUE
DO 130 L=1,6
I=JC(L)
IF(I.EQ.0) GO TO 130
DO 135 K=1,6
J=JC(K)
IF(J.EQ.0.OR.J.LT.I) GO TO 135
JJ=J-I+1
KB(I,JJ)=KB(I,JJ)+KE(L,K)
135 CONTINUE
130 CONTINUE
105 CONTINUE
N1=N-1
DO 140 K=1,N1
IM=K+NW-1
IF(N.LT.IM) IM=N
I1=K+1
DO 145 I=I1,IM
L=I-K+1
C=KB(K,L)/KB(K,1)
JM=NW-L+1
DO 150 J=1,JM
JJ=J+I-K
150 KB(I,J)=KB(I,J)-C*KB(K,JJ)
145 P(I)=P(I)-C*P(K)
140 CONTINUE
P(N)=P(N)/KB(N,1)
DO 155 K=1,N1
I=N-K
JM=K+1
IF(NW.LT.JM) JM=NW
DO 160 J=2,JM
L=J+I-1
160 P(I)=P(I)-KB(I,J)*P(L)
155 P(I)=P(I)/KB(I,1)
WRITE(*,165)
165 FORMAT(/7X,'NODE',10X,'U',14X,'V',11X,'CATA')
DO 170 I=1,NJ
DO 175 J=1,3
D(J)=0.D0
L=JN(J,I)
IF(L.EQ.0) GO TO 175
D(J)=P(L)
175 CONTINUE
WRITE(*,180) I,D(1),D(2),D(3)
180 FORMAT(1X,I10,3E15.6)
170 CONTINUE
WRITE(*,200)
200 FORMAT(/4X,'ELEMENT',13X,'N',17X,'Q',17X,'M')
DO 205 M=1,NE
CALL SCL(M,NE,NJ,BL,SI,CO,JE,X,Y)
CALL ESM(M,NE,BL,EA,EI,KD)
CALL CTM(SI,CO,T)
CALL EJC(M,NE,NJ,JE,JN,JC)
DO 210 I=1,6
L=JC(I)
D(I)=0.D0
IF(L.EQ.0) GO TO 210
D(I)=P(L)
210 CONTINUE
DO 220 I=1,6
F(I)=0.D0
DO 230 J=1,6
DO 240 K=1,6
240 F

(I)=F(I)+KD(I,J)*T(J,K)*D(K)
230 CONTINUE
220 CONTINUE
IF(NPF.EQ.0) GO TO 270
DO 250 I=1,NPF
L=PF(1,I)
IF(M.NE.L) GO TO 250
CALL EFX(I,NPF,BL,PF,F0)
DO 260 J=1,6
260 F(J)=F(J)+F0(J)
250 CONTINUE
270 WRITE(*,280) M,(F(I),I=1,6)
280 FORMAT(/1X,I10,3X,'N1=',F12.4,3X,'Q1=',F12.4,3X,'M1=',F12.4/14X,'N2=',F12.4,3X,'Q2=',F12.4,3X,'M2=',F12.4)
205 CONTINUE
CLOSE(5)
STOP
END
SUBROUTINE EJC(M,NE,NJ,JE,JN,JC)
DIMENSION JE(2,NE),JN(3,NJ),JC(6)
J1=JE(1,M)
J2=JE(2,M)
DO 10 I=1,3
JC(I)=JN(I,J1)
10 JC(I+3)=JN(I,J2)
RETURN
END
SUBROUTINE SCL(M,NE,NJ,BL,SI,CO,JE,X,Y)
DIMENSION JE(2,NE),X(NJ),Y(NJ)
REAL*8 BL,SI,CO,DX,DY
J1=JE(1,M)
J2=JE(2,M)
DX=X(J2)-X(J1)
DY=Y(J2)-Y(J1)
BL=DSQRT(DX*DX+DY*DY)
SI=DY/BL
CO=DX/BL
RETURN
END
SUBROUTINE ESM(M,NE,BL,EA,EI,KD)
DIMENSION EA(NE),EI(NE)
REAL*8 KD(6,6),BL,S,G,G1,G2,G3
G=EA(M)/BL
G1=2.D0*EI(M)/BL
G2=3.D0*G1/BL
G3=2.D0*G2/BL
DO 10 I=1,6
DO 10 J=1,6
10 KD(I,J)=0.D0
KD(1,1)=G
KD(1,4)=-G
KD(4,4)=G
KD(2,2)=G3
KD(5,5)=G3
KD(2,5)=-G3
KD(2,3)=-G2
KD(2,6)=-G2
KD(3,5)=G2
KD(5,6)=G2
KD(3,3)=2.D0*G1
KD(6,6)=2.D0*G1
KD(3,6)=G1
DO 20 I=1,5
I1=I+1
DO 30 J=I1,6
30 KD(J,I)=KD(I,J)
20 CONTINUE
RETURN
END
SUBROUTINE CTM(SI,CO,T)
REAL*8 T(6,6),SI,CO
DO 10 I=1,6
DO 10 J=1,6
10 T(I,J)=0.D0
T(1,1)=CO
T(1,2)=SI
T(2,1)=-SI
T(2,2)=CO
T(3,3)=1.D0
DO 20 I=1,3
DO 20 J=1,3
20 T(I+3,J+3)=T(I,J)
RETURN
END
SUBROUTINE EFX(I,NPF,BL,PF,F0)
DIMENSION PF(4,NPF)
REAL*8 F0(6),A,B,C,G,Q,S,BL
IND=PF(2,I)
A=PF(3,I)
Q=PF(4,I)
C=A/BL
G=C*C
B=BL-A
DO 5 J=1,6
5 F0(J)=0.D0
GO TO (10,20,30,40,50,60,70),IND
10 S=Q*A*0.5D0
F0(2)=-S*(2.D0-2.D0*G+C*G)
F0(5)=-S*G*(2.D0-C)
S=S*A/6.D0
F0(3)=S*(6.D0-8.D0*C+3.D0*G)
F0(6)=-S*C*(4.D0-3.D0*C)
GO TO 100
20 S=B/BL
F0(2)=-Q*S*S*(1.D0+2.D0*C)
F0(5)=-Q*G*(1.D0+2.D0*S)
F0(3)=Q*S*S*A
F0(6)=-Q*B*G
GOTO 100
30 S=B/BL
F0(2)=-6.D0*Q*C*S/BL
F0(5)=-F0(2)
F0(3)=Q*S*(2.D0-3.D0*S)
F0(6)=Q*C*(2.D0-3.D0*C)
GO TO 100
40 S=Q*A*0.25D0
F0(2)=-S*(2.D0-3.D0*G+1.6D0*G*C)
F0(5)=-S*G*(3.D0-1.6D0*C)
S=S*A
F0(3)=S*(2.D0-3.D0*C+1.2D0*G)/1.5D0
F0(6)=-S*C*(1.D0-0.8D0*C)
GO TO 100
50 F0(1)=-Q*A*(1.D0-0.5D0*C)
F0(4)=-0.5D0*Q*C*A
GOTO 100
60 F0(1)=-Q*B/BL
F0(4)=-Q*C
GO TO 100
70 S=B/BL
F0(2)=-Q*G*(3.D0*S+C)
F0(5)=-F0(2)
S=S*B/BL
F0(3)=-Q*S*A
F0(6)=Q*G*B
100 RETURN
END