[下载]随机有限元程序及算例
本程序由FORTRAN语言编写,可以在DOS、POWER STATION等环境下运行。输入的数据及输出的结果采用数据文件格式。! 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
!4FORMAT(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
4FORMAT(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
20FORMAT(//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
! 30FORMAT(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)
! 54FORMAT(5X,'CH(I)=',4F11.3)
55READ(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)
122FORMAT(25X,'THE MEAN VALUES ===EE'//(10X,5F12.3))
123FORMAT(25X,'THE STANDARD DEVIATIONS ===SS'//(10X,5F12.3))
CALL CONDA
WRITE(6,111)
111FORMAT(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 CALLMR
CALLFORMMA(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
10JR(J,I)=1
IF(NR) 20,20,30
30READ(5,*)(JC(I),I=1,NR)
WRITE(6,90)(JC(I),I=1,NR)
90FORMAT(/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
50JR(2,J)=M
20N=0
DO 80 I=1,NP
DO 80 J=1,2
IF(JR(J,I)) 80, 80, 70
70N=N+1
JR(J,I)=N
80CONTINUE
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)
7FORMAT (25X,'DISTRIBUTIONPARAMETER**********'//(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
120FORMAT(10X,'CONTRAL ERROR **EOB=',F13.5)
CALL COV (NV,A,SQ,EOB)
130DO 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
95CONTINUE
MX=MX+1
NH=MA(N)
WRITE(6,115)N,NH,MX
115FORMAT(/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)
47EP(I)=EE(I)
GOTO 333
111WRITE(6,222)
222FORMAT(/20X,'TOTAL STORAGE IS GREATER THAN 100000')
STOP
333DO 5 I=1,NV
SD(I)=SP(I)
ED(I)=EP(I)
J1=JA(I)
JJ=J1-2
IF (JJ) 5,3,4
3VV=(SS(I)/EE(I))**2
SL=ALOG(1+VV)
IF (XA(I).GE.0.0) GOTO 6
XA(I)=ABS(XA(I))
6SP(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
5CONTINUE
DO 14 I=1,NM
! K=2*(I+1)
K=3*(I-1)+4
AE(1,I)=XA(K)
14AE(3,I)=XA(K-1)
JDM=JDM+1
IF (JDM.GT.6) GOTO 75
WRITE(6,23) JDM
23FORMAT(10X,'ITERATE NUMBER',5X,'JDM=',I5)
WRITE(6,37)(SD(I),I=1,NV)
WRITE(6,40)(ED(I),I=1,NV)
37FORMAT(25X,'THE LAST STANDARD DEVIATIONS ---SD'//(10X,5F12.3))
40FORMAT(25X,'THE LAST MEAN VALUES -----------ED'//(10X,5F12.3))
WRITE(6,39)(SP(I),I=1,NV)
39FORMAT(25X,'THE NEXT STANDARD DEVIATIONS----SP'//(10X,5F12.3))
WRITE(6,41)(EP(I),I=1,NV)
41FORMAT(25X,'THE NEXT MEAN VALUES------------EP'//(10X,5F12.3))
WRITE(6,38)AE
38FORMAT(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
85CALL ERFAC(AE,X,Y,MEO,NM,NP,NE)
65CALL 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
75RETURN
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)
!2FORMAT(/5X,5I5,5X,4F13.2)
S=(BI(2)*CI(3)-CI(2)*BI(3))/2.0
IF(S) 40,40,50
40WRITE(6,222) JJ
222FORMAT(/20X,'ZERO OR NEGATIVE AREA',10X,'ELEMENT NUMBER',I5)
STOP
50EO=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
10SK(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
30CONTINUE
DO 40 I=1,3
J2=ME(I)
DO 40 J=1,2
J3=2*(I-1)+J
NN(J3)=JR(J,J2)
40CONTINUE
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)
60CONTINUE
20CONTINUE
WRITE(6,555)NN1
555FORMAT(5X,'NN1=',I3)
IF(NN1.EQ.0) GOTO 100
WRITE(6,444)(SK(I),I=1,NH)
444FORMAT(5X,'SK=',5F13.4)
100CONTINUE
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
10R(I)=0.0
IF(NG) 70,70,30
30DO 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
40R(J3)=R(J3)-T*W*S/3.0
50CONTINUE
20CONTINUE
70IF(NL) 200,200,80
80READ(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)
25FORMAT(//20X,'NODES OF APPLIED LOAD**NF='//(10X,10I6))
35FORMAT(//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)
100CONTINUE
200IF(NQ) 90,90,210
210DO 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
220DO 218 I=1,3
J2=ME(I)
J3=JR(1,J2)
IF(J3) 218,218,222
222R(J3)=R(J3)+0.643*T*S/3.0
218CONTINUE
GOTO 250
230DO 228 I=1,3
J2=ME(I)
J3=JR(1,J2)
IF(J3) 228,228,232
232R(J3)=R(J3)+0.421*T*S/3.0
228CONTINUE
GOTO 250
215IF(Y(JK)-38.0) 240,260,270
270DO 275 I=1,3
J2=ME(I)
J3=JR(1,J2)
IF(J3) 275,275,280
280R(J3)=R(J3)+0.272*T*S/3.0
275CONTINUE
GOTO 250
260DO 264 I=1,3
J2=ME(I)
J3=JR(1,J2)
IF(J3) 264,264,268
268R(J3)=R(J3)+0.223*T*S/3.0
264CONTINUE
GOTO 250
240IF(Y(JK)-24.0) 290,310,320
320DO 325 I=1,3
J2=ME(I)
J3=JR(1,J2)
IF(J3) 325,325,328
328R(J3)=R(J3)+0.186*T*S/3.0
325CONTINUE
GOTO 250
310DO 315 I=1,3
J2=ME(I)
J3=JR(1,J2)
IF(J3) 315,315,318
318R(J3)=R(J3)+0.148*T*S/3.0
315CONTINUE
GOTO 250
290DO 330 I=1,3
J2=ME(I)
J3=JR(1,J2)
IF(J3) 330,330,335
335R(J3)=R(J3)+0.124*T*S/3.0
330CONTINUE
250CONTINUE
90IF(NA) 120,120,130
130DO 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
135C1=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
132IF (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
133C6=XA(1)-Y(IK)
C7=X(IK)-X(JK)
IF(LK.EQ.0) GOTO 134
R(LK)=R(LK)-0.5*C6*C7
134IF(KK.EQ.0) GOTO 150
R(KK)=R(KK)-0.5*C6*C7
GOTO 150
140K1=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
147A1=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
145IF (Y(K1).EQ.Y(K2)) GOTO 146
B1=XA(2)-Y(K1)
B2=XA(2)-Y(K2)
B3=Y(K2)-Y(K1)
R(K3)=R(K3)-(B2/6.0+B1/3.0)*B3
R(K4)=R(K4)-(B2/6.0+B1/3.0)*B3*TA2
R(K5)=R(K5)-(B1/6.0+B2/3.0)*B3
R(K6)=R(K6)-(B1/6.0+B2/3.0)*B3*TA2
GOTO 150
146E1=XA(2)-Y(K1)
E2=X(K1)-X(K2)
IF(K4.EQ.0) GOTO 148
R(K4)=R(K4)-0.5*E1*E2
148IF(K6.EQ.0) GOTO 150
R(K6)=R(K6)-0.5*E1*E2
GOTO 150
160IF(NB.EQ.3) GOTO 150
! 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 170
! IF(XA(1).LT.Y(IK).AND.XA(1).GT.Y(JK)) GOTO 180
! GOTO 150
! 170D1=XA(1)-Y(IK)
! D2=Y(IK)-Y(JK)
! D3=XA(1)-Y(JK)
! R(MK)=R(MK)+(D1/3.0+D3/6.0)*D2
! R(LK)=R(LK)-(D1/3.0+D3/6.0)*D2*TA1
! R(NK)=R(NK)+(D1/6.0+D3/3.0)*D2
! R(KK)=R(KK)-(D1/6.0+D3/3.0)*D2*TA1
! GOTO 150
! 180D4=XA(1)-Y(JK)
! D5=Y(IK)-Y(JK)
! D6=D4*D4*D4
! D7=D4*D4*D5
! R(MK)=R(MK)+D6/6.0/D5
! R(LK)=R(LK)-D6*TA1/6.0/D5
! R(NK)=R(NK)+(3.0*D7-D6)/6.0/D5
! R(KK)=R(KK)-(3.0*D7-D6)*TA1/6.0/D5
150CONTINUE
120RETURN
END
! ******************************************************************
SUBROUTINE OUTPUT
DIMENSION JR(2,800),R(1600)
COMMON/CA/NP,NE,NM,NR,NI,NL,NG,ND,NC/CC/N,NH,JR,R
DO 100 I=1,NP
L=JR(1,I)
IF(L) 10,20,30
30S=R(L)
GOTO 10
20S=0.0
10L=JR(2,I)
IF(L) 40,50,60
60SS=R(L)
GOTO 40
50SS=0.0
40WRITE(6,75) I,S,SS
75FORMAT(8X,I4,5X,F20.8,5X,F20.8)
100CONTINUE
RETURN
END
! *****************************************************************
SUBROUTINE DECOMP(SK,MA,KH,KN)
DIMENSION SK(KH),MA(KN)
COMMON/CC/N,NH,JR(2,800),R(1600)
DO 130 I=2,N
L=MA(I-1)+I-MA(I)+1
K=I-1
DO 280 J=L,K
IF(J-L) 20,280,20
20JP=MA(I)-I+J
M=MA(J-1)+J-MA(J)+1
IF(L.GT.M) M=L
MP=J-1
DO 230 IP=M,MP
JJ=MA(I)-I+IP
JK=MA(J)-J+IP
SK(JP)=SK(JP)-SK(JJ)*SK(JK)
230CONTINUE
280CONTINUE
DO 400 IP=L,K
JI=MA(I)-I+IP
JL=MA(IP)
SK(JI)=SK(JI)/SK(JL)
JN=MA(I)
SK(JN)=SK(JN)-SK(JI)*SK(JI)*SK(JL)
400CONTINUE
130CONTINUE
RETURN
END
! *****************************************************************
SUBROUTINE FOBA(SK,MA,KH,KN)
DIMENSION SK(KH),MA(KN),JR(2,800),R(1600)
COMMON/CC/N,NH,JR,R/CA/NP,NE,NM,NR,NI,NL,NG,ND,NC
DO 210 I=2,N
JJ=MA(I)
JK=I-1
JL=MA(JK)+I-JJ+1
DO 210 J=JL,JK
JP=MA(I)-I+J
R(I)=R(I)-SK(JP)*R(J)
210CONTINUE
DO 220 I=1,N
JJ=MA(I)
R(I)=R(I)/SK(JJ)
220CONTINUE
DO 230 J4=2,N
I=2+N-J4
JJ=MA(I-1)+I-MA(I)+1
M=MA(I)-I
JP=I-1
DO 240 J=JJ,JP
JM=J+M
R(J)=R(J)-SK(JM)*R(I)
240CONTINUE
230CONTINUE
RETURN
END
! ******************************************************************
SUBROUTINE CES(AE,X,Y,MEO,KM,KP,KE)
DIMENSION AE(4,KM),X(KP),Y(KP),MEO(2,KE),R(1600),JR(2,800), &
ME(3),BI(3),CI(3),B(6),CC(3)
COMMON/CA/NP,NE,NM,NR,NI,NL,NG,ND,NC
COMMON/CC/N,NH,JR,R/CB/EO,VO,W,T,S,H11,H12,H21,H22,ME,BI,CI, &
ER,TA1,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
DO 6 I=1,6
A1(I)=0.0
6CONTINUE
DO 10 IE=1,NE
CALL DIV (IE,AE,X,Y,MEO,NM,NP,NE)
V1=VO/(1.0-VO)
V2=EO/(1.0-VO*VO)
ET=V2/(1.0-V1*V1)/S/2.0
DO 55 I=1,3
J2=ME(I)
I2=JR(1,J2)
I3=JR(2,J2)
IF(I2) 50,60,70
70B(2*I-1)=R(I2)
GOTO 50
60B(2*I-1)=0.0
50IF(I3) 55,65,75
75B(2*I)=R(I3)
GOTO 55
65B(2*I)=0.0
55CONTINUE
H1=0.0
H2=0.0
H3=0.0
DO 100 I=1,3
H1=H1+BI(I)*B(2*I-1)
H2=H2+CI(I)*B(2*I)
H3=H3+BI(I)*B(2*I)+CI(I)*B(2*I-1)
100CONTINUE
AA1=ET*(H1+V1*H2)
A4=ET*(H2+V1*H1)
A3=ET*(1.0-V1)*H3/2.0
H1=AA1+A4
H2=SQRT((AA1-A4)*(AA1-A4)+4.0*A3*A3)
B(4)=(H1+H2)/2.0
B(5)=(H1-H2)/2.0
IF(ABS(A3).GT.1E-4) GOTO 400
B(6)=90.0
GOTO 450
400B(6)=ATAN((B(4)-AA1)/A3)*57.29578
450B(1)=AA1
B(2)=A4
B(3)=A3
IF (IE.EQ.2)GOTO 1001
GOTO 1002
1001KL=IE
DO 1000 I=1,6
A1(I)=B(I)
1000CONTINUE
! IF(IE.GE.KL1.AND.IE.LE.KL2) GOTO 220
! IF(IE.GE.KL1) GOTO 220
! GOTO 445
! 220IQ=IE-KL1+1
! DO 230 I=1,6
! A2(I,IQ)=B(I)
! 230CONTINUE
1002 WRITE(6,555)IE,(B(I),I=1,6)
555 FORMAT(4X,I4,2X,F8.3,2X,F8.3,2X,F8.3,3X,F10.4,2X,F10.4,2X,F8.3)
! 445CONTINUE
10CONTINUE
WRITE(6,999)KL,(A1(I),I=1,6)
999 FORMAT(4X,'KL=',I4,2X,F8.3,2X,F8.3,2X,F8.3,3X,F10.4,2X,F10.4, &
2X,F8.3)
! IF (NI.LT.0) GOTO 411
! C1=(CH(4)-CH(2))*(CH(4)-CH(3))/(CH(1)-CH(2))/(CH(1)-CH(3))
! C2=(CH(4)-CH(1))*(CH(4)-CH(3))/(CH(2)-CH(1))/(CH(2)-CH(3))
! C3=(CH(4)-CH(1))*(CH(4)-CH(2))/(CH(3)-CH(1))/(CH(3)-CH(2))
! DO 150 I=1,3
! DO 160 J=KL1,KL2,2
! P=J-KL1+1
! P=J-1+1
! L=P+1
! M=(J-KL1)/2+1
! CC(M)=(A2(I,P)+A2(I,L))/2.0
!160CONTINUE
! A1(I)=C1*CC(1)+C2*CC(2)+C3*CC(3)
!150CONTINUE
! C4=(A1(1)+A1(2))/2.0
! C5=(A1(1)-A1(2))/2.0
! A1(4)=C4+SQRT(C5*C5+A1(3)*A1(3))
! A1(5)=C4-SQRT(C5*C5+A1(3)*A1(3))
! WRITE(6,610) A1
!610FORMAT(5X,'A1(I)**=',6(2X,F8.3))
!411IP=KL2-KL1+1
! WRITE(6,421) ((A2(I,J),I=1,6),J=1,IP)
! 421FORMAT(2X,'A2(KL1---KL2 STRESSES)'/(5X,6F12.4))
RETURN
END
! *****************************************************************
SUBROUTINE ERFAC (AE,X,Y,MEO,KM,KP,KE)
DIMENSION NCI(20),NCE(4,20),ME(3),BI(3),R(1600),JR(2,800), &
AE(4,KM),X(KP),Y(KP),MEO(2,KE),CI(3)
COMMON /CA/NP,NE,NM,NR,NI,NL,NG,ND,NC,NA,NN1
COMMON /CC/N,NH,JR,R
COMMON/CB/EO,VO,W,T,S,H11,H12,H21,H22,ME,BI,CI
READ(5,*) (NCI(J),J=1,NC)
READ (5,*) ((NCE(I,J),I=1,4),J=1,NC)
WRITE(6,35)(NCI(J),J=1,NC)
WRITE(6,45)((NCE(I,J),I=1,4),J=1,NC)
35FORMAT(//30X,'NODN-NAME***NCI='//(10X,10I6))
45FORMAT(//30X,'ELEMENT-NAME***NCE='//(10X,12I5))
WRITE(6,999)
999FORMAT(//30X,'NODAL REACTIONS'//15X,'NODE',13X'X-COMP', &
13X'Y-COMP')
DO 20 JJ=1,NC
FX=0.0
FY=0.0
L=NCI(JJ)
DO 80 M=1,4
IF(NCE(M,JJ)) 80,80,70
70IE=NCE(M,JJ)
CALL DIV(IE,AE,X,Y,MEO,NM,NP,NE)
DO 100 IM=1,3
K=IM
IF(L-ME(IM)) 100,110,100
100CONTINUE
WRITE(6,555) L
555FORMAT(//20X,'ERROR OF ELEMENT MESSAGE******NODE NUMBER',I5)
110DO 400 IP=1,3
CALL KRS(BI(K),BI(IP),CI(K),CI(IP))
NL=ME(IP)
JI=JR(1,NL)
JP=JR(2,NL)
IF(JI) 210,210,220
210S=0.0
GOTO 200
220S=R(JI)
200IF(JP) 310,310,320
310SS=0.0
GOTO 300
320SS=R(JP)
300FX=FX+H11*S+H12*SS
FY=FY+H21*S+H22*SS
400CONTINUE
80CONTINUE
WRITE(6,444) L,FX,FY
444FORMAT(14X,I6,10X,F12.4,10X,F12.4)
20CONTINUE
RETURN
END
! ****************************************************************
SUBROUTINE KRS(BR,BS,CR,CS)
COMMON/CB/EO,VO,W,T,S,H11,H12,H21,H22,ME(3),BI(3),CI(3)
ET=EO*T*(1.0-VO)/4.0/S/(1.0+VO)/(1.0-2.0*VO)
V1=VO/(1.0-VO)
V2=(1.0-2.0*VO)/2.0/(1.0-VO)
H11=ET*(BR*BS+V2*CR*CS)
H12=ET*(V1*BR*CS+V2*BS*CR)
H21=ET*(V1*CR*BS+V2*BR*CS)
H22=ET*(CR*CS+V2*BR*BS)
RETURN
END
! ****************************************************************
SUBROUTINE TREAT(SK,MA,KH,KN)
DIMENSION SK(KH),MA(KN),R(1600),NDI(10),DV(2,10),JR(2,800)
COMMON/CA/NP,NE,NM,NR,NI,NL,NG,ND,NC
COMMON/CC/N,NH,JR,R
READ(5,*)(NDI(J),J=1,ND)
READ(5,*)((DV(I,J),I=1,2),J=1,ND)
WRITE(6,35)(NDI(J),J=1,ND)
WRITE(6,45)((DV(I,J),I=1,2),J=1,ND)
35FORMAT(//25X,'NODE-NAME**NDI='//(10X,10I6))
45FORMAT(//20X,'DISPLACEMENT-VALUES**DV'//(10X,6F10.4))
DO 120 I=1,ND
DO 120 J=1,2
IF(DV(I,J)) 70,120,70
70JJ=NDI(I)
L=JR(J,JJ)
IF(L.EQ.0) GOTO 120
JN=MA(L)
SK(JN)=1E15
R(L)=DV(I,J)*1E15
120CONTINUE
RETURN
END
! **************************************************************
SUBROUTINE RIC(GY,AE,X,Y,MEO,MA,SK,KH,KN,KP,KE,KM,RD,R1,KV, &
A,SQ,AS)
DIMENSION FM(6,20),RD(KN,KV),FK(6),R1(KN), &
YA(20),GA(20),A(KV,KV),SQ(KV,KV),AS(KV,KV), &
XA(20),AE(4,KM),X(KP),Y(KP),MEO(2,KE),BI(3), &
CI(3),B(6,6),NN(6),SK(KH),MA(KN)
COMMON/CA/NP,NE,NM,NR,NI,NL,NG,ND,NC,NA,NN1,DH,NV,NS
COMMON/CC/N,NH,JR(2,800),R(1600)
COMMON/CB/EO,VO,W,T,S,H11,H12,H21,H22,ME(3),BI,CI,ER,TA1, &
TA2,NB,L
COMMON/CD/EE(20),SS(20),A1(6),XA,YA,JA(20),SD(20), &
ED(20),SP(20),EP(20),A2(6,50),KL1,KL2,CH(4),KL
DO 3 I=1,N
DO 3 J=1,NV
RD(I,J)=0.0
3 CONTINUE
DO 10 IE=1,NE
DO 5 I=1,6
DO 5 J=1,NV
FM(I,J)=0.0
5 CONTINUE
CALL DIV(IE,AE,X,Y,MEO,NM,NP,NE)
LL=2*L+1
FM(2,LL)=-S*T/3.0
FM(4,LL)=FM(2,LL)
FM(6,LL)=FM(2,LL)
J1=ME(1)
J2=ME(2)
Y1=Y(J1)
Y2=Y(J2)
X1=X(J1)
X2=X(J2)
NQ=NB-2
IF (NQ) 11,12,13
11IF (XA(1).GT.Y1) GOTO 15
IF (XA(1).EQ.Y1) GOTO 38
IF (XA(1).LT.Y1.AND.XA(1).GT.Y2) GOTO 20
! GOTO 40
GOTO 13
15IF(Y1.EQ.Y2) GOTO 35
! FM(1,1)=(Y1-Y2)/2.0
! FM(3,1)=(Y1-Y2)/2.0
FM(1,1)=(Y1-Y2)/2.0
FM(3,1)=(Y1-Y2)/2.0
FM(2,3)=-S*T/3.0
FM(4,3)=FM(2,3)
FM(6,3)=FM(2,3)
GOTO 40
!35FM(2,1)=(X2-X1)*0.5
! FM(4,1)=FM(2,1)
35FM(2,1)=(X2-X1)*0.5
FM(4,1)=FM(2,1)
FM(2,6)=-S*T/3.0
FM(4,6)=FM(2,6)
FM(6,6)=FM(2,6)
GOTO 40
!38FM(1,1)=(Y1-Y2)/6.0
! FM(3,1)=(Y1-Y2)/3.0
38FM(1,1)=(Y1-Y2)/6.0
FM(3,1)=(Y1-Y2)/3.0
FM(2,3)=-S*T/3.0
FM(4,3)=FM(2,3)
FM(6,3)=FM(2,3)
GOTO 40
!20FM(1,1)=3.*(XA(1)-Y2)**2/6./(Y1-Y2)
! FM(3,1)=(6.*(XA(1)-Y2)*(Y1-Y2)-3.*(XA(1)-Y2)**2)/6./(Y1-Y2)
20FM(1,1)=3.*(XA(1)-Y2)**2/6./(Y1-Y2)
FM(3,1)=(6.*(XA(1)-Y2)*(Y1-Y2)-3.*(XA(1)-Y2)**2)/6./(Y1-Y2)
FM(2,3)=-S*T/3.0
FM(4,3)=FM(2,3)
FM(6,3)=FM(2,3)
GOTO 40
12IF (XA(2).GT.Y2) GOTO 22
IF (XA(2).EQ.Y2) GOTO 26
IF (XA(2).LT.Y2.AND.XA(2).GT.Y1) GOTO 24
! GOTO 40
GOTO 13
22IF (Y1.EQ.Y2) GOTO 23
! FM(1,2)=-(Y2-Y1)/2.0
! FM(2,2)=FM(1,2)*TA2
! FM(3,2)=FM(1,2)
! FM(4,2)=FM(2,2)
FM(1,2)=-(Y2-Y1)/2.0
FM(2,2)=FM(1,2)*TA2
FM(3,2)=FM(1,2)
FM(4,2)=FM(2,2)
FM(2,3)=-S*T/3.0
FM(4,3)=FM(2,3)
FM(6,3)=FM(2,3)
GOTO 40
!23FM(2,2)=0.5*(X2-X1)
! FM(4,2)=FM(2,2)
23FM(2,2)=0.5*(X2-X1)
FM(4,2)=FM(2,2)
FM(2,6)=-S*T/3.0
FM(4,6)=FM(2,6)
FM(6,6)=FM(2,6)
GOTO 40
!26FM(1,2)=-(Y2-Y1)/3.0
! FM(2,2)=FM(1,2)*TA2
! FM(3,2)=-(Y2-Y1)/6.0
! FM(4,2)=FM(3,2)*TA2
26FM(1,2)=-(Y2-Y1)/3.0
FM(2,2)=FM(1,2)*TA2
FM(3,2)=-(Y2-Y1)/6.0
FM(4,2)=FM(3,2)*TA2
FM(2,3)=-S*T/3.0
FM(4,3)=FM(2,3)
FM(6,3)=FM(2,3)
GOTO 40
!24FM(1,2)=-(6.*(XA(2)-Y1)*(Y2-Y1)-3.*(XA(2)-Y1)**2)/6./(Y2-Y1)
! FM(2,2)=-(6.*(XA(2)-Y1)*(Y2-Y1)-3.*(XA(2)-Y1)**2)*TA2/6./(Y2-Y1)
! FM(3,2)=-3.*(XA(2)-Y1)**2/6./(Y2-Y1)
! FM(4,2)=-3.*(XA(2)-Y1)**2*TA2/6./(Y2-Y1)
24FM(1,2)=-(6.*(XA(2)-Y1)*(Y2-Y1)-3.*(XA(2)-Y1)**2)/6./(Y2-Y1)
FM(2,2)=-(6.*(XA(2)-Y1)*(Y2-Y1)-3.*(XA(2)-Y1)**2)*TA2/6./(Y2-Y1)
FM(3,2)=-3.*(XA(2)-Y1)**2/6./(Y2-Y1)
FM(4,2)=-3.*(XA(2)-Y1)**2*TA2/6./(Y2-Y1)
FM(2,3)=-S*T/3.0
FM(4,3)=FM(2,3)
FM(6,3)=FM(2,3)
GOTO 40
13FM(2,3)=-S*T/3.0
FM(4,3)=FM(2,3)
FM(6,3)=FM(2,3)
IF(NB.EQ.3) GOTO 40
! IF(XA(1).GT.Y1) GOTO 19
! IF(XA(1).EQ.Y1) GOTO 21
! IF(XA(1).LT.Y1.AND.XA(1).GT.Y2) GOTO 33
! GOTO 40
!19FM(1,1)=0.5*(Y1-Y2)
! FM(2,1)=-FM(1,1)*TA1
! FM(2,1)=-FM(1,1)
! FM(3,1)=FM(1,1)
! FM(4,1)=FM(2,1)
! GOTO 40
!21FM(1,1)=(Y1-Y2)/6.
! FM(2,1)=-FM(1,1)*TA1
! FM(2,1)=-FM(1,1)
! FM(3,1)=FM(1,1)*2.
! FM(4,1)=FM(2,1)*2.
! GOTO 40
!33D8=XA(1)-Y2
! D9=Y1-Y2
! FM(1,1)=D8*D8/D9/2.
! FM(2,1)=-FM(1,1)*TA1
! FM(3,1)=(6.*D9*D8-3.*D8*D8)/D9/6.
! FM(4,1)=-FM(3,1)*TA1
40CONTINUE
! WRITE(6,16)((FM(I,J),J=1,NV),I=1,6)
!16FORMAT(25X,'FM ********=='//(10X,5F12.4))
DO 45 K=1,NV
DO 45 I=1,3
J2=ME(I)
DO 45 J3=1,2
K1=2*(I-1)+J3
J4=JR(J3,J2)
IF (J4.GT.0) RD(J4,K)=RD(J4,K)+FM(K1,K)
45CONTINUE
! LL=2*(L+1)
LL=3*(L-1)+4
DO 50 I=1,3
DO 50 J=1,3
CALL KRS(BI(I),BI(J),CI(I),CI(J))
B(2*I-1,2*J-1)=H11/XA(LL)
B(2*I-1,2*J)=H12/XA(LL)
B(2*I,2*J-1)=H21/XA(LL)
B(2*I,2*J)=H22/XA(LL)
50CONTINUE
DO 52 I=1,6
FK(I)=0.0
52CONTINUE
DO 55 I=1,3
J2=ME(I)
DO 55 J=1,2
J3=2*(I-1)+J
NN(J3)=JR(J,J2)
55CONTINUE
DO 60 I=1,6
DO 60 J=1,6
NJ=NN(J)
IF (NJ.EQ.0) GOTO 60
FK(I)=FK(I)+B(I,J)*R(NJ)
60CONTINUE
! WRITE(6,61) (FK(I),I=1,6)
!61FORMAT (25X,'FK**='//5X,3E16.6)
DO 65 I=1,3
J2=ME(I)
DO 65 J=1,2
K1=2*(I-1)+J
NJ=JR(J,J2)
IF (NJ.GT.0) RD(NJ,LL)=RD(NJ,LL)-FK(K1)
65CONTINUE
10CONTINUE
! WRITE(6,64) ((RD(I,J),J=1,NV),I=1,N)
!64FORMAT(25X,'RD********='//(5X,5E14.6))
DO 18 I=1,N
R1(I)=R(I)
18CONTINUE
DO 70 I=1,NV
DO 75 J=1,N
R(J)=RD(J,I)
75CONTINUE
CALL FOBA(SK,MA,NH,N)
DO 80 J=1,N
RD(J,I)=R(J)
80CONTINUE
70CONTINUE
! WRITE(6,71)((RD(I,J),J=1,NV),I=1,N)
!71FORMAT(25X,'RD*******J='//(5X,5E14.6))
IF(NI.LT.0) GOTO 444
CALL RI1(GY,AE,X,Y,MEO,RD,GA,A,SQ,AS,N,NP,NE,NM,NV)
GOTO 555
444CALL RI2(GY,AE,X,Y,MEO,RD,GA,A,SQ,AS,N,NP,NE,NM,NV)
555BB=0.0
DO 175 I=1,NV
BB=BB+YA(I)*YA(I)
175CONTINUE
BB=SQRT(BB)
WRITE(6,177) BB
177FORMAT(25X,'RELIABILITY-INDEX'/25X,'*****************'/25X, &
'BB=',F16.8)
WRITE(6,185) (XA(I),I=1,NV)
185FORMAT(10X,'XXXXX'//(10X,5F13.4))
! IF (XA(1).LE.DH) GOTO 182
! XA(1)=DH
! WRITE(6,181) XA(1)
! 181FORMAT(10X,'XA(1)=DH***',F15.5)
WRITE(6,190) (YA(I),I=1,NV)
190FORMAT(10X,'YYYYY'//(10X,5F13.4))
WRITE(6,192) (GA(I),I=1,NV)
192FORMAT(10X,'GGGAAA'//(10X,5F13.4))
WRITE(6,195) GY
195FORMAT(20X,'GY=',F13.6)
RETURN
END
!****************************************************************
SUBROUTINE RI1(GY,AE,X,Y,MEO,RD,GA,A,SQ,AS,KN,KP,KE,KM,KV)
DIMENSION AE(4,KM),X(KP),Y(KP),MEO(2,KE),RD(KN,KV),A(KV,KV), &
SQ(KV,KV),AS(KV,KV),RP(20),RE(20,6),S1(6,3),RS(20,3), &
DCG(3),RC(20),GA(KV),XA(20),YA(20),BI(3),CI(3),B1(3), &
RSS(20,3)
COMMON/CA/NP,NE,NM,NR,NI,NL,NG,ND,NC,NA,NN1,DH,NV,NS
COMMON/CC/N,NH,JR(2,800),R(1600)
COMMON/CB/EO,VO,W,T,S,H11,H12,H21,H22,ME(3),BI,CI,ER,TA1, &
TA2,NB,L
COMMON/CD/EE(20),SS(20),A1(6),XA,YA,JA(20),SD(20), &
ED(20),SP(20),EP(20),A2(6,50),KL1,KL2,CH(4),KL
! B1(1)=(CH(4)-CH(2))*(CH(4)-CH(3))/(CH(1)-CH(2))/(CH(1)-CH(3))
! B1(2)=(CH(4)-CH(1))*(CH(4)-CH(3))/(CH(2)-CH(1))/(CH(2)-CH(3))
! B1(3)=(CH(4)-CH(1))*(CH(4)-CH(2))/(CH(3)-CH(1))/(CH(3)-CH(2))
DO 20 I=1,NV
DO 20 J=1,3
RS(I,J)=0.0
20CONTINUE
! DO 115 IE=KL1,KL2
! JJ=(IE-KL1)/2+1
! JJ=(IE-1)/2+1
! CALL DIV(IE,AE,X,Y,MEO,NM,NP,NE)
CALL DIV(2,AE,X,Y,MEO,NM,NP,NE)
DO 85 J1=1,NV
DO 90 I=1,3
J2=ME(I)
DO 90 K=1,2
K1=2*(I-1)+K
NJ=JR(K,J2)
IF(NJ.GT.0)GOTO 92
GOTO 94
92RE(J1,K1)=RD(NJ,J1)
GOTO 90
94RE(J1,K1)=0.0
90CONTINUE
85CONTINUE
WRITE(6,82)((RE(I,J),J=1,6),I=1,NV)
82FORMAT(25X,'RE(I,J)*****='//(5X,5E14.6))
V1=VO/(1.0-VO)
V2=EO/(1.0-VO*VO)
A11=V2/2.0/S/(1.0-V1*V1)
A22=(1.0-V1)/2.0
DO 100 J=1,3
K1=2*J-1
K2=2*J
S1(K1,1)=A11*BI(J)
S1(K2,1)=A11*CI(J)*V1
S1(K1,2)=A11*BI(J)*V1
S1(K2,2)=A11*CI(J)
S1(K1,3)=A11*CI(J)*A22
S1(K2,3)=A11*BI(J)*A22
100CONTINUE
WRITE(6,102)((S1(I,J),J=1,3),I=1,6)
102FORMAT(25X,'S1********='//(5X,5E14.6))
DO 105 I=1,NV
DO 105 K=1,3
! RSS(I,K)=0.0
RS(I,K)=0.0
DO 110 J=1,6
! RSS(I,K)=RSS(I,K)+RE(I,J)*S1(J,K)
RS(I,K)=RS(I,K)+RE(I,J)*S1(J,K)
110CONTINUE
105CONTINUE
! LL=2*(L+1)
! J=IE-KL1+1
DO 125 I=1,3
! RSS(LL,I)=RSS(LL,I)+A2(I,J)/XA(LL)
RS(4,I)=RS(4,I)+A1(I)/XA(4)
125CONTINUE
! DO 128 I=1,NV
! DO 128 J=1,3
! RS(I,J)=RS(I,J)+RSS(I,J)*B1(JJ)/2.0
! 128CONTINUE
! 115CONTINUE
WRITE(6,116)((RS(I,J),J=1,3),I=1,NV)
116FORMAT(25X,'RS**='//(5X,5F12.4))
IF (NI.GT.0) GOTO 120
DCG(1)=0.5+(A1(1)-A1(2))/4/SQRT(((A1(1)-A1(2))/2)**2+A1(3)**2)
DCG(2)=0.5+(A1(2)-A1(1))/4/SQRT(((A1(1)-A1(2))/2)**2+A1(3)**2)
DCG(3)=A1(3)/SQRT(((A1(1)-A1(2))/2)**2+A1(3)**2)
GOTO 130
120DCG(1)=0.5-(A1(1)-A1(2))/4/SQRT(((A1(1)-A1(2))/2)**2+A1(3)**2)
DCG(2)=0.5-(A1(2)-A1(1))/4/SQRT(((A1(1)-A1(2))/2)**2+A1(3)**2)
DCG(3)=-A1(3)/SQRT(((A1(1)-A1(2))/2)**2+A1(3)**2)
130DO 145 I=1,NV
RC(I)=0.0
DO 150 J=1,3
RC(I)=RC(I)+RS(I,J)*DCG(J)
150CONTINUE
145CONTINUE
IF (NI.GT.0) GOTO 154
DO 155 I=1,NV
RC(I)=-RC(I)
155CONTINUE
! RC(NV)=1.0+RC(NV)
RC(5)=1.0+RC(5)
GOTO 158
! 154RC(NV)=RC(NV)+1.0
154RC(5)=RC(5)+1.0
158WRITE(6,156)(RC(I),I=1,NV)
156FORMAT(2X,'PARSURE XX'/(5X,5F12.4))
GG=0.0
IF(NI.GT.0) GOTO 162
GY=XA(NV)-A1(4)
GOTO 161
! 162GY=XA(NV)+A1(5)
162GY=XA(5)+A1(5)
161IF(NS.NE.0) GOTO 400
CALL RI3(GY,RC,GA)
GOTO 440
400DO 160 I=1,NV
RC(I)=RC(I)*SP(I)
160CONTINUE
DO 200 I=1,NV
DO 200 J=1,NV
AS(I,J)=0.0
AS(I,J)=SQRT(A(I,I))*SQ(J,I)
200 CONTINUE
DO 220 I=1,NV
RP(I)=0.0
220 CONTINUE
DO 300 I=1,NV
DO 300 J=1,NV
RP(I)=RP(I)+AS(I,J)*RC(J)
300 CONTINUE
WRITE(6,301)(RP(I),I=1,NV)
301 FORMAT(2X,'PARSURE ZZ'/(5X,5F12.4))
DO 330 I=1,NV
GG=GG+RP(I)*RP(I)
330 CONTINUE
GG=SQRT(GG)
DO 164 I=1,NV
GA(I)=-RP(I)/GG
YA(I)=SD(I)*YA(I)/SP(I)+(ED(I)-EP(I))/SP(I)
164CONTINUE
WRITE(6,190) (YA(I),I=1,NV)
190FORMAT(10X,'YYYYY'//(10X,5F13.4))
Y1=0.0
DO 165 I=1,NV
Y1=Y1+YA(I)*GA(I)
165CONTINUE
Y1=Y1+GY/GG
DO 170 I=1,NV
YA(I)=GA(I)*Y1
170CONTINUE
DO 350 I=1,NV
RC(I)=0.0
DO 350 J=1,NV
RC(I)=RC(I)+AS(J,I)*YA(J)
350CONTINUE
DO 180 I=1,NV
XA(I)=RC(I)*SP(I)+EP(I)
180CONTINUE
440CONTINUE
RETURN
END
! ****************************************************************
SUBROUTINE RI2(GY,AE,X,Y,MEO,RD,GA,A,SQ,AS,KN,KP,KE,KM,KV)
DIMENSION AE(4,KM),X(KP),Y(KP),MEO(2,KE),RD(KN,KV),A(KV,KV), &
SQ(KV,KV),AS(KV,KV),RP(20),RE(20,6),S1(6,3),RS(20,3), &
RC(20),GA(KV),XA(20),YA(20),BI(3),CI(3)
COMMON/CA/NP,NE,NM,NR,NI,NL,NG,ND,NC,NA,NN1,DH,NV,NS
COMMON/CC/N,NH,JR(2,800),R(1600)
COMMON/CB/EO,VO,W,T,S,H11,H12,H21,H22,ME(3),BI,CI,ER,TA1, &
TA2,NB,L
COMMON/CD/EE(20),SS(20),A1(6),XA,YA,JA(20),SD(20), &
ED(20),SP(20),EP(20),A2(6,50),KL1,KL2
DO 240 I=1,NV
240RC(I)=0.0
GY=0.0
NX=NV-2
DO 260 II=KL1,KL2,2
IQ=II-KL1+1
IF(A2(2,IQ)) 50,260,260
50 CALL DIV(II,AE,X,Y,MEO,NM,NP,NE)
DO 85 J1=1,NX
DO 90 I=1,3
J2=ME(I)
DO 90 K=1,2
K1=2*(I-1)+K
NJ=JR(K,J2)
IF(NJ.GT.0) GOTO 92
GOTO 94
92RE(J1,K1)=RD(NJ,J1)
GOTO 90
94RE(J1,K1)=0.0
90CONTINUE
85CONTINUE
! WRITE(6,82) II,((RE(I,J),J=1,6),I=1,NX)
! 82FORMAT(20X,'II=',I3,5X,'RE(I,J)====='//(5X,5E14.6))
V1=VO/(1.0-VO)
V2=EO/(1.0-VO*VO)
A11=V2/2.0/S/(1.0-V1*V1)
A22=(1.0-V1)/2.0
DO 100 J=1,3
K1=2*J-1
K2=2*J
S1(K1,1)=A11*BI(J)
S1(K2,1)=A11*CI(J)*V1
S1(K1,2)=A11*BI(J)*V1
S1(K2,2)=A11*CI(J)
S1(K1,3)=A11*CI(J)*A22
S1(K2,3)=A11*BI(J)*A22
100CONTINUE
! WRITE(6,102)((S1(I,J),J=1,3),I=1,6)
!102FORMAT(25X,'S1************='//(5X,5F12.4))
DO 105 I=1,NX
DO 105 K=1,3
RS(I,K)=0.0
DO 110 J=1,6
RS(I,K)=RS(I,K)+RE(I,J)*S1(J,K)
110CONTINUE
105CONTINUE
LL=2*(L+1)
DO 115 I=1,3
RS(LL,I)=RS(LL,I)+A2(I,IQ)/XA(LL)
115CONTINUE
! WRITE(6,116)((RS(I,J),J=1,3),I=1,NX)
!116FORMAT(25X,'RS(I,J)=*****'//(5X,5F12.4))
CB=ABS(CI(1))
IF(NB.NE.2) GOTO 55
CB=ABS(CI(2))
55 DO 216 I=1,NX
RC(I)=RC(I)-CB*(XA(NX+1)*RS(I,2)+RS(I,3))
216CONTINUE
RC(NX+1)=RC(NX+1)+ABS(A2(2,IQ))*CB
RC(NX+2)=RC(NX+2)+CB
GY=GY+CB*(XA(NX+1)*ABS(A2(2,IQ))+XA(NX+2)-A2(3,IQ))
! WRITE (6,266) (RC(I),I=1,NV)
!266FORMAT(2X,'PARSURE XX'/(5X,6F11.4))
260CONTINUE
IF(NS.NE.0) GOTO 122
CALL RI3(GY,RC,GA)
GOTO 550
122DO 160 I=1,NV
160RC(I)=RC(I)*SP(I)
DO 200 I=1,NV
DO 200 J=1,NV
AS(I,J)=0.0
AS(I,J)=SQRT(A(I,I))*SQ(J,I)
200CONTINUE
DO 220 I=1,NV
RP(I)=0.0
220CONTINUE
DO 300 I=1,NV
DO 300 J=1,NV
RP(I)=RP(I)+AS(I,J)*RC(J)
300CONTINUE
! WRITE(6,262) (RP(I),I=1,NV)
!262FORMAT(2X,'PARSURE ZZ'/(5X,6F11.4))
GG=0.0
DO 330 I=1,NV
GG=GG+RP(I)*RP(I)
330CONTINUE
GG=SQRT(GG)
DO 164 I=1,NV
GA(I)=-RP(I)/GG
YA(I)=SD(I)*YA(I)/SP(I)+(ED(I)-EP(I))/SP(I)
164CONTINUE
WRITE(6,190) (YA(I),I=1,NV)
190FORMAT(10X,'YYYYY'//(10X,5F13.4))
Y1=0.0
DO 165 I=1,NV
Y1=Y1+YA(I)*GA(I)
165CONTINUE
Y1=Y1+GY/GG
DO 170 I=1,NV
YA(I)=GA(I)*Y1
170CONTINUE
DO 350 I=1,NV
RC(I)=0.0
DO 350 J=1,NV
RC(I)=RC(I)+AS(J,I)*YA(J)
350CONTINUE
DO 180 I=1,NV
XA(I)=RC(I)*SP(I)+EP(I)
180CONTINUE
550CONTINUE
RETURN
END
! *****************************************************************
SUBROUTINE RI3(GY,RC,GA)
DIMENSION RC(20),GA(20),XA(20),YA(20)
COMMON/CA/NP,NE,NM,NR,NI,NL,NG,ND,NC,NA,NN1,DH,NV,NS
COMMON/CD/EE(20),SS(20),A1(6),XA,YA,JA(20),SD(20), &
ED(20),SP(20),EP(20)
GG=0.0
DO 160 I=1,NV
RC(I)=RC(I)*SP(I)
GG=GG+RC(I)**2
160CONTINUE
GG=SQRT(GG)
DO 164 I=1,NV
GA(I)=-RC(I)/GG
YA(I)=SD(I)*YA(I)/SP(I)+(ED(I)-EP(I))/SP(I)
164CONTINUE
WRITE(6,190)(YA(I),I=1,NV)
190FORMAT(10X,'YYYYY'//(10X,5F13.4))
Y1=0.0
DO 165 I=1,NV
Y1=Y1+YA(I)*GA(I)
165CONTINUE
Y1=Y1+GY/GG
DO 170 I=1,NV
YA(I)=GA(I)*Y1
170CONTINUE
DO 180 I=1,NV
XA(I)=YA(I)*SP(I)+EP(I)
180CONTINUE
RETURN
END
! *****************************************************************
SUBROUTINE COV(N1,A,SQ,EO)
DIMENSION A(N1,N1),SQ(N1,N1)
READ(5,*)((A(I,J),I=1,N1),J=1,N1)
WRITE(6,6)((A(I,J),I=1,N1),J=1,N1)
6FORMAT(10X,'A(I,J)='/(5X,5F13.5))
DO 10 I=1,N1
DO 10 J=1,N1
SQ(I,J)=0.0
10 CONTINUE
DO 15 I=1,N1
SQ(I,I)=1.0
15 CONTINUE
G=0.0
DO 40 I=2,N1
I1=I-1
DO 40 J=1,I1
G=G+2.0*A(I,J)**2
40 CONTINUE
! WRITE(6,42) G
!42 FORMAT(5X,'G=',F13.5)
S1=SQRT(G)
FN=FLOAT(N1)
S2=EO*S1/FN
! WRITE (6,45) S2
S3=S1
!45 FORMAT (10X,'S2=',F13.5)
L=0
50 S3=S3/FN
! WRITE (6,55) S3
!55 FORMAT (10X,'S3=',F13.5)
60 DO 130 IQ=2,N1
IQ1=IQ-1
DO 130 IP=1,IQ1
IF(ABS(A(IP,IQ)).LT.S3) GOTO 130
L=1
V1=A(IP,IP)
V2=A(IP,IQ)
V3=A(IQ,IQ)
U=0.5*(V1-V3)
! WRITE(6,65) V1,V2,V3,U
!65 FORMAT(10X,'V1*V2*V3*U=',4F13.5)
IF(U) 70,80,90
70 B=-1.0
GOTO 100
90 B=1.0
100 G=-B*V2/SQRT(V2*V2+U*U)
GOTO 105
80 IF(V2.GT.S3) GOTO 85
G=1.0
GOTO 105
85 G=-1.0
105ST=G/SQRT(2.0*(1.0+SQRT(1.0-G*G)))
CT=SQRT(1.0-ST*ST)
! WRITE(6,107) ST,CT
!107FORMAT(10X,'ST*CT=',2F13.5)
DO 110 I=1,N1
G=A(I,IP)*CT-A(I,IQ)*ST
A(I,IQ)=A(I,IP)*ST+A(I,IQ)*CT
A(I,IP)=G
G=SQ(I,IP)*CT-SQ(I,IQ)*ST
SQ(I,IQ)=SQ(I,IP)*ST+SQ(I,IQ)*CT
SQ(I,IP)=G
110CONTINUE
DO 120 I=1,N1
A(IP,I)=A(I,IP)
A(IQ,I)=A(I,IQ)
120CONTINUE
G=2.0*V2*ST*CT
A(IP,IP)=V1*CT*CT+V3*ST*ST-G
A(IQ,IQ)=V1*ST*ST+V3*CT*CT+G
A(IP,IQ)=(V1-V3)*ST*CT+V2*(CT*CT-ST*ST)
A(IQ,IP)=A(IP,IQ)
130CONTINUE
IF(L-1) 150,140,150
140L=0
GOTO 60
150WRITE(6,145)((A(I,J),I=1,N1),J=1,N1)
WRITE(6,142)((SQ(I,J),I=1,N1),J=1,N1)
142FORMAT(20X,'SQ(I,J)='/(5X,5F13.5))
145FORMAT(20X,'A(I,J)='/(5X,5F13.5))
IF(S3.GT.S2) GOTO 50
RETURN
END
例子：INPUT.dat
6,4,1,3,0,0,1,0,0,1,1,5,0
0.01,1,40
30,5,1.4,2000000,10
3,0.5,0.028,200000,2.5
2000000,0.2,1.4,1
0,0,20,0,20,40
40,20,20,0,0,0
001002,003011,002004,005011,002005,003013,006003,005012
00400,00500,00600
1,1,1,1,2
0
4,5,6
0,0,0,2,0,2,3,4,0,0,0,4 请问楼主：该随机有限元程序是不是 武清玺编著的那本书上的程序？<BR> 谢谢先，
有以下疑问：
（1）没有有关程序注释的readme文件。
（2）随时.90后缀的f90文件，却完全是f77标准的程序。
（3）敬请楼主说明程序的配套书籍^_^
再次感谢。 谢谢，能不能说得详细点儿。 就是武清玺那本书上的程序 <P>查了一下，图书馆没有；超星也没有。<BR>不知道哪位可以提供一下电子版的，<BR>非常感谢！！</P> 是啊，急需电子版的，哪位好心人做点好事啊 谢谢 汗！看来要靠自己写注释了 谢谢楼主,辛苦了 ou的神！终于弄到武老师的这个程序啦！拜谢！ 怎么运行时出现问题啊？哪位运行通过了的啊？ 好些年前的程式, 个人以为不见得编译仍相同, 可能需稍为修改下吧! 不能编译，而且是fortran77的代码，看起来太累了··
要是fortran90的就太好了 要是有电子版的就好了
页:
[1]