📄 程序.f90
字号:
! main program for SFEM
COMMON/CA/NP,NE,NM,NR,NI,NL,NG,ND,NC,NA,NN1,DH,NV,NS,NQ
COMMON/CC/N,NH,JR(2,800),R(1600)
COMMON /CB/EO,VO,W,T,S,H11,H12,H21,H22,ME(3),BI(3),CI(3),ER,TA2,NB,L
COMMON/CD/EE(20),SS(20),A1(6),XA(20),YA(20),JA(20),SD(20),ED(20),SP(20),EP(20),A2(6,50),KL1,KL2,CH(4),KL
open (5,FILE='D:\Str.及SFEM\程序及算例\INPUT.DAT',status='OLD',ACCESS='SEQUENTIAL')
open (6,FILE='D:\Str.及SFEM\程序及算例\OUTPUT.DAT',status='UNKNOWN',ACCESS='SEQUENTIAL')
READ(5,*)NP,NE,NM,NR,NI,NL,NG,ND,NC,NA,NN1,NV,NQ
! READ(5,*)ER,TA1,TA2,DH
! WRITE(6,4)ER,TA1,TA2,DH
! 4 FORMAT(5X,'ER=',F13.5,3X,'TA1=',F11.3,3X,'TA2=',F11.3,3X,'DH=',F11.3)
READ(5,*)ER,TA2,DH
WRITE(6,4)ER,TA2,DH
4 FORMAT(5X,'ER=',F13.5,3X,'TA2=',F11.3,3X,'DH=',F11.3)
WRITE(6,20)NP,NE,NM,NR,NI,NQ,NL,NG,ND,NC,NA,NV
20 FORMAT(//9X,'NP=',I5,3X,'NE=',I5,3X,'NM=',I5,3X,'NR=',I5,3X, &
'NI=',I5,3X,'NQ=',I5/9X,'NL=',I5,3X,'NG=',I5,3X, &
'ND=',I5,3X,'NC=',I5,3X,'NA=',I5,3X,'NV=',I5)
! READ(5,*) KL1,KL2
! WRITE(6,30)KL1,KL2
! 30 FORMAT(10X,'KL1=',I5,5X,'KL2=',I5)
! IF(NI.LT.0) GOTO 55
! READ(5,*)(CH(I),I=1,4)
! WRITE(6,54)(CH(I),I=1,4)
! 54 FORMAT(5X,'CH(I)=',4F11.3)
55 READ(5,*)(EE(I),I=1,NV)
READ(5,*)(SS(I),I=1,NV)
WRITE(6,122) (EE(I),I=1,NV)
WRITE(6,123) (SS(I),I=1,NV)
122 FORMAT(25X,'THE MEAN VALUES ===EE'//(10X,5F12.3))
123 FORMAT(25X,'THE STANDARD DEVIATIONS ===SS'//(10X,5F12.3))
CALL CONDA
WRITE(6,111)
111 FORMAT(5X,'END')
STOP
END
! *********************************************************
SUBROUTINE CONDA
DIMENSION X(800),Y(800),MEO(2,1000),AE(4,5),A(20,20), &
SQ(20,20),AS(20,20)
COMMON /CA/NP,NE,NM,NR,NI,NL,NG,ND,NC,NA,NN1,DH,NV
COMMON /CC/N,NH,JR(2,800),R(1600)
COMMON /CB/EO,VO,W,T,S,H11,H12,H21,H22,ME(3),BI(3),CI(3),ER, &
TA1,TA2,NB,L
COMMON /CD/EE(20),SS(20)
READ(5,*)((AE(I,J),I=1,4),J=1,NM)
READ(5,*)(X(I),I=1,NP)
READ(5,*)(Y(I),I=1,NP)
READ(5,*)((MEO(I,J),I=1,2),J=1,NE)
WRITE(6,81)((AE(I,J),I=1,4),J=1,NM)
WRITE(6,82)(X(I),I=1,NP)
WRITE(6,83)(Y(I),I=1,NP)
WRITE(6,84)((MEO(I,J),I=1,2),J=1,NE)
81 FORMAT(//25X,'EO*VO*W*T**'//(5X,4F15.4))
82 FORMAT(//20X,'X--COORDINATE'//(5X,10F7.2))
83 FORMAT(//20X,'Y--COORDINATE'//(5X,10F7.2))
84 FORMAT(//20X,'ELEMENT-DATA'//(5X,10I7))
90 CALL MR
CALL FORMMA(AE,X,Y,MEO,NM,NP,NE,NV,A,SQ,AS)
RETURN
END
! ****************************************************************
SUBROUTINE MR
DIMENSION JC(50),JR(2,800)
COMMON/CA/NP,NE,NM,NR,NI,NL,NG,ND,NC
COMMON/CC/N,NH,JR,R(1600)
DO 10 I=1,NP
DO 10 J=1,2
10 JR(J,I)=1
IF(NR) 20,20,30
30 READ(5,*)(JC(I),I=1,NR)
WRITE(6,90)(JC(I),I=1,NR)
90 FORMAT(/20X,'CONSTRINED MESSAGE-JC='//(5X,10I7))
DO 50 I=1,NR
K=JC(I)
J=K/100
L=(K-J*100)/10
M=K-J*100-L*10
JR(1,J)=L
50 JR(2,J)=M
20 N=0
DO 80 I=1,NP
DO 80 J=1,2
IF(JR(J,I)) 80, 80, 70
70 N=N+1
JR(J,I)=N
80 CONTINUE
RETURN
END
! **************************************************************
SUBROUTINE FORMMA(AE,X,Y,MEO,KM,KP,KE,KV,A,SQ,AS)
DIMENSION MA(2000),X(KP),Y(KP),MEO(2,KE),AE(4,KM),R(1600), &
JR(2,800),A(KV,KV),SQ(KV,KV),AS(KV,KV),ME(3),NN(6), &
SK(100000),RD(1600,20),R1(1600)
COMMON/CA/NP,NE,NM,NR,NI,NL,NG,ND,NC,NA,NN1,DH,NV,NS,NQ
COMMON/CC/N,NH,JR,R/CB/EO,VO,W,T,S,H11,H12,H21,H22,ME,BI(3), &
CI(3),ER,TA1,TA2,NB,L/CD/EE(20),SS(20),A1(6),XA(20), &
YA(20),JA(20),SD(20),ED(20),SP(20),EP(20),A2(6,50),KL1, &
KL2,CH(4)
! DATA (XA(I),I=1,7)/88.,15.5,2.7,1.15E6,2.85,3.2E6,200.0/
READ(5,*)(JA(I),I=1,NV)
WRITE(6,7) (JA(I),I=1,NV)
7 FORMAT (25X,'DISTRIBUTION PARAMETER**********'//(10X,10I5))
READ(5,*) NS
WRITE(6,2) NS
2 FORMAT(5X,'CORRELATED VARIABLE MESSAGE **NS=',I5)
IF(NS.EQ.0) GOTO 130
READ(5,*) EOB
WRITE(6,120) EOB
120 FORMAT(10X,'CONTRAL ERROR **EOB=',F13.5)
CALL COV (NV,A,SQ,EOB)
130 DO 8 I=1,NV
XA(I)=0.0
XA(I)=EE(I)
YA(I)=0.0
8 CONTINUE
DO 10 I=1,N
10 MA(I)=0
DO 20 IE=1,NE
CALL DIV(IE,AE,X,Y,MEO,NM,NP,NE)
DO 30 I=1,3
JB=ME(I)
DO 30 M=1,2
JJ=2*(I-1)+M
NN(JJ)=JR(M,JB)
30 CONTINUE
L=N
DO 80 I=1,6
IF(NN(I)) 80,80,60
60 IF(NN(I)-L) 70,80,80
70 L=NN(I)
80 CONTINUE
DO 15 M=1,6
JP=NN(M)
IF(JP.EQ.0) GOTO 15
IF(JP-L.GT.MA(JP))MA(JP)=JP-L
15 CONTINUE
20 CONTINUE
MX=0
MA(1)=1
DO 95 I=2,N
IF(MA(I).GT.MX) MX=MA(I)
MA(I)=MA(I)+MA(I-1)+1
95 CONTINUE
MX=MX+1
NH=MA(N)
WRITE(6,115)N,NH,MX
115 FORMAT(/22X,'T0TAL DEGREES OF FREEDOM N=',I5/25X,'TOTAL- &
STORAGE','NH=',I6/24X,'MAX-SEMI-BANDWIDTH MX=',I5)
IF(NH.GT.100000) GOTO 111
JDM=0
DO 47 I=1,NV
SP(I)=SS(I)
47 EP(I)=EE(I)
GOTO 333
111 WRITE(6,222)
222 FORMAT(/20X,'TOTAL STORAGE IS GREATER THAN 100000')
STOP
333 DO 5 I=1,NV
SD(I)=SP(I)
ED(I)=EP(I)
J1=JA(I)
JJ=J1-2
IF (JJ) 5,3,4
3 VV=(SS(I)/EE(I))**2
SL=ALOG(1+VV)
IF (XA(I).GE.0.0) GOTO 6
XA(I)=ABS(XA(I))
6 SP(I)=XA(I)*SQRT(SL)
EP(I)=XA(I)*(1.0+ALOG(EE(I))-ALOG(XA(I)*SQRT(1.0+VV)))
GOTO 5
4 AR=1.28255/SS(I)
AK=EE(I)-0.5772/AR
Q=EXP(-AR*(XA(I)-AK))
AF=EXP(-Q)
AGF=AR*Q*EXP(-Q)
AA=1.0-AF
B0=1.570796288
B1=3.706987906E-2
B2=-8.364353589E-4
B3=-2.250947176E-4
B4=6.841218299E-6
B5=5.824238515E-6
B6=-1.04527497E-6
B7=8.360937017E-8
B8=-3.231081277E-9
B9=3.657763036E-11
B10=6.936233982E-13
IF(AA-0.5) 230,240,250
230 BB=AA
GOTO 255
250 BB=1.0-AA
255 YY=-ALOG(4.0*BB*(1.0-BB))
U1=YY*(B5+YY*(B6+YY*(B7+YY*(B8+YY*(B9+YY*B10)))))
UU=YY*(B0+YY*(B1+YY*(B2+YY*(B3+YY*(B4+U1)))))
UB=SQRT(UU)
IF(AA.LT.0.5) GOTO 260
UB=-UB
GOTO 260
240 UB=0.0
260 SP(I)=EXP(-UB*UB/2.0)/SQRT(2.0*3.1416)/AGF
EP(I)=XA(I)-SP(I)*UB
5 CONTINUE
DO 14 I=1,NM
! K=2*(I+1)
K=3*(I-1)+4
AE(1,I)=XA(K)
14 AE(3,I)=XA(K-1)
JDM=JDM+1
IF (JDM.GT.6) GOTO 75
WRITE(6,23) JDM
23 FORMAT(10X,'ITERATE NUMBER',5X,'JDM=',I5)
WRITE(6,37)(SD(I),I=1,NV)
WRITE(6,40)(ED(I),I=1,NV)
37 FORMAT(25X,'THE LAST STANDARD DEVIATIONS ---SD'//(10X,5F12.3))
40 FORMAT(25X,'THE LAST MEAN VALUES -----------ED'//(10X,5F12.3))
WRITE(6,39)(SP(I),I=1,NV)
39 FORMAT(25X,'THE NEXT STANDARD DEVIATIONS----SP'//(10X,5F12.3))
WRITE(6,41)(EP(I),I=1,NV)
41 FORMAT(25X,'THE NEXT MEAN VALUES------------EP'//(10X,5F12.3))
WRITE(6,38)AE
38 FORMAT(5X,'AE='//(4F13.4))
CALL MGK(AE,X,Y,MEO,MA,SK,NM,NP,NE,N,NH)
CALL LOAD(AE,X,Y,MEO,NM,NP,NE)
WRITE(6,444)
444 FORMAT(/35X,'NODAL FORCES'//15X,'NODE',13X,'X-COMP',13X,'Y-COMP')
CALL OUTPUT
IF(ND.GT.0) CALL TREAT(SK,MA,NH,N)
CALL DECOMP(SK,MA,NH,N)
CALL FOBA(SK,MA,NH,N)
WRITE(6,555)
555 FORMAT(/35X,'NODAL DISPLACEMENTS'//15X,'NODE',13X, &
'X-COMP',13X,'Y-COMP')
CALL OUTPUT
WRITE(6,666)
666 FORMAT(//35X,'ELEMENT STRESSES'//1X,'ELEMENT',2X,'X-STRESS',2X, &
'Y-STRESS',2X,'XY-STRESS',2X,'MAX-STRESS',2X,'MIN-STRESS' &
,2X,'ANGLE'//)
CALL CES(AE,X,Y,MEO,NM,NP,NE)
IF(NC) 65,65,85
85 CALL ERFAC(AE,X,Y,MEO,NM,NP,NE)
65 CALL RIC(GY,AE,X,Y,MEO,MA,SK,NH,N,NP,NE,NM,RD,R1,NV,A,SQ,AS)
GY=ABS(GY)
IF (GY.GT.ER) GOTO 333
75 RETURN
END
! ******************************************************************
SUBROUTINE DIV(JJ,AE,X,Y,MEO,KM,KP,KE)
DIMENSION ME(3),BI(3),CI(3),X(KP),Y(KP),AE(4,KM),MEO(2,KE)
COMMON/CB/EO,VO,W,T,S,A(4),ME,BI,CI,ER,TA1,TA2,NB,L
II=MEO(1,JJ)
I=II/1000
ME(1)=I
J=II-I*1000
ME(2)=J
IJ=MEO(2,JJ)
M=IJ/1000
ME(3)=M
IK=IJ-M*1000
L=IK/10
NB=IJ-M*1000-L*10
BI(1)=Y(J)-Y(M)
CI(1)=X(M)-X(J)
BI(2)=Y(M)-Y(I)
CI(2)=X(I)-X(M)
BI(3)=Y(I)-Y(J)
CI(3)=X(J)-X(I)
! WRITE(6,2) I,J,M,L,NB,CI(2),BI(2),CI(3),BI(3)
! 2 FORMAT(/5X,5I5,5X,4F13.2)
S=(BI(2)*CI(3)-CI(2)*BI(3))/2.0
IF(S) 40,40,50
40 WRITE(6,222) JJ
222 FORMAT(/20X,'ZERO OR NEGATIVE AREA',10X,'ELEMENT NUMBER',I5)
STOP
50 EO=AE(1,L)
VO=AE(2,L)
W=AE(3,L)
T=AE(4,L)
RETURN
END
! *****************************************************************
SUBROUTINE MGK(AE,X,Y,MEO,MA,SK,KM,KP,KE,KN,KH)
DIMENSION AE(4,KM),X(KP),Y(KP),MEO(2,KE),MA(KN),SK(KH), &
JR(2,800),NN(6),ME(3),BI(3),CI(3),B(6,6)
COMMON /CA/NP,NE,NM,NR,NI,NL,NG,ND,NC,NA,NN1
COMMON /CB/EO,VO,W,T,S,H11,H12,H21,H22,ME,BI,CI,ER,TA1,TA2,NB
COMMON/CC/N,NH,JR,R(1600)/CD/EE(20),SS(20),A1(6), &
XA(20),YA(20),JA(20)
DO 10 I=1,NH
10 SK(I)=0.0
DO 20 IE=1,NE
CALL DIV(IE,AE,X,Y,MEO,NM,NP,NE)
DO 30 I=1,3
DO 30 J=1,3
CALL KRS(BI(I),BI(J),CI(I),CI(J))
B(2*I-1,2*J-1)=H11
B(2*I-1,2*J)=H12
B(2*I,2*J-1)=H21
B(2*I,2*J)=H22
30 CONTINUE
DO 40 I=1,3
J2=ME(I)
DO 40 J=1,2
J3=2*(I-1)+J
NN(J3)=JR(J,J2)
40 CONTINUE
DO 60 I=1,6
DO 60 J=1,6
IF(NN(J).EQ.0.OR.NN(I).LT.NN(J)) GOTO 60
JJ=NN(I)
JK=NN(J)
JL=MA(JJ)
JM=JJ-JK
JN=JL-JM
SK(JN)=SK(JN)+B(I,J)
60 CONTINUE
20 CONTINUE
WRITE(6,555)NN1
555 FORMAT(5X,'NN1=',I3)
IF(NN1.EQ.0) GOTO 100
WRITE(6,444)(SK(I),I=1,NH)
444 FORMAT(5X,'SK=',5F13.4)
100 CONTINUE
RETURN
END
! *****************************************************************
SUBROUTINE LOAD(AE,X,Y,MEO,KM,KP,KE)
DIMENSION AE(4,KM),X(KP),Y(KP),MEO(2,KE),ME(3),R(1600),NF(50), &
FV(2,50),JR(2,800)
COMMON/CA/NP,NE,NM,NR,NI,NL,NG,ND,NC,NA,NN1,DH,NV,NS,NQ
COMMON/CC/N,NH,JR,R/CB/EO,VO,W,T,S,A(4),ME,BI(3),CI(3),ER,TA1, &
TA2,NB,L
COMMON/CD/EE(20),SS(20),AA(6),XA(20)
DO 10 I=1,N
10 R(I)=0.0
IF(NG) 70,70,30
30 DO 20 IE=1,NE
CALL DIV(IE,AE,X,Y,MEO,NM,NP,NE)
DO 50 I=1,3
J2=ME(I)
J3=JR(2,J2)
IF(J3) 50,50,40
40 R(J3)=R(J3)-T*W*S/3.0
50 CONTINUE
20 CONTINUE
70 IF(NL) 200,200,80
80 READ(5,*)(NF(I),I=1,NL)
READ(5,*)((FV(I,J),I=1,2),J=1,NL)
WRITE(6,25)(NF(I),I=1,NL)
WRITE(6,35)((FV(I,J),I=1,2),J=1,NL)
25 FORMAT(//20X,'NODES OF APPLIED LOAD**NF='//(10X,10I6))
35 FORMAT(//20X,'LUMPED-LOADS**FV='//(10X,6F10.4))
DO 100 I=1,NL
JJ=NF(I)
J=JR(1,JJ)
M=JR(2,JJ)
IF(J.GT.0)R(J)=R(J)+FV(1,I)
IF(M.GT.0) R(M)=R(M)+FV(2,I)
100 CONTINUE
200 IF(NQ) 90,90,210
210 DO 250 IE=1,NE
CALL DIV(IE,AE,X,Y,MEO,NM,NP,NE)
IF(L.GT.2) GOTO 250
JK=ME(2)
IF(Y(JK)-62.0) 215,230,220
220 DO 218 I=1,3
J2=ME(I)
J3=JR(1,J2)
IF(J3) 218,218,222
222 R(J3)=R(J3)+0.643*T*S/3.0
218 CONTINUE
GOTO 250
230 DO 228 I=1,3
J2=ME(I)
J3=JR(1,J2)
IF(J3) 228,228,232
232 R(J3)=R(J3)+0.421*T*S/3.0
228 CONTINUE
GOTO 250
215 IF(Y(JK)-38.0) 240,260,270
270 DO 275 I=1,3
J2=ME(I)
J3=JR(1,J2)
IF(J3) 275,275,280
280 R(J3)=R(J3)+0.272*T*S/3.0
275 CONTINUE
GOTO 250
260 DO 264 I=1,3
J2=ME(I)
J3=JR(1,J2)
IF(J3) 264,264,268
268 R(J3)=R(J3)+0.223*T*S/3.0
264 CONTINUE
GOTO 250
240 IF(Y(JK)-24.0) 290,310,320
320 DO 325 I=1,3
J2=ME(I)
J3=JR(1,J2)
IF(J3) 325,325,328
328 R(J3)=R(J3)+0.186*T*S/3.0
325 CONTINUE
GOTO 250
310 DO 315 I=1,3
J2=ME(I)
J3=JR(1,J2)
IF(J3) 315,315,318
318 R(J3)=R(J3)+0.148*T*S/3.0
315 CONTINUE
GOTO 250
290 DO 330 I=1,3
J2=ME(I)
J3=JR(1,J2)
IF(J3) 330,330,335
335 R(J3)=R(J3)+0.124*T*S/3.0
330 CONTINUE
250 CONTINUE
90 IF(NA) 120,120,130
130 DO 150 IE=1,NE
CALL DIV(IE,AE,X,Y,MEO,NM,NP,NE)
IF(NB.GT.2) GOTO 160
IF(NB.EQ.2) GOTO 140
IK=ME(1)
JK=ME(2)
MK=JR(1,IK)
NK=JR(1,JK)
LK=JR(2,IK)
KK=JR(2,JK)
IF(XA(1).GE.Y(IK)) GOTO 132
IF(XA(1).LT.Y(IK).AND.XA(1).GT.Y(JK)) GOTO 135
GOTO 150
135 C1=XA(1)-Y(JK)
C2=Y(IK)-Y(JK)
R(MK)=R(MK)+C1*C1*C1/(6.0*C2)
R(NK)=R(NK)+(3.0*C1*C1*C2-C1*C1*C1)/(6.0*C2)
GOTO 150
132 IF (Y(IK).EQ.Y(JK)) GOTO 133
C3=XA(1)-Y(IK)
C4=XA(1)-Y(JK)
C5=Y(IK)-Y(JK)
R(MK)=R(MK)+(C4*C5)/6.0+(C3*C5)/3.0
R(NK)=R(NK)+(C3*C5)/6.0+(C4*C5)/3.0
GOTO 150
133 C6=XA(1)-Y(IK)
C7=X(IK)-X(JK)
IF(LK.EQ.0) GOTO 134
R(LK)=R(LK)-0.5*C6*C7
134 IF(KK.EQ.0) GOTO 150
R(KK)=R(KK)-0.5*C6*C7
GOTO 150
140 K1=ME(1)
K2=ME(2)
K3=JR(1,K1)
K4=JR(2,K1)
K5=JR(1,K2)
K6=JR(2,K2)
IF(XA(2).GE.Y(K2)) GOTO 145
IF(XA(2).GT.Y(K1).AND.XA(2).LT.Y(K2)) GOTO 147
GOTO 150
147 A1=XA(2)-Y(K1)
A2=Y(K2)-Y(K1)
A3=3*A1*A1*A2-A1*A1*A1
R(K3)=R(K3)-A3/(6.0*A2)
R(K4)=R(K4)-A3*TA2/(6.0*A2)
R(K5)=R(K5)-A1*A1*A1/(6.0*A2)
R(K6)=R(K6)-A1*A1*A1*TA2/(6.0*A2)
GOTO 150
145 IF (Y(K1).EQ.Y(K2)) GOTO 146
B1=XA(2)-Y(K1)
B2=XA(2)-Y(K2)
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -