DIMENSION JE(2,100),JN(3,100),JC(6),EA(100),EI(100),X(100),Y(100),FJ(2,50),FF(4,100)
REAL*8 KE(6,6),KD(6,6),T(6,6),FP(300),KB(200,20),F(6),FO(6),D(6),BL,SI,CO,S,C
OPEN(5,FILE='PMGJ.TXT')
OPEN(6,FILE='PMGJ.OUT')
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,*)((FJ(I,J),I=1,2),J=1,NPJ)
IF(NPF.NE.0)READ(5,*)((FF(I,J),I=1,4),J=1,NPF)
WRITE(*,10)NE,NJ,N,NW,NPJ,NPF
WRITE(6,10)NE,NJ,N,NW,NPJ,NPF
WRITE(*,20)(J,X(J),Y(J),(JN(I,J),I=1,3),J=1,NJ)
WRITE(6,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)
WRITE(6,30)(J,(JE(I,J),I=1,2),EA(J),EI(J),J=1,NE)
IF(NPJ.NE.0)WRITE(*,40)((FJ(I,J),I=1,2),J=1,NPJ)
IF(NPJ.NE.0)WRITE(6,40)((FJ(I,J),I=1,2),J=1,NPJ)
IF(NPJ.NE.0)WRITE(*,50)((FJ(I,J),I=1,4),J=1,NPF)
IF(NPJ.NE.0)WRITE(6,50)((FJ(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,'FX-FY-FM'/(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 FP(I)=0.D0
IF(NPJ.EQ.0)GO TO 65
DO 60 I=1,NPJ
L=FJ(1,I)
60 FP(L)=FJ(2,I)
65 IF(NPF.EQ.0)GO TO 90
DO 70 I=1,NPF
M=FF(1,I)
CALL SCL(M,NE,NJ,BL,SI,CO,JE,X,Y)
CALL EFX(I,M,NE,NPF,BL,FF,FO,EA,EI)
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)*FO(K)
F(L)=S
75 CONTINUE
DO 85 J=1,6
L=JC(J)
IF(L.EQ.0)GO TO 85
FP(L)=FP(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 ESM(M,NE,BL,EA,EI,KD)
CALL CTM(SI,CO,T)
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 FP(I)=FP(I)-C*FP(K)
140 CONTINUE
FP(N)=FP(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 FP(I)=FP(I)-KB(I,J)*FP(L)
155 FP(I)=FP(I)/KB(I,1)
WRITE(*,165)
WRITE(6,165)
165 FORMAT(/7X,'NODE',10X,'U',14X,'V',11X,'FEI')
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)=FP(L)
175 CONTINUE
WRITE(*,180)I,D(1),D(2),D(3)
WRITE(6,180)I,D(1),D(2),D(3)
180 FORMAT(1X,I10,3E15.6)
170 CONTINUE
WRITE(*,200)
WRITE(6,200)
200 FORMAT(/4X,'ELEMENT',12X,'FN',16X,'FS',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)=FP(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=FF(1,I)
IF(M.NE.L)GO TO 250
CALL EFX(I,M,NE,NPF,BL,FF,FO,EA,EI)
DO 260 J=1,6
260 F(J)=F(J)+FO(J)
250 CONTINUE
270 WRITE(*,280)M,(F(I),I=1,6)
WRITE(6,280)M,(F(I),I=1,6)
280 FORMAT(/1X,I10,2X,'FN1=',F12.4,2X,'FS1=',F12.4,3X,'M1=',F12.4/13X,'FN2',F12.4,2X,'FS2=',F12.4,3X,'M2=',F12.4)
205 CONTINUE
CLOSE(5)
CLOSE(6)
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
10 JC(I)=JN(I,J1)
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,M,NE,NPF,BL,FF,FO,EA,EI)
DIMENSION FF(4,NPF),EA(NE),EI(NE)
REAL*8 FO(6),A,B,C,G,Q,S,BL
IND=FF(2,1)
A=FF(3,I)
Q=FF(4,I)
C=A/BL
G=C*C
B=BL-A
DO 5 J=1,6
5 FO(J)=0.D0
GO TO(10,20,30,40,50,60,70,75,80,85,90,95),IND
10 S=Q*A*0.5D0
FO(2)=-S*(2.D0*G+C*G)
FO(5)=-S*G*(2.D0-C)
S=S*A/6.D0
FO(3)=S*(6.D0-8.D0*C+3.D0*G)
FO(6)=-S*C*(4.D0-3.D0*C)
GO TO 100
20 S=B/BL
FO(2)=-Q*S*S*(1.D0+2.D0*C)
FO(5)=-Q*G*(1.D0+2.D0*S)
FO(3)=Q*S*S*A
FO(6)=-Q*B*G
GO TO 100
30 S=B/BL
FO(2)=-6.D0*Q*C*S/BL
FO(5)=-FO(2)
FO(3)=Q*S*(2.D0-3.D0*S)
FO(6)=Q*C*(2.D0-3.D0*C)
GO TO 100
40 S=Q*A*0.25D0
FO(2)=-S*(2.D0-3.D0*G+1.6D0*G*C)
FO(5)=-S*G*(3.D0-1.6D0*C)
S=S*A
FO(3)=S*(2.D0-3.D0*C+1.2D0*G)/1.5D0
FO(6)=-S*C*(1.D0-0.8D0*C)
GO TO 100
50 FO(1)=-Q*A*(1.D0-0.5D0*C)
FO(4)=-0.5D0*Q*C*A
GO TO 100
60 FO(1)=-Q*B/BL
FO(4)=-Q*C
GO TO 100
70 S=B/BL
FO(2)=-Q*G*(3.D0*S+C)
FO(5)=-FO(2)
S=S*B/BL
FO(3)=-Q*S*S*A
FO(6)=Q*G*B
GO TO 100
75 L1=IDINT(A)
S=Q*EA(M)/BL
FO(L1)=S
IF(L1.EQ.1)FO(4)=-S
IF(L1.EQ.4)FO(1)=-S
GO TO 100
80 L1=IDINT(A)
FO(L1)=12.D0*EI(M)*Q/(BL*BL*BL)
IF(L1.EQ.2) FO(5)=-FO(2)
IF(L1.EQ.5) FO(2)=-FO(5)
FO(3)=-0.5D0*BL*FO(2)
FO(6)=FO(3)
GO TO 100
85 L1=IDINT(A)
S=2.D0*EI(M)*Q/BL
FO(L1)=2.D0*S
IF(L1.EQ.3) FO(6)=S
IF(L1.EQ.6) FO(3)=S
FO(5)=3.D0*S/BL
FO(2)=-FO(5)
GO TO 100
90 FO(3)=2.D0*EI(M)*A*Q
FO(6)=-FO(3)
GO TO 100
95 FO(1)=EA(M)*A*Q
FO(4)=-FO(1)
100 RETURN
END