[下载]随机有限元程序及算例

yoyoxo

本程序由FORTRAN语言编写,可以在DOS、POWER STATION等环境下运行。输入的数据及输出的结果采用数据文件格式。

1. !        main program for SFEM          
   
2.           COMMON/CA/NP,NE,NM,NR,NI,NL,NG,ND,NC,NA,NN1,DH,NV,NS,NQ
   
3.       COMMON/CC/N,NH,JR(2,800),R(1600)
   
4.       COMMON /CB/EO,VO,W,T,S,H11,H12,H21,H22,ME(3),BI(3),CI(3),ER,TA2,NB,L
   
5.       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
   
6.       open (5,FILE='D:\Str.及SFEM\程序及算例\INPUT.DAT',status='OLD',ACCESS='SEQUENTIAL')
   
7.           open (6,FILE='D:\Str.及SFEM\程序及算例\OUTPUT.DAT',status='UNKNOWN',ACCESS='SEQUENTIAL')
   
8.           READ(5,*)NP,NE,NM,NR,NI,NL,NG,ND,NC,NA,NN1,NV,NQ
   
9. !     READ(5,*)ER,TA1,TA2,DH
   
10. !     WRITE(6,4)ER,TA1,TA2,DH
    
11. !  4  FORMAT(5X,'ER=',F13.5,3X,'TA1=',F11.3,3X,'TA2=',F11.3,3X,'DH=',F11.3)
    
12.           READ(5,*)ER,TA2,DH
    
13.       WRITE(6,4)ER,TA2,DH
    
14.    4  FORMAT(5X,'ER=',F13.5,3X,'TA2=',F11.3,3X,'DH=',F11.3)
    
15.       WRITE(6,20)NP,NE,NM,NR,NI,NQ,NL,NG,ND,NC,NA,NV
    
16.   20  FORMAT(//9X,'NP=',I5,3X,'NE=',I5,3X,'NM=',I5,3X,'NR=',I5,3X, &
    
17.               'NI=',I5,3X,'NQ=',I5/9X,'NL=',I5,3X,'NG=',I5,3X, &
    
18.                         'ND=',I5,3X,'NC=',I5,3X,'NA=',I5,3X,'NV=',I5)
    
19. !     READ(5,*) KL1,KL2
    
20. !     WRITE(6,30)KL1,KL2
    
21. ! 30  FORMAT(10X,'KL1=',I5,5X,'KL2=',I5)
    
22. !          IF(NI.LT.0) GOTO 55
    
23. !     READ(5,*)(CH(I),I=1,4)
    
24. !     WRITE(6,54)(CH(I),I=1,4)
    
25. ! 54  FORMAT(5X,'CH(I)=',4F11.3)
    
26.   55  READ(5,*)(EE(I),I=1,NV)
    
27.       READ(5,*)(SS(I),I=1,NV)
    
28.       WRITE(6,122) (EE(I),I=1,NV)
    
29.       WRITE(6,123) (SS(I),I=1,NV)
    
30. 122  FORMAT(25X,'THE   MEAN   VALUES ===EE'//(10X,5F12.3))
    
31. 123  FORMAT(25X,'THE STANDARD DEVIATIONS ===SS'//(10X,5F12.3))
    
32.       CALL CONDA
    
33.       WRITE(6,111)
    
34. 111  FORMAT(5X,'END')
    
35.       STOP
    
36.       END
    
37. 
    
38.           !     *********************************************************
    
39.       SUBROUTINE CONDA
    
40.       DIMENSION X(800),Y(800),MEO(2,1000),AE(4,5),A(20,20), &
    
41.                   SQ(20,20),AS(20,20)
    
42.       COMMON /CA/NP,NE,NM,NR,NI,NL,NG,ND,NC,NA,NN1,DH,NV
    
43.       COMMON /CC/N,NH,JR(2,800),R(1600)
    
44.       COMMON /CB/EO,VO,W,T,S,H11,H12,H21,H22,ME(3),BI(3),CI(3),ER, &
    
45.                    TA1,TA2,NB,L
    
46.           COMMON /CD/EE(20),SS(20)
    
47.       READ(5,*)((AE(I,J),I=1,4),J=1,NM)
    
48.       READ(5,*)(X(I),I=1,NP)
    
49.       READ(5,*)(Y(I),I=1,NP)
    
50.       READ(5,*)((MEO(I,J),I=1,2),J=1,NE)
    
51.       WRITE(6,81)((AE(I,J),I=1,4),J=1,NM)
    
52.       WRITE(6,82)(X(I),I=1,NP)
    
53.       WRITE(6,83)(Y(I),I=1,NP)
    
54.       WRITE(6,84)((MEO(I,J),I=1,2),J=1,NE)
    
55. 81   FORMAT(//25X,'EO*VO*W*T**'//(5X,4F15.4))
    
56. 82   FORMAT(//20X,'X--COORDINATE'//(5X,10F7.2))
    
57. 83   FORMAT(//20X,'Y--COORDINATE'//(5X,10F7.2))
    
58. 84   FORMAT(//20X,'ELEMENT-DATA'//(5X,10I7))
    
59. 90   CALL  MR
    
60.       CALL  FORMMA(AE,X,Y,MEO,NM,NP,NE,NV,A,SQ,AS)
    
61.       RETURN
    
62.       END
    
63. !     ****************************************************************
    
64.       SUBROUTINE MR
    
65.       DIMENSION JC(50),JR(2,800)
    
66.       COMMON/CA/NP,NE,NM,NR,NI,NL,NG,ND,NC
    
67.       COMMON/CC/N,NH,JR,R(1600)
    
68.       DO 10 I=1,NP
    
69.       DO 10 J=1,2
    
70.   10  JR(J,I)=1
    
71.       IF(NR) 20,20,30
    
72.   30  READ(5,*)(JC(I),I=1,NR)
    
73.       WRITE(6,90)(JC(I),I=1,NR)
    
74.   90  FORMAT(/20X,'CONSTRINED MESSAGE-JC='//(5X,10I7))
    
75.       DO 50 I=1,NR
    
76.       K=JC(I)
    
77.       J=K/100
    
78.       L=(K-J*100)/10
    
79.       M=K-J*100-L*10
    
80.       JR(1,J)=L
    
81.   50  JR(2,J)=M
    
82.   20  N=0
    
83.       DO 80 I=1,NP
    
84.       DO 80 J=1,2
    
85.       IF(JR(J,I)) 80, 80, 70
    
86.   70  N=N+1
    
87.       JR(J,I)=N
    
88.   80  CONTINUE
    
89.       RETURN
    
90.       END
    
91. !     **************************************************************
    
92.       SUBROUTINE FORMMA(AE,X,Y,MEO,KM,KP,KE,KV,A,SQ,AS)
    
93.       DIMENSION MA(2000),X(KP),Y(KP),MEO(2,KE),AE(4,KM),R(1600), &
    
94.                              JR(2,800),A(KV,KV),SQ(KV,KV),AS(KV,KV),ME(3),NN(6), &
    
95.                 SK(100000),RD(1600,20),R1(1600)
    
96.       
    
97.         COMMON/CA/NP,NE,NM,NR,NI,NL,NG,ND,NC,NA,NN1,DH,NV,NS,NQ
    
98.       COMMON/CC/N,NH,JR,R/CB/EO,VO,W,T,S,H11,H12,H21,H22,ME,BI(3), &
    
99.              CI(3),ER,TA1,TA2,NB,L/CD/EE(20),SS(20),A1(6),XA(20),        &
    
100.              YA(20),JA(20),SD(20),ED(20),SP(20),EP(20),A2(6,50),KL1, &
     
101.              KL2,CH(4)
     
102. !     DATA (XA(I),I=1,7)/88.,15.5,2.7,1.15E6,2.85,3.2E6,200.0/
     
103.       READ(5,*)(JA(I),I=1,NV)
     
104.       WRITE(6,7) (JA(I),I=1,NV)
     
105.    7  FORMAT (25X,'DISTRIBUTION  PARAMETER**********'//(10X,10I5))
     
106.       READ(5,*) NS
     
107.       WRITE(6,2) NS
     
108.   2   FORMAT(5X,'CORRELATED VARIABLE MESSAGE  **NS=',I5)
     
109.       IF(NS.EQ.0) GOTO 130
     
110.       READ(5,*) EOB
     
111.       WRITE(6,120) EOB
     
112. 120  FORMAT(10X,'CONTRAL ERROR    **EOB=',F13.5)
     
113.       CALL COV (NV,A,SQ,EOB)
     
114. 130  DO 8 I=1,NV
     
115.       XA(I)=0.0
     
116.       XA(I)=EE(I)
     
117.       YA(I)=0.0
     
118. 8    CONTINUE
     
119.       DO 10 I=1,N
     
120. 10   MA(I)=0
     
121.       DO 20 IE=1,NE
     
122.       CALL DIV(IE,AE,X,Y,MEO,NM,NP,NE)
     
123.       DO 30 I=1,3
     
124.       JB=ME(I)
     
125.       DO 30 M=1,2
     
126.       JJ=2*(I-1)+M
     
127.       NN(JJ)=JR(M,JB)
     
128. 30   CONTINUE
     
129.       L=N
     
130.       DO 80 I=1,6
     
131.       IF(NN(I)) 80,80,60
     
132. 60   IF(NN(I)-L) 70,80,80
     
133. 70   L=NN(I)
     
134. 80   CONTINUE
     
135.       DO 15 M=1,6
     
136.       JP=NN(M)
     
137.       IF(JP.EQ.0) GOTO 15
     
138.       IF(JP-L.GT.MA(JP))MA(JP)=JP-L
     
139. 15   CONTINUE
     
140. 20   CONTINUE
     
141.       MX=0
     
142.       MA(1)=1
     
143.       DO 95 I=2,N
     
144.       IF(MA(I).GT.MX) MX=MA(I)
     
145.       MA(I)=MA(I)+MA(I-1)+1
     
146.   95  CONTINUE
     
147.       MX=MX+1
     
148.       NH=MA(N)
     
149.       WRITE(6,115)N,NH,MX
     
150. 115  FORMAT(/22X,'T0TAL DEGREES OF FREEDOM N=',I5/25X,'TOTAL- &
     
151.             STORAGE','NH=',I6/24X,'MAX-SEMI-BANDWIDTH    MX=',I5)
     
152.       IF(NH.GT.100000) GOTO 111
     
153.       JDM=0
     
154.       DO 47 I=1,NV
     
155.       SP(I)=SS(I)
     
156.   47  EP(I)=EE(I)
     
157.       GOTO 333
     
158. 111  WRITE(6,222)
     
159. 222  FORMAT(/20X,'TOTAL STORAGE IS GREATER THAN 100000')
     
160.       STOP
     
161. 333  DO 5 I=1,NV
     
162.       SD(I)=SP(I)
     
163.       ED(I)=EP(I)
     
164.       J1=JA(I)
     
165.       JJ=J1-2
     
166.       IF (JJ) 5,3,4
     
167.    3  VV=(SS(I)/EE(I))**2
     
168.       SL=ALOG(1+VV)
     
169.       IF (XA(I).GE.0.0) GOTO 6
     
170.       XA(I)=ABS(XA(I))
     
171.    6  SP(I)=XA(I)*SQRT(SL)
     
172.       EP(I)=XA(I)*(1.0+ALOG(EE(I))-ALOG(XA(I)*SQRT(1.0+VV)))
     
173.       GOTO 5
     
174.   4   AR=1.28255/SS(I)
     
175.       AK=EE(I)-0.5772/AR
     
176.       Q=EXP(-AR*(XA(I)-AK))
     
177.       AF=EXP(-Q)
     
178.       AGF=AR*Q*EXP(-Q)
     
179.       AA=1.0-AF
     
180.       B0=1.570796288
     
181.       B1=3.706987906E-2
     
182.       B2=-8.364353589E-4
     
183.       B3=-2.250947176E-4
     
184.       B4=6.841218299E-6
     
185.       B5=5.824238515E-6
     
186.       B6=-1.04527497E-6
     
187.       B7=8.360937017E-8
     
188.       B8=-3.231081277E-9
     
189.       B9=3.657763036E-11
     
190.       B10=6.936233982E-13
     
191.       IF(AA-0.5) 230,240,250
     
192. 230   BB=AA
     
193.       GOTO 255
     
194. 250   BB=1.0-AA
     
195. 255   YY=-ALOG(4.0*BB*(1.0-BB))
     
196.       U1=YY*(B5+YY*(B6+YY*(B7+YY*(B8+YY*(B9+YY*B10)))))
     
197.       UU=YY*(B0+YY*(B1+YY*(B2+YY*(B3+YY*(B4+U1)))))
     
198.       UB=SQRT(UU)
     
199.       IF(AA.LT.0.5) GOTO 260
     
200.       UB=-UB
     
201.       GOTO 260
     
202. 240   UB=0.0
     
203. 260   SP(I)=EXP(-UB*UB/2.0)/SQRT(2.0*3.1416)/AGF
     
204.       EP(I)=XA(I)-SP(I)*UB
     
205.    5  CONTINUE
     
206.       DO 14 I=1,NM
     
207. !      K=2*(I+1)
     
208.           K=3*(I-1)+4
     
209.       AE(1,I)=XA(K)
     
210.   14  AE(3,I)=XA(K-1)
     
211.       JDM=JDM+1
     
212.       IF (JDM.GT.6) GOTO 75
     
213.       WRITE(6,23) JDM
     
214.   23  FORMAT(10X,'ITERATE NUMBER',5X,'JDM=',I5)
     
215.       WRITE(6,37)(SD(I),I=1,NV)
     
216.       WRITE(6,40)(ED(I),I=1,NV)
     
217.   37  FORMAT(25X,'THE LAST STANDARD DEVIATIONS ---SD'//(10X,5F12.3))
     
218.   40  FORMAT(25X,'THE LAST MEAN VALUES -----------ED'//(10X,5F12.3))
     
219.       WRITE(6,39)(SP(I),I=1,NV)
     
220.   39  FORMAT(25X,'THE NEXT STANDARD DEVIATIONS----SP'//(10X,5F12.3))
     
221.       WRITE(6,41)(EP(I),I=1,NV)
     
222.   41  FORMAT(25X,'THE NEXT MEAN VALUES------------EP'//(10X,5F12.3))
     
223.       WRITE(6,38)AE
     
224.   38  FORMAT(5X,'AE='//(4F13.4))
     
225.       CALL MGK(AE,X,Y,MEO,MA,SK,NM,NP,NE,N,NH)
     
226.       CALL LOAD(AE,X,Y,MEO,NM,NP,NE)
     
227.       WRITE(6,444)
     
228. 444   FORMAT(/35X,'NODAL FORCES'//15X,'NODE',13X,'X-COMP',13X,'Y-COMP')
     
229.       CALL OUTPUT
     
230.       IF(ND.GT.0) CALL TREAT(SK,MA,NH,N)
     
231.       CALL DECOMP(SK,MA,NH,N)
     
232.       CALL FOBA(SK,MA,NH,N)
     
233.          WRITE(6,555)
     
234. 555   FORMAT(/35X,'NODAL DISPLACEMENTS'//15X,'NODE',13X, &
     
235.          'X-COMP',13X,'Y-COMP')
     
236.          CALL OUTPUT
     
237.          WRITE(6,666)
     
238. 666    FORMAT(//35X,'ELEMENT STRESSES'//1X,'ELEMENT',2X,'X-STRESS',2X, &
     
239.                  'Y-STRESS',2X,'XY-STRESS',2X,'MAX-STRESS',2X,'MIN-STRESS' &
     
240.                  ,2X,'ANGLE'//)
     
241.       CALL CES(AE,X,Y,MEO,NM,NP,NE)
     
242.       IF(NC) 65,65,85
     
243.   85  CALL ERFAC(AE,X,Y,MEO,NM,NP,NE)
     
244.   65  CALL RIC(GY,AE,X,Y,MEO,MA,SK,NH,N,NP,NE,NM,RD,R1,NV,A,SQ,AS)
     
245.       GY=ABS(GY)
     
246.       IF (GY.GT.ER) GOTO 333
     
247.   75  RETURN
     
248.       END
     
249. !     ******************************************************************
     
250.       SUBROUTINE DIV(JJ,AE,X,Y,MEO,KM,KP,KE)
     
251.       DIMENSION ME(3),BI(3),CI(3),X(KP),Y(KP),AE(4,KM),MEO(2,KE)
     
252.       COMMON/CB/EO,VO,W,T,S,A(4),ME,BI,CI,ER,TA1,TA2,NB,L
     
253.       II=MEO(1,JJ)
     
254.       I=II/1000
     
255.       ME(1)=I
     
256.       J=II-I*1000
     
257.       ME(2)=J
     
258.       IJ=MEO(2,JJ)
     
259.       M=IJ/1000
     
260.       ME(3)=M
     
261.       IK=IJ-M*1000
     
262.       L=IK/10
     
263.       NB=IJ-M*1000-L*10
     
264.       BI(1)=Y(J)-Y(M)
     
265.       CI(1)=X(M)-X(J)
     
266.       BI(2)=Y(M)-Y(I)
     
267.       CI(2)=X(I)-X(M)
     
268.       BI(3)=Y(I)-Y(J)
     
269.       CI(3)=X(J)-X(I)
     
270. !     WRITE(6,2) I,J,M,L,NB,CI(2),BI(2),CI(3),BI(3)
     
271. !  2  FORMAT(/5X,5I5,5X,4F13.2)
     
272.       S=(BI(2)*CI(3)-CI(2)*BI(3))/2.0
     
273.       IF(S) 40,40,50
     
274.   40  WRITE(6,222) JJ
     
275. 222  FORMAT(/20X,'ZERO OR NEGATIVE AREA',10X,'ELEMENT NUMBER',I5)
     
276.       STOP
     
277.   50  EO=AE(1,L)
     
278.       VO=AE(2,L)
     
279.       W=AE(3,L)
     
280.       T=AE(4,L)
     
281.       RETURN
     
282.       END
     
283. !     *****************************************************************
     
284.       SUBROUTINE MGK(AE,X,Y,MEO,MA,SK,KM,KP,KE,KN,KH)
     
285.       DIMENSION AE(4,KM),X(KP),Y(KP),MEO(2,KE),MA(KN),SK(KH),        &
     
286.                 JR(2,800),NN(6),ME(3),BI(3),CI(3),B(6,6)
     
287.       COMMON /CA/NP,NE,NM,NR,NI,NL,NG,ND,NC,NA,NN1
     
288.       COMMON /CB/EO,VO,W,T,S,H11,H12,H21,H22,ME,BI,CI,ER,TA1,TA2,NB
     
289.       COMMON/CC/N,NH,JR,R(1600)/CD/EE(20),SS(20),A1(6), &
     
290.                 XA(20),YA(20),JA(20)
     
291.       DO 10 I=1,NH
     
292.   10  SK(I)=0.0
     
293.       DO 20 IE=1,NE
     
294.       CALL DIV(IE,AE,X,Y,MEO,NM,NP,NE)
     
295.       DO 30 I=1,3
     
296.       DO 30 J=1,3
     
297.       CALL KRS(BI(I),BI(J),CI(I),CI(J))
     
298.       B(2*I-1,2*J-1)=H11
     
299.       B(2*I-1,2*J)=H12
     
300.       B(2*I,2*J-1)=H21
     
301.       B(2*I,2*J)=H22
     
302.   30  CONTINUE
     
303.       DO 40 I=1,3
     
304.       J2=ME(I)
     
305.       DO 40 J=1,2
     
306.       J3=2*(I-1)+J
     
307.       NN(J3)=JR(J,J2)
     
308.   40  CONTINUE
     
309.       DO 60 I=1,6
     
310.       DO 60 J=1,6
     
311.       IF(NN(J).EQ.0.OR.NN(I).LT.NN(J)) GOTO 60
     
312.       JJ=NN(I)
     
313.       JK=NN(J)
     
314.       JL=MA(JJ)
     
315.       JM=JJ-JK
     
316.       JN=JL-JM
     
317.       SK(JN)=SK(JN)+B(I,J)
     
318.   60  CONTINUE
     
319.   20  CONTINUE
     
320.       WRITE(6,555)NN1
     
321. 555  FORMAT(5X,'NN1=',I3)
     
322.       IF(NN1.EQ.0) GOTO 100
     
323.       WRITE(6,444)(SK(I),I=1,NH)
     
324. 444  FORMAT(5X,'SK=',5F13.4)
     
325. 100  CONTINUE
     
326.       RETURN
     
327.       END
     
328. !     *****************************************************************
     
329.       SUBROUTINE LOAD(AE,X,Y,MEO,KM,KP,KE)
     
330.       DIMENSION AE(4,KM),X(KP),Y(KP),MEO(2,KE),ME(3),R(1600),NF(50), &
     
331.                 FV(2,50),JR(2,800)
     
332.       COMMON/CA/NP,NE,NM,NR,NI,NL,NG,ND,NC,NA,NN1,DH,NV,NS,NQ
     
333.       COMMON/CC/N,NH,JR,R/CB/EO,VO,W,T,S,A(4),ME,BI(3),CI(3),ER,TA1, &
     
334.                 TA2,NB,L
     
335.       COMMON/CD/EE(20),SS(20),AA(6),XA(20)
     
336.       DO 10 I=1,N
     
337.   10  R(I)=0.0
     
338.       IF(NG) 70,70,30
     
339.   30  DO 20 IE=1,NE
     
340.       CALL DIV(IE,AE,X,Y,MEO,NM,NP,NE)
     
341.       DO 50 I=1,3
     
342.       J2=ME(I)
     
343.       J3=JR(2,J2)
     
344.       IF(J3) 50,50,40
     
345.   40  R(J3)=R(J3)-T*W*S/3.0
     
346.   50  CONTINUE
     
347.   20  CONTINUE
     
348.   70  IF(NL) 200,200,80
     
349.   80  READ(5,*)(NF(I),I=1,NL)
     
350.       READ(5,*)((FV(I,J),I=1,2),J=1,NL)
     
351.       WRITE(6,25)(NF(I),I=1,NL)
     
352.       WRITE(6,35)((FV(I,J),I=1,2),J=1,NL)
     
353.   25  FORMAT(//20X,'NODES OF APPLIED LOAD**NF='//(10X,10I6))
     
354.   35  FORMAT(//20X,'LUMPED-LOADS**FV='//(10X,6F10.4))
     
355.       DO 100 I=1,NL
     
356.       JJ=NF(I)
     
357.       J=JR(1,JJ)
     
358.       M=JR(2,JJ)
     
359.       IF(J.GT.0)R(J)=R(J)+FV(1,I)
     
360.       IF(M.GT.0) R(M)=R(M)+FV(2,I)
     
361. 100  CONTINUE
     
362. 200  IF(NQ) 90,90,210
     
363. 210  DO 250 IE=1,NE
     
364.       CALL DIV(IE,AE,X,Y,MEO,NM,NP,NE)
     
365.       IF(L.GT.2) GOTO 250
     
366.       JK=ME(2)
     
367.       IF(Y(JK)-62.0) 215,230,220
     
368. 220  DO 218 I=1,3
     
369.       J2=ME(I)
     
370.       J3=JR(1,J2)
     
371.       IF(J3) 218,218,222
     
372. 222  R(J3)=R(J3)+0.643*T*S/3.0
     
373. 218  CONTINUE
     
374.       GOTO 250
     
375. 230  DO 228 I=1,3
     
376.       J2=ME(I)
     
377.       J3=JR(1,J2)
     
378.       IF(J3) 228,228,232
     
379. 232  R(J3)=R(J3)+0.421*T*S/3.0
     
380. 228  CONTINUE
     
381.       GOTO 250
     
382. 215  IF(Y(JK)-38.0) 240,260,270
     
383. 270  DO 275 I=1,3
     
384.       J2=ME(I)
     
385.       J3=JR(1,J2)
     
386.       IF(J3) 275,275,280
     
387. 280  R(J3)=R(J3)+0.272*T*S/3.0
     
388. 275  CONTINUE
     
389.       GOTO 250
     
390. 260  DO 264 I=1,3
     
391.       J2=ME(I)
     
392.       J3=JR(1,J2)
     
393.       IF(J3) 264,264,268
     
394. 268  R(J3)=R(J3)+0.223*T*S/3.0
     
395. 264  CONTINUE
     
396.       GOTO 250
     
397. 240  IF(Y(JK)-24.0) 290,310,320
     
398. 320  DO 325 I=1,3
     
399.       J2=ME(I)
     
400.       J3=JR(1,J2)
     
401.       IF(J3) 325,325,328
     
402. 328  R(J3)=R(J3)+0.186*T*S/3.0
     
403. 325  CONTINUE
     
404.       GOTO 250
     
405. 310  DO 315 I=1,3
     
406.       J2=ME(I)
     
407.       J3=JR(1,J2)
     
408.       IF(J3) 315,315,318
     
409. 318  R(J3)=R(J3)+0.148*T*S/3.0
     
410. 315  CONTINUE
     
411.       GOTO 250
     
412. 290  DO 330 I=1,3
     
413.       J2=ME(I)
     
414.       J3=JR(1,J2)
     
415.       IF(J3) 330,330,335
     
416. 335  R(J3)=R(J3)+0.124*T*S/3.0
     
417. 330  CONTINUE
     
418. 250  CONTINUE
     
419.   90  IF(NA) 120,120,130
     
420. 130  DO 150 IE=1,NE
     
421.       CALL DIV(IE,AE,X,Y,MEO,NM,NP,NE)
     
422.       IF(NB.GT.2) GOTO 160
     
423.       IF(NB.EQ.2) GOTO 140
     
424.       IK=ME(1)
     
425.       JK=ME(2)
     
426.       MK=JR(1,IK)
     
427.       NK=JR(1,JK)
     
428.       LK=JR(2,IK)
     
429.       KK=JR(2,JK)
     
430.       IF(XA(1).GE.Y(IK)) GOTO 132
     
431.       IF(XA(1).LT.Y(IK).AND.XA(1).GT.Y(JK)) GOTO 135
     
432.       GOTO 150
     
433. 135  C1=XA(1)-Y(JK)
     
434.       C2=Y(IK)-Y(JK)
     
435.       R(MK)=R(MK)+C1*C1*C1/(6.0*C2)
     
436.       R(NK)=R(NK)+(3.0*C1*C1*C2-C1*C1*C1)/(6.0*C2)
     
437.       GOTO 150
     
438. 132  IF (Y(IK).EQ.Y(JK)) GOTO 133
     
439.       C3=XA(1)-Y(IK)
     
440.       C4=XA(1)-Y(JK)
     
441.       C5=Y(IK)-Y(JK)
     
442.       R(MK)=R(MK)+(C4*C5)/6.0+(C3*C5)/3.0
     
443.       R(NK)=R(NK)+(C3*C5)/6.0+(C4*C5)/3.0
     
444.       GOTO 150
     
445. 133  C6=XA(1)-Y(IK)
     
446.       C7=X(IK)-X(JK)
     
447.       IF(LK.EQ.0) GOTO 134
     
448.       R(LK)=R(LK)-0.5*C6*C7
     
449. 134  IF(KK.EQ.0) GOTO 150
     
450.       R(KK)=R(KK)-0.5*C6*C7
     
451.       GOTO 150
     
452. 140  K1=ME(1)
     
453.       K2=ME(2)
     
454.       K3=JR(1,K1)
     
455.       K4=JR(2,K1)
     
456.       K5=JR(1,K2)
     
457.       K6=JR(2,K2)
     
458.       IF(XA(2).GE.Y(K2)) GOTO 145
     
459.       IF(XA(2).GT.Y(K1).AND.XA(2).LT.Y(K2)) GOTO 147
     
460.       GOTO 150
     
461. 147  A1=XA(2)-Y(K1)
     
462.       A2=Y(K2)-Y(K1)
     
463.       A3=3*A1*A1*A2-A1*A1*A1
     
464.       R(K3)=R(K3)-A3/(6.0*A2)
     
465.       R(K4)=R(K4)-A3*TA2/(6.0*A2)
     
466.       R(K5)=R(K5)-A1*A1*A1/(6.0*A2)
     
467.       R(K6)=R(K6)-A1*A1*A1*TA2/(6.0*A2)
     
468.       GOTO 150
     
469. 145  IF (Y(K1).EQ.Y(K2)) GOTO 146
     
470.       B1=XA(2)-Y(K1)
     
471.       B2=XA(2)-Y(K2)
     
472.       B3=Y(K2)-Y(K1)
     
473.       R(K3)=R(K3)-(B2/6.0+B1/3.0)*B3
     
474.       R(K4)=R(K4)-(B2/6.0+B1/3.0)*B3*TA2
     
475.       R(K5)=R(K5)-(B1/6.0+B2/3.0)*B3
     
476.       R(K6)=R(K6)-(B1/6.0+B2/3.0)*B3*TA2
     
477.       GOTO 150
     
478. 146  E1=XA(2)-Y(K1)
     
479.       E2=X(K1)-X(K2)
     
480.       IF(K4.EQ.0) GOTO 148
     
481.       R(K4)=R(K4)-0.5*E1*E2
     
482. 148  IF(K6.EQ.0) GOTO 150
     
483.       R(K6)=R(K6)-0.5*E1*E2
     
484.       GOTO 150
     
485. 160  IF(NB.EQ.3) GOTO 150
     
486. !      IK=ME(1)
     
487. !      JK=ME(2)
     
488. !      MK=JR(1,IK)
     
489. !      NK=JR(1,JK)
     
490. !      LK=JR(2,IK)
     
491. !      KK=JR(2,JK)
     
492. !      IF(XA(1).GE.Y(IK)) GOTO 170
     
493. !      IF(XA(1).LT.Y(IK).AND.XA(1).GT.Y(JK)) GOTO 180
     
494. !      GOTO 150
     
495. ! 170  D1=XA(1)-Y(IK)
     
496. !      D2=Y(IK)-Y(JK)
     
497. !      D3=XA(1)-Y(JK)
     
498. !      R(MK)=R(MK)+(D1/3.0+D3/6.0)*D2
     
499. !      R(LK)=R(LK)-(D1/3.0+D3/6.0)*D2*TA1
     
500. !      R(NK)=R(NK)+(D1/6.0+D3/3.0)*D2
     
501. !     R(KK)=R(KK)-(D1/6.0+D3/3.0)*D2*TA1
     
502. !      GOTO 150
     
503. ! 180  D4=XA(1)-Y(JK)
     
504. !      D5=Y(IK)-Y(JK)
     
505. !      D6=D4*D4*D4
     
506. !      D7=D4*D4*D5
     
507. !      R(MK)=R(MK)+D6/6.0/D5
     
508. !      R(LK)=R(LK)-D6*TA1/6.0/D5
     
509. !          R(NK)=R(NK)+(3.0*D7-D6)/6.0/D5
     
510. !      R(KK)=R(KK)-(3.0*D7-D6)*TA1/6.0/D5
     
511. 150  CONTINUE
     
512. 120  RETURN
     
513.       END
     
514. !     ******************************************************************
     
515.       SUBROUTINE OUTPUT
     
516.       DIMENSION JR(2,800),R(1600)
     
517.       COMMON/CA/NP,NE,NM,NR,NI,NL,NG,ND,NC/CC/N,NH,JR,R
     
518.       DO 100 I=1,NP
     
519.       L=JR(1,I)
     
520.       IF(L) 10,20,30
     
521.   30  S=R(L)
     
522.       GOTO 10
     
523.   20  S=0.0
     
524.   10  L=JR(2,I)
     
525.       IF(L) 40,50,60
     
526.   60  SS=R(L)
     
527.       GOTO 40
     
528.   50  SS=0.0
     
529.   40  WRITE(6,75) I,S,SS
     
530.   75  FORMAT(8X,I4,5X,F20.8,5X,F20.8)
     
531. 100  CONTINUE
     
532.       RETURN
     
533.       END
     
534. !     *****************************************************************
     
535.       SUBROUTINE DECOMP(SK,MA,KH,KN)
     
536.       DIMENSION SK(KH),MA(KN)
     
537.       COMMON/CC/N,NH,JR(2,800),R(1600)
     
538.       DO 130 I=2,N
     
539.       L=MA(I-1)+I-MA(I)+1
     
540.       K=I-1
     
541.       DO 280 J=L,K
     
542.       IF(J-L) 20,280,20
     
543.   20  JP=MA(I)-I+J
     
544.       M=MA(J-1)+J-MA(J)+1
     
545.       IF(L.GT.M) M=L
     
546.       MP=J-1
     
547.       DO 230 IP=M,MP
     
548.       JJ=MA(I)-I+IP
     
549.       JK=MA(J)-J+IP
     
550.       SK(JP)=SK(JP)-SK(JJ)*SK(JK)
     
551. 230  CONTINUE
     
552. 280  CONTINUE
     
553.       DO 400 IP=L,K
     
554.       JI=MA(I)-I+IP
     
555.       JL=MA(IP)
     
556.       SK(JI)=SK(JI)/SK(JL)
     
557.       JN=MA(I)
     
558.       SK(JN)=SK(JN)-SK(JI)*SK(JI)*SK(JL)
     
559. 400  CONTINUE
     
560. 130  CONTINUE
     
561.       RETURN
     
562.       END
     
563. !     *****************************************************************
     
564.       SUBROUTINE FOBA(SK,MA,KH,KN)
     
565.       DIMENSION SK(KH),MA(KN),JR(2,800),R(1600)
     
566.       COMMON/CC/N,NH,JR,R/CA/NP,NE,NM,NR,NI,NL,NG,ND,NC
     
567.       DO 210 I=2,N
     
568.       JJ=MA(I)
     
569.       JK=I-1
     
570.       JL=MA(JK)+I-JJ+1
     
571.       DO 210 J=JL,JK
     
572.       JP=MA(I)-I+J
     
573.       R(I)=R(I)-SK(JP)*R(J)
     
574. 210  CONTINUE
     
575.       DO 220 I=1,N
     
576.       JJ=MA(I)
     
577.       R(I)=R(I)/SK(JJ)
     
578. 220  CONTINUE
     
579.       DO 230 J4=2,N
     
580.       I=2+N-J4
     
581.       JJ=MA(I-1)+I-MA(I)+1
     
582.       M=MA(I)-I
     
583.       JP=I-1
     
584.       DO 240 J=JJ,JP
     
585.       JM=J+M
     
586.       R(J)=R(J)-SK(JM)*R(I)
     
587. 240  CONTINUE
     
588. 230  CONTINUE
     
589.       RETURN
     
590.       END
     
591. !     ******************************************************************
     
592.       SUBROUTINE CES(AE,X,Y,MEO,KM,KP,KE)
     
593.       DIMENSION AE(4,KM),X(KP),Y(KP),MEO(2,KE),R(1600),JR(2,800), &
     
594.                 ME(3),BI(3),CI(3),B(6),CC(3)
     
595.       COMMON/CA/NP,NE,NM,NR,NI,NL,NG,ND,NC
     
596.       COMMON/CC/N,NH,JR,R/CB/EO,VO,W,T,S,H11,H12,H21,H22,ME,BI,CI, &
     
597.                 ER,TA1,TA2,NB,L
     
598.       COMMON/CD/EE(20),SS(20),A1(6),XA(20),YA(20),JA(20),SD(20), &
     
599.                 ED(20),SP(20),EP(20),A2(6,50),KL1,KL2,CH(4),KL
     
600.       DO 6 I=1,6
     
601.       A1(I)=0.0
     
602.    6  CONTINUE
     
603.       DO 10 IE=1,NE
     
604.       CALL DIV (IE,AE,X,Y,MEO,NM,NP,NE)
     
605.       V1=VO/(1.0-VO)
     
606.       V2=EO/(1.0-VO*VO)
     
607.       ET=V2/(1.0-V1*V1)/S/2.0
     
608.       DO 55 I=1,3
     
609.       J2=ME(I)
     
610.       I2=JR(1,J2)
     
611.       I3=JR(2,J2)
     
612.       IF(I2) 50,60,70
     
613.   70  B(2*I-1)=R(I2)
     
614.       GOTO 50
     
615.   60  B(2*I-1)=0.0
     
616.   50  IF(I3) 55,65,75
     
617.   75  B(2*I)=R(I3)
     
618.       GOTO 55
     
619.   65  B(2*I)=0.0
     
620.   55  CONTINUE
     
621.       H1=0.0
     
622.       H2=0.0
     
623.       H3=0.0
     
624.       DO 100 I=1,3
     
625.       H1=H1+BI(I)*B(2*I-1)
     
626.       H2=H2+CI(I)*B(2*I)
     
627.       H3=H3+BI(I)*B(2*I)+CI(I)*B(2*I-1)
     
628. 100  CONTINUE
     
629.       AA1=ET*(H1+V1*H2)
     
630.       A4=ET*(H2+V1*H1)
     
631.       A3=ET*(1.0-V1)*H3/2.0
     
632.       H1=AA1+A4
     
633.       H2=SQRT((AA1-A4)*(AA1-A4)+4.0*A3*A3)
     
634.       B(4)=(H1+H2)/2.0
     
635.       B(5)=(H1-H2)/2.0
     
636.       IF(ABS(A3).GT.1E-4) GOTO 400
     
637.       B(6)=90.0
     
638.       GOTO 450
     
639. 400  B(6)=ATAN((B(4)-AA1)/A3)*57.29578
     
640. 450  B(1)=AA1
     
641.       B(2)=A4
     
642.       B(3)=A3
     
643.           IF (IE.EQ.2)GOTO 1001
     
644.           GOTO 1002       
     
645. 1001  KL=IE
     
646.           DO 1000 I=1,6
     
647.                    A1(I)=B(I)
     
648. 1000  CONTINUE
     
649. !     IF(IE.GE.KL1.AND.IE.LE.KL2) GOTO 220
     
650. !          IF(IE.GE.KL1) GOTO 220
     
651. !          GOTO 445
     
652. ! 220  IQ=IE-KL1+1
     
653. !      DO 230 I=1,6
     
654. !      A2(I,IQ)=B(I)
     
655. !          230  CONTINUE
     
656. 1002      WRITE(6,555)IE,(B(I),I=1,6)
     
657. 555   FORMAT(4X,I4,2X,F8.3,2X,F8.3,2X,F8.3,3X,F10.4,2X,F10.4,2X,F8.3)
     
658. ! 445  CONTINUE
     
659.   10  CONTINUE
     
660.        WRITE(6,999)KL,(A1(I),I=1,6)
     
661. 999   FORMAT(4X,'KL=',I4,2X,F8.3,2X,F8.3,2X,F8.3,3X,F10.4,2X,F10.4, &
     
662.               2X,F8.3)
     
663. 
     
664. !     IF (NI.LT.0) GOTO 411
     
665. !     C1=(CH(4)-CH(2))*(CH(4)-CH(3))/(CH(1)-CH(2))/(CH(1)-CH(3))
     
666. !     C2=(CH(4)-CH(1))*(CH(4)-CH(3))/(CH(2)-CH(1))/(CH(2)-CH(3))
     
667. !     C3=(CH(4)-CH(1))*(CH(4)-CH(2))/(CH(3)-CH(1))/(CH(3)-CH(2))
     
668. !     DO 150 I=1,3
     
669. !     DO 160 J=KL1,KL2,2
     
670. !     P=J-KL1+1
     
671. !           P=J-1+1
     
672. !      L=P+1
     
673. !     M=(J-KL1)/2+1
     
674. !      CC(M)=(A2(I,P)+A2(I,L))/2.0
     
675. !160  CONTINUE
     
676. !     A1(I)=C1*CC(1)+C2*CC(2)+C3*CC(3)
     
677. !150  CONTINUE
     
678. !     C4=(A1(1)+A1(2))/2.0
     
679. !     C5=(A1(1)-A1(2))/2.0
     
680. !     A1(4)=C4+SQRT(C5*C5+A1(3)*A1(3))
     
681. !     A1(5)=C4-SQRT(C5*C5+A1(3)*A1(3))
     
682. !     WRITE(6,610) A1
     
683. !610  FORMAT(5X,'A1(I)**=',6(2X,F8.3))
     
684. !411  IP=KL2-KL1+1
     
685. !     WRITE(6,421) ((A2(I,J),I=1,6),J=1,IP)
     
686. ! 421  FORMAT(2X,'A2(KL1---KL2 STRESSES)'/(5X,6F12.4))
     
687.       RETURN
     
688.       END
     
689. !     *****************************************************************
     
690.       SUBROUTINE ERFAC (AE,X,Y,MEO,KM,KP,KE)
     
691.       DIMENSION NCI(20),NCE(4,20),ME(3),BI(3),R(1600),JR(2,800), &
     
692.                  AE(4,KM),X(KP),Y(KP),MEO(2,KE),CI(3)
     
693.       COMMON /CA/NP,NE,NM,NR,NI,NL,NG,ND,NC,NA,NN1
     
694.       COMMON /CC/N,NH,JR,R
     
695.       COMMON/CB/EO,VO,W,T,S,H11,H12,H21,H22,ME,BI,CI
     
696.       READ(5,*) (NCI(J),J=1,NC)
     
697.       READ (5,*) ((NCE(I,J),I=1,4),J=1,NC)
     
698.       WRITE(6,35)(NCI(J),J=1,NC)
     
699.       WRITE(6,45)((NCE(I,J),I=1,4),J=1,NC)
     
700.   35  FORMAT(//30X,'NODN-NAME***NCI='//(10X,10I6))
     
701.   45  FORMAT(//30X,'ELEMENT-NAME***NCE='//(10X,12I5))
     
702.       WRITE(6,999)
     
703. 999  FORMAT(//30X,'NODAL REACTIONS'//15X,'NODE',13X'X-COMP', &
     
704.              13X'Y-COMP')
     
705.       DO 20 JJ=1,NC
     
706.       FX=0.0
     
707.       FY=0.0
     
708.       L=NCI(JJ)
     
709.       DO 80 M=1,4
     
710.       IF(NCE(M,JJ)) 80,80,70
     
711.   70  IE=NCE(M,JJ)
     
712.       CALL DIV(IE,AE,X,Y,MEO,NM,NP,NE)
     
713.       DO 100 IM=1,3
     
714.       K=IM
     
715.       IF(L-ME(IM)) 100,110,100
     
716. 100  CONTINUE
     
717.       WRITE(6,555) L
     
718. 555  FORMAT(//20X,'ERROR OF ELEMENT MESSAGE******NODE NUMBER',I5)
     
719. 110  DO 400 IP=1,3
     
720.       CALL KRS(BI(K),BI(IP),CI(K),CI(IP))
     
721.       NL=ME(IP)
     
722.       JI=JR(1,NL)
     
723.       JP=JR(2,NL)
     
724.       IF(JI) 210,210,220
     
725. 210  S=0.0
     
726.       GOTO 200
     
727. 220  S=R(JI)
     
728. 200  IF(JP) 310,310,320
     
729. 310  SS=0.0
     
730.       GOTO 300
     
731. 320  SS=R(JP)
     
732. 300  FX=FX+H11*S+H12*SS
     
733.       FY=FY+H21*S+H22*SS
     
734. 400  CONTINUE
     
735.   80  CONTINUE
     
736.       WRITE(6,444) L,FX,FY
     
737. 444  FORMAT(14X,I6,10X,F12.4,10X,F12.4)
     
738.   20  CONTINUE
     
739.       RETURN
     
740.       END
     
741. !     ****************************************************************
     
742.       SUBROUTINE KRS(BR,BS,CR,CS)
     
743.       COMMON/CB/EO,VO,W,T,S,H11,H12,H21,H22,ME(3),BI(3),CI(3)
     
744.       ET=EO*T*(1.0-VO)/4.0/S/(1.0+VO)/(1.0-2.0*VO)
     
745.       V1=VO/(1.0-VO)
     
746.       V2=(1.0-2.0*VO)/2.0/(1.0-VO)
     
747.       H11=ET*(BR*BS+V2*CR*CS)
     
748.       H12=ET*(V1*BR*CS+V2*BS*CR)
     
749.       H21=ET*(V1*CR*BS+V2*BR*CS)
     
750.       H22=ET*(CR*CS+V2*BR*BS)
     
751.       RETURN
     
752.       END
     
753. !     ****************************************************************
     
754.       SUBROUTINE TREAT(SK,MA,KH,KN)
     
755.       DIMENSION SK(KH),MA(KN),R(1600),NDI(10),DV(2,10),JR(2,800)
     
756.       COMMON/CA/NP,NE,NM,NR,NI,NL,NG,ND,NC
     
757.       COMMON/CC/N,NH,JR,R
     
758.       READ(5,*)(NDI(J),J=1,ND)
     
759.       READ(5,*)((DV(I,J),I=1,2),J=1,ND)
     
760.       WRITE(6,35)(NDI(J),J=1,ND)
     
761.       WRITE(6,45)((DV(I,J),I=1,2),J=1,ND)
     
762.   35  FORMAT(//25X,'NODE-NAME**NDI='//(10X,10I6))
     
763.   45  FORMAT(//20X,'DISPLACEMENT-VALUES**DV'//(10X,6F10.4))
     
764.       DO 120 I=1,ND
     
765.       DO 120 J=1,2
     
766.       IF(DV(I,J)) 70,120,70
     
767.   70  JJ=NDI(I)
     
768.       L=JR(J,JJ)
     
769.       IF(L.EQ.0) GOTO 120
     
770.       JN=MA(L)
     
771.       SK(JN)=1E15
     
772.       R(L)=DV(I,J)*1E15
     
773. 120  CONTINUE
     
774.       RETURN
     
775.       END
     
776. !     **************************************************************
     
777.       SUBROUTINE RIC(GY,AE,X,Y,MEO,MA,SK,KH,KN,KP,KE,KM,RD,R1,KV, &
     
778.                      A,SQ,AS)
     
779.       DIMENSION FM(6,20),RD(KN,KV),FK(6),R1(KN), &
     
780.                 YA(20),GA(20),A(KV,KV),SQ(KV,KV),AS(KV,KV),        &
     
781.                 XA(20),AE(4,KM),X(KP),Y(KP),MEO(2,KE),BI(3), &
     
782.                 CI(3),B(6,6),NN(6),SK(KH),MA(KN)
     
783.       COMMON/CA/NP,NE,NM,NR,NI,NL,NG,ND,NC,NA,NN1,DH,NV,NS
     
784.       COMMON/CC/N,NH,JR(2,800),R(1600)
     
785.       COMMON/CB/EO,VO,W,T,S,H11,H12,H21,H22,ME(3),BI,CI,ER,TA1, &
     
786.                 TA2,NB,L
     
787.       COMMON/CD/EE(20),SS(20),A1(6),XA,YA,JA(20),SD(20), &
     
788.                 ED(20),SP(20),EP(20),A2(6,50),KL1,KL2,CH(4),KL
     
789.       DO 3 I=1,N
     
790.       DO 3 J=1,NV
     
791.       RD(I,J)=0.0
     
792.   3   CONTINUE
     
793.       DO 10 IE=1,NE
     
794.       DO 5 I=1,6
     
795.       DO 5 J=1,NV
     
796.       FM(I,J)=0.0
     
797.   5   CONTINUE
     
798.       CALL DIV(IE,AE,X,Y,MEO,NM,NP,NE)
     
799.       LL=2*L+1
     
800.       FM(2,LL)=-S*T/3.0
     
801.       FM(4,LL)=FM(2,LL)
     
802.       FM(6,LL)=FM(2,LL)
     
803.       J1=ME(1)
     
804.       J2=ME(2)
     
805.       Y1=Y(J1)
     
806.       Y2=Y(J2)
     
807.       X1=X(J1)
     
808.       X2=X(J2)
     
809.       NQ=NB-2
     
810.       IF (NQ) 11,12,13
     
811.   11  IF (XA(1).GT.Y1) GOTO 15
     
812.       IF (XA(1).EQ.Y1) GOTO 38
     
813.       IF (XA(1).LT.Y1.AND.XA(1).GT.Y2) GOTO 20
     
814. !      GOTO 40
     
815.           GOTO 13
     
816.   15  IF(Y1.EQ.Y2) GOTO 35
     
817. !      FM(1,1)=(Y1-Y2)/2.0
     
818. !      FM(3,1)=(Y1-Y2)/2.0
     
819.           FM(1,1)=(Y1-Y2)/2.0
     
820.       FM(3,1)=(Y1-Y2)/2.0
     
821.           FM(2,3)=-S*T/3.0
     
822.           FM(4,3)=FM(2,3)
     
823.           FM(6,3)=FM(2,3)
     
824.       GOTO 40
     
825. !  35  FM(2,1)=(X2-X1)*0.5
     
826. !      FM(4,1)=FM(2,1)
     
827.    35  FM(2,1)=(X2-X1)*0.5
     
828.        FM(4,1)=FM(2,1)
     
829.            FM(2,6)=-S*T/3.0
     
830.            FM(4,6)=FM(2,6)
     
831.            FM(6,6)=FM(2,6)
     
832.       GOTO 40
     
833. !  38  FM(1,1)=(Y1-Y2)/6.0
     
834. !      FM(3,1)=(Y1-Y2)/3.0
     
835.   38  FM(1,1)=(Y1-Y2)/6.0
     
836.       FM(3,1)=(Y1-Y2)/3.0
     
837.           FM(2,3)=-S*T/3.0
     
838.           FM(4,3)=FM(2,3)
     
839.           FM(6,3)=FM(2,3)
     
840.       GOTO 40
     
841. !  20  FM(1,1)=3.*(XA(1)-Y2)**2/6./(Y1-Y2)
     
842. !      FM(3,1)=(6.*(XA(1)-Y2)*(Y1-Y2)-3.*(XA(1)-Y2)**2)/6./(Y1-Y2)
     
843.   20  FM(1,1)=3.*(XA(1)-Y2)**2/6./(Y1-Y2)
     
844.       FM(3,1)=(6.*(XA(1)-Y2)*(Y1-Y2)-3.*(XA(1)-Y2)**2)/6./(Y1-Y2)
     
845.           FM(2,3)=-S*T/3.0
     
846.           FM(4,3)=FM(2,3)
     
847.           FM(6,3)=FM(2,3)
     
848.       GOTO 40
     
849.   12  IF (XA(2).GT.Y2) GOTO 22
     
850.       IF (XA(2).EQ.Y2) GOTO 26
     
851.       IF (XA(2).LT.Y2.AND.XA(2).GT.Y1) GOTO 24
     
852. !      GOTO 40
     
853.           GOTO 13
     
854.   22  IF (Y1.EQ.Y2) GOTO 23
     
855. !      FM(1,2)=-(Y2-Y1)/2.0
     
856. !      FM(2,2)=FM(1,2)*TA2
     
857. !      FM(3,2)=FM(1,2)
     
858. !      FM(4,2)=FM(2,2)
     
859.           FM(1,2)=-(Y2-Y1)/2.0
     
860.       FM(2,2)=FM(1,2)*TA2
     
861.       FM(3,2)=FM(1,2)
     
862.       FM(4,2)=FM(2,2)
     
863.           FM(2,3)=-S*T/3.0
     
864.           FM(4,3)=FM(2,3)
     
865.           FM(6,3)=FM(2,3)
     
866.       GOTO 40
     
867. !  23  FM(2,2)=0.5*(X2-X1)
     
868. !      FM(4,2)=FM(2,2)
     
869.   23  FM(2,2)=0.5*(X2-X1)
     
870.       FM(4,2)=FM(2,2)
     
871.           FM(2,6)=-S*T/3.0
     
872.           FM(4,6)=FM(2,6)
     
873.           FM(6,6)=FM(2,6)
     
874.       GOTO 40
     
875. !  26  FM(1,2)=-(Y2-Y1)/3.0
     
876. !      FM(2,2)=FM(1,2)*TA2
     
877. !      FM(3,2)=-(Y2-Y1)/6.0
     
878. !      FM(4,2)=FM(3,2)*TA2
     
879.   26  FM(1,2)=-(Y2-Y1)/3.0
     
880.       FM(2,2)=FM(1,2)*TA2
     
881.       FM(3,2)=-(Y2-Y1)/6.0
     
882.       FM(4,2)=FM(3,2)*TA2
     
883.           FM(2,3)=-S*T/3.0
     
884.           FM(4,3)=FM(2,3)
     
885.           FM(6,3)=FM(2,3)
     
886.       GOTO 40
     
887. !  24  FM(1,2)=-(6.*(XA(2)-Y1)*(Y2-Y1)-3.*(XA(2)-Y1)**2)/6./(Y2-Y1)
     
888. !      FM(2,2)=-(6.*(XA(2)-Y1)*(Y2-Y1)-3.*(XA(2)-Y1)**2)*TA2/6./(Y2-Y1)
     
889. !      FM(3,2)=-3.*(XA(2)-Y1)**2/6./(Y2-Y1)
     
890. !      FM(4,2)=-3.*(XA(2)-Y1)**2*TA2/6./(Y2-Y1)
     
891.   24  FM(1,2)=-(6.*(XA(2)-Y1)*(Y2-Y1)-3.*(XA(2)-Y1)**2)/6./(Y2-Y1)
     
892.       FM(2,2)=-(6.*(XA(2)-Y1)*(Y2-Y1)-3.*(XA(2)-Y1)**2)*TA2/6./(Y2-Y1)
     
893.       FM(3,2)=-3.*(XA(2)-Y1)**2/6./(Y2-Y1)
     
894.       FM(4,2)=-3.*(XA(2)-Y1)**2*TA2/6./(Y2-Y1)
     
895.           FM(2,3)=-S*T/3.0
     
896.           FM(4,3)=FM(2,3)
     
897.           FM(6,3)=FM(2,3)
     
898.       GOTO 40
     
899.   13  FM(2,3)=-S*T/3.0
     
900.           FM(4,3)=FM(2,3)
     
901.           FM(6,3)=FM(2,3)
     
902.           IF(NB.EQ.3) GOTO 40
     
903. !      IF(XA(1).GT.Y1) GOTO 19
     
904. !      IF(XA(1).EQ.Y1) GOTO 21
     
905. !      IF(XA(1).LT.Y1.AND.XA(1).GT.Y2) GOTO 33
     
906. !      GOTO 40
     
907. !  19  FM(1,1)=0.5*(Y1-Y2)
     
908. !      FM(2,1)=-FM(1,1)*TA1
     
909. !          FM(2,1)=-FM(1,1)
     
910. !      FM(3,1)=FM(1,1)
     
911. !      FM(4,1)=FM(2,1)
     
912. !      GOTO 40
     
913. !  21  FM(1,1)=(Y1-Y2)/6.
     
914. !      FM(2,1)=-FM(1,1)*TA1
     
915. !          FM(2,1)=-FM(1,1)
     
916. !      FM(3,1)=FM(1,1)*2.
     
917. !      FM(4,1)=FM(2,1)*2.
     
918. !      GOTO 40
     
919. !  33  D8=XA(1)-Y2
     
920. !      D9=Y1-Y2
     
921. !     FM(1,1)=D8*D8/D9/2.
     
922. !      FM(2,1)=-FM(1,1)*TA1
     
923. !      FM(3,1)=(6.*D9*D8-3.*D8*D8)/D9/6.
     
924. !      FM(4,1)=-FM(3,1)*TA1
     
925.   40  CONTINUE
     
926. !      WRITE(6,16)((FM(I,J),J=1,NV),I=1,6)
     
927. !  16  FORMAT(25X,'FM ********=='//(10X,5F12.4))
     
928.       DO 45 K=1,NV
     
929.       DO 45 I=1,3
     
930.       J2=ME(I)
     
931.       DO 45 J3=1,2
     
932.       K1=2*(I-1)+J3
     
933.       J4=JR(J3,J2)
     
934.       IF (J4.GT.0) RD(J4,K)=RD(J4,K)+FM(K1,K)
     
935.   45  CONTINUE
     
936. !      LL=2*(L+1)
     
937.           LL=3*(L-1)+4
     
938.       DO 50 I=1,3
     
939.       DO 50 J=1,3
     
940.       CALL KRS(BI(I),BI(J),CI(I),CI(J))
     
941.       B(2*I-1,2*J-1)=H11/XA(LL)
     
942.       B(2*I-1,2*J)=H12/XA(LL)
     
943.       B(2*I,2*J-1)=H21/XA(LL)
     
944.       B(2*I,2*J)=H22/XA(LL)
     
945.   50  CONTINUE
     
946.       DO 52 I=1,6
     
947.       FK(I)=0.0
     
948.   52  CONTINUE
     
949.       DO 55 I=1,3
     
950.       J2=ME(I)
     
951.       DO 55 J=1,2
     
952.       J3=2*(I-1)+J
     
953.       NN(J3)=JR(J,J2)
     
954.   55  CONTINUE
     
955.       DO 60 I=1,6
     
956.       DO 60 J=1,6
     
957.       NJ=NN(J)
     
958.       IF (NJ.EQ.0) GOTO 60
     
959.       FK(I)=FK(I)+B(I,J)*R(NJ)
     
960.   60  CONTINUE
     
961. !      WRITE(6,61) (FK(I),I=1,6)
     
962. !  61  FORMAT (25X,'FK**='//5X,3E16.6)
     
963.       DO 65 I=1,3
     
964.       J2=ME(I)
     
965.       DO 65 J=1,2
     
966.       K1=2*(I-1)+J
     
967.       NJ=JR(J,J2)
     
968.       IF (NJ.GT.0) RD(NJ,LL)=RD(NJ,LL)-FK(K1)
     
969.   65  CONTINUE
     
970.   10  CONTINUE
     
971. !      WRITE(6,64) ((RD(I,J),J=1,NV),I=1,N)
     
972. !  64  FORMAT(25X,'RD********='//(5X,5E14.6))
     
973.       DO 18 I=1,N
     
974.       R1(I)=R(I)
     
975.   18  CONTINUE
     
976.       DO 70 I=1,NV
     
977.       DO 75 J=1,N
     
978.       R(J)=RD(J,I)
     
979.   75  CONTINUE
     
980.       CALL FOBA(SK,MA,NH,N)
     
981.       DO 80 J=1,N
     
982.       RD(J,I)=R(J)
     
983.   80  CONTINUE
     
984.   70  CONTINUE
     
985. !      WRITE(6,71)((RD(I,J),J=1,NV),I=1,N)
     
986. !  71  FORMAT(25X,'RD*******J='//(5X,5E14.6))
     
987.       IF(NI.LT.0) GOTO 444
     
988.       CALL RI1(GY,AE,X,Y,MEO,RD,GA,A,SQ,AS,N,NP,NE,NM,NV)
     
989.       GOTO 555
     
990. 444  CALL RI2(GY,AE,X,Y,MEO,RD,GA,A,SQ,AS,N,NP,NE,NM,NV)
     
991. 555  BB=0.0
     
992.       DO 175 I=1,NV
     
993.       BB=BB+YA(I)*YA(I)
     
994. 175  CONTINUE
     
995.       BB=SQRT(BB)
     
996.       WRITE(6,177) BB
     
997. 177  FORMAT(25X,'RELIABILITY-INDEX'/25X,'*****************'/25X, &
     
998.               'BB=',F16.8)
     
999.       WRITE(6,185) (XA(I),I=1,NV)
     
1000. 185  FORMAT(10X,'XXXXX'//(10X,5F13.4))
      
1001. !      IF (XA(1).LE.DH) GOTO 182
      
1002. !      XA(1)=DH
      
1003. !      WRITE(6,181) XA(1)
      
1004. ! 181  FORMAT(10X,'XA(1)=DH***',F15.5)
      
1005.       WRITE(6,190) (YA(I),I=1,NV)
      
1006. 190  FORMAT(10X,'YYYYY'//(10X,5F13.4))
      
1007.       WRITE(6,192) (GA(I),I=1,NV)
      
1008. 192  FORMAT(10X,'GGGAAA'//(10X,5F13.4))
      
1009.       WRITE(6,195) GY
      
1010. 195  FORMAT(20X,'GY=',F13.6)
      
1011.       RETURN
      
1012.       END
      
1013. 
      
1014. !  ****************************************************************
      
1015.       SUBROUTINE RI1(GY,AE,X,Y,MEO,RD,GA,A,SQ,AS,KN,KP,KE,KM,KV)
      
1016.       DIMENSION AE(4,KM),X(KP),Y(KP),MEO(2,KE),RD(KN,KV),A(KV,KV), &
      
1017.                 SQ(KV,KV),AS(KV,KV),RP(20),RE(20,6),S1(6,3),RS(20,3), &
      
1018.                 DCG(3),RC(20),GA(KV),XA(20),YA(20),BI(3),CI(3),B1(3), &
      
1019.                 RSS(20,3)
      
1020.       COMMON/CA/NP,NE,NM,NR,NI,NL,NG,ND,NC,NA,NN1,DH,NV,NS
      
1021.       COMMON/CC/N,NH,JR(2,800),R(1600)
      
1022.       COMMON/CB/EO,VO,W,T,S,H11,H12,H21,H22,ME(3),BI,CI,ER,TA1,        &
      
1023.                 TA2,NB,L
      
1024.       COMMON/CD/EE(20),SS(20),A1(6),XA,YA,JA(20),SD(20), &
      
1025.                 ED(20),SP(20),EP(20),A2(6,50),KL1,KL2,CH(4),KL
      
1026. !     B1(1)=(CH(4)-CH(2))*(CH(4)-CH(3))/(CH(1)-CH(2))/(CH(1)-CH(3))
      
1027. !     B1(2)=(CH(4)-CH(1))*(CH(4)-CH(3))/(CH(2)-CH(1))/(CH(2)-CH(3))
      
1028. !     B1(3)=(CH(4)-CH(1))*(CH(4)-CH(2))/(CH(3)-CH(1))/(CH(3)-CH(2))
      
1029.       DO 20 I=1,NV
      
1030.       DO 20 J=1,3
      
1031.       RS(I,J)=0.0
      
1032.   20  CONTINUE
      
1033. !     DO 115 IE=KL1,KL2
      
1034. !         JJ=(IE-KL1)/2+1
      
1035. !        JJ=(IE-1)/2+1
      
1036. !    CALL DIV(IE,AE,X,Y,MEO,NM,NP,NE)
      
1037.           CALL DIV(2,AE,X,Y,MEO,NM,NP,NE)
      
1038.       DO 85 J1=1,NV
      
1039.       DO 90 I=1,3
      
1040.       J2=ME(I)
      
1041.       DO 90 K=1,2
      
1042.       K1=2*(I-1)+K
      
1043.       NJ=JR(K,J2)
      
1044.       IF(NJ.GT.0)  GOTO 92
      
1045.       GOTO 94
      
1046.   92  RE(J1,K1)=RD(NJ,J1)
      
1047.       GOTO 90
      
1048.   94  RE(J1,K1)=0.0
      
1049.   90  CONTINUE
      
1050.   85  CONTINUE
      
1051.      WRITE(6,82)((RE(I,J),J=1,6),I=1,NV)
      
1052. 82  FORMAT(25X,'RE(I,J)*****='//(5X,5E14.6))
      
1053.       V1=VO/(1.0-VO)
      
1054.       V2=EO/(1.0-VO*VO)
      
1055.       A11=V2/2.0/S/(1.0-V1*V1)
      
1056.       A22=(1.0-V1)/2.0
      
1057.       DO 100 J=1,3
      
1058.       K1=2*J-1
      
1059.       K2=2*J
      
1060.       S1(K1,1)=A11*BI(J)
      
1061.       S1(K2,1)=A11*CI(J)*V1
      
1062.       S1(K1,2)=A11*BI(J)*V1
      
1063.       S1(K2,2)=A11*CI(J)
      
1064.       S1(K1,3)=A11*CI(J)*A22
      
1065.       S1(K2,3)=A11*BI(J)*A22
      
1066. 100  CONTINUE
      
1067.      WRITE(6,102)((S1(I,J),J=1,3),I=1,6)
      
1068. 102  FORMAT(25X,'S1********='//(5X,5E14.6))
      
1069.       DO 105 I=1,NV
      
1070.       DO 105 K=1,3
      
1071. !      RSS(I,K)=0.0
      
1072.           RS(I,K)=0.0
      
1073.       DO 110 J=1,6
      
1074. !      RSS(I,K)=RSS(I,K)+RE(I,J)*S1(J,K)
      
1075.           RS(I,K)=RS(I,K)+RE(I,J)*S1(J,K)
      
1076. 110  CONTINUE
      
1077. 105  CONTINUE
      
1078. !      LL=2*(L+1)
      
1079. !      J=IE-KL1+1
      
1080.      DO 125 I=1,3
      
1081. !      RSS(LL,I)=RSS(LL,I)+A2(I,J)/XA(LL)
      
1082.           RS(4,I)=RS(4,I)+A1(I)/XA(4)
      
1083. 125  CONTINUE
      
1084. !      DO 128 I=1,NV
      
1085. !      DO 128 J=1,3
      
1086. !      RS(I,J)=RS(I,J)+RSS(I,J)*B1(JJ)/2.0
      
1087. ! 128  CONTINUE
      
1088. ! 115  CONTINUE
      
1089.      WRITE(6,116)((RS(I,J),J=1,3),I=1,NV)
      
1090. 116  FORMAT(25X,'RS**='//(5X,5F12.4))
      
1091.       IF (NI.GT.0) GOTO 120
      
1092.       DCG(1)=0.5+(A1(1)-A1(2))/4/SQRT(((A1(1)-A1(2))/2)**2+A1(3)**2)
      
1093.       DCG(2)=0.5+(A1(2)-A1(1))/4/SQRT(((A1(1)-A1(2))/2)**2+A1(3)**2)
      
1094.       DCG(3)=A1(3)/SQRT(((A1(1)-A1(2))/2)**2+A1(3)**2)
      
1095.       GOTO 130
      
1096. 120  DCG(1)=0.5-(A1(1)-A1(2))/4/SQRT(((A1(1)-A1(2))/2)**2+A1(3)**2)
      
1097.       DCG(2)=0.5-(A1(2)-A1(1))/4/SQRT(((A1(1)-A1(2))/2)**2+A1(3)**2)
      
1098.       DCG(3)=-A1(3)/SQRT(((A1(1)-A1(2))/2)**2+A1(3)**2)
      
1099. 130  DO 145 I=1,NV
      
1100.      RC(I)=0.0
      
1101.       DO 150 J=1,3
      
1102.       RC(I)=RC(I)+RS(I,J)*DCG(J)
      
1103. 150  CONTINUE
      
1104. 145  CONTINUE
      
1105.       IF (NI.GT.0) GOTO 154
      
1106.       DO 155 I=1,NV
      
1107.       RC(I)=-RC(I)
      
1108. 155  CONTINUE
      
1109. !      RC(NV)=1.0+RC(NV)
      
1110.           RC(5)=1.0+RC(5)
      
1111.       GOTO 158
      
1112. ! 154  RC(NV)=RC(NV)+1.0
      
1113. 154  RC(5)=RC(5)+1.0
      
1114. 158  WRITE(6,156)(RC(I),I=1,NV)
      
1115. 156  FORMAT(2X,'PARSURE XX'/(5X,5F12.4))
      
1116.       GG=0.0
      
1117.       IF(NI.GT.0) GOTO 162
      
1118.       GY=XA(NV)-A1(4)
      
1119.       GOTO 161
      
1120. ! 162  GY=XA(NV)+A1(5)
      
1121.   162  GY=XA(5)+A1(5)
      
1122. 161  IF(NS.NE.0) GOTO 400
      
1123.       CALL RI3(GY,RC,GA)
      
1124.       GOTO 440
      
1125. 400  DO 160 I=1,NV
      
1126.       RC(I)=RC(I)*SP(I)
      
1127. 160  CONTINUE
      
1128.       DO 200 I=1,NV
      
1129.       DO 200 J=1,NV
      
1130.       AS(I,J)=0.0
      
1131.       AS(I,J)=SQRT(A(I,I))*SQ(J,I)
      
1132. 200   CONTINUE
      
1133.       DO 220 I=1,NV
      
1134.       RP(I)=0.0
      
1135. 220   CONTINUE
      
1136.       DO 300 I=1,NV
      
1137.       DO 300 J=1,NV
      
1138.       RP(I)=RP(I)+AS(I,J)*RC(J)
      
1139. 300   CONTINUE
      
1140.       WRITE(6,301)(RP(I),I=1,NV)
      
1141. 301   FORMAT(2X,'PARSURE ZZ'/(5X,5F12.4))
      
1142.       DO 330 I=1,NV
      
1143.       GG=GG+RP(I)*RP(I)
      
1144. 330   CONTINUE
      
1145.       GG=SQRT(GG)
      
1146.       DO 164 I=1,NV
      
1147.       GA(I)=-RP(I)/GG
      
1148.       YA(I)=SD(I)*YA(I)/SP(I)+(ED(I)-EP(I))/SP(I)
      
1149. 164  CONTINUE
      
1150.       WRITE(6,190) (YA(I),I=1,NV)
      
1151. 190  FORMAT(10X,'YYYYY'//(10X,5F13.4))
      
1152.       Y1=0.0
      
1153.       DO 165 I=1,NV
      
1154.       Y1=Y1+YA(I)*GA(I)
      
1155. 165  CONTINUE
      
1156.       Y1=Y1+GY/GG
      
1157.       DO 170 I=1,NV
      
1158.       YA(I)=GA(I)*Y1
      
1159. 170  CONTINUE
      
1160.       DO 350 I=1,NV
      
1161.       RC(I)=0.0
      
1162.       DO 350 J=1,NV
      
1163.       RC(I)=RC(I)+AS(J,I)*YA(J)
      
1164. 350  CONTINUE
      
1165.       DO 180 I=1,NV
      
1166.       XA(I)=RC(I)*SP(I)+EP(I)
      
1167. 180  CONTINUE
      
1168. 440  CONTINUE
      
1169.       RETURN
      
1170.       END
      
1171. !     ****************************************************************
      
1172.       SUBROUTINE RI2(GY,AE,X,Y,MEO,RD,GA,A,SQ,AS,KN,KP,KE,KM,KV)
      
1173.       DIMENSION AE(4,KM),X(KP),Y(KP),MEO(2,KE),RD(KN,KV),A(KV,KV), &
      
1174.                 SQ(KV,KV),AS(KV,KV),RP(20),RE(20,6),S1(6,3),RS(20,3), &
      
1175.                 RC(20),GA(KV),XA(20),YA(20),BI(3),CI(3)
      
1176.       COMMON/CA/NP,NE,NM,NR,NI,NL,NG,ND,NC,NA,NN1,DH,NV,NS
      
1177.       COMMON/CC/N,NH,JR(2,800),R(1600)
      
1178.       COMMON/CB/EO,VO,W,T,S,H11,H12,H21,H22,ME(3),BI,CI,ER,TA1,        &
      
1179.                 TA2,NB,L
      
1180.       COMMON/CD/EE(20),SS(20),A1(6),XA,YA,JA(20),SD(20), &
      
1181.                 ED(20),SP(20),EP(20),A2(6,50),KL1,KL2
      
1182.       DO 240 I=1,NV
      
1183. 240  RC(I)=0.0
      
1184.       GY=0.0
      
1185.       NX=NV-2
      
1186.       DO 260 II=KL1,KL2,2
      
1187.       IQ=II-KL1+1
      
1188.       IF(A2(2,IQ)) 50,260,260
      
1189. 50   CALL DIV(II,AE,X,Y,MEO,NM,NP,NE)
      
1190.       DO 85 J1=1,NX
      
1191.       DO 90 I=1,3
      
1192.       J2=ME(I)
      
1193.       DO 90 K=1,2
      
1194.       K1=2*(I-1)+K
      
1195.       NJ=JR(K,J2)
      
1196.       IF(NJ.GT.0) GOTO 92
      
1197.       GOTO 94
      
1198.   92  RE(J1,K1)=RD(NJ,J1)
      
1199.       GOTO 90
      
1200.   94  RE(J1,K1)=0.0
      
1201.   90  CONTINUE
      
1202.   85  CONTINUE
      
1203. !     WRITE(6,82) II,((RE(I,J),J=1,6),I=1,NX)
      
1204. ! 82  FORMAT(20X,'II=',I3,5X,'RE(I,J)====='//(5X,5E14.6))
      
1205.       V1=VO/(1.0-VO)
      
1206.       V2=EO/(1.0-VO*VO)
      
1207.       A11=V2/2.0/S/(1.0-V1*V1)
      
1208.       A22=(1.0-V1)/2.0
      
1209.       DO 100 J=1,3
      
1210.       K1=2*J-1
      
1211.       K2=2*J
      
1212.       S1(K1,1)=A11*BI(J)
      
1213.       S1(K2,1)=A11*CI(J)*V1
      
1214.       S1(K1,2)=A11*BI(J)*V1
      
1215.       S1(K2,2)=A11*CI(J)
      
1216.       S1(K1,3)=A11*CI(J)*A22
      
1217.       S1(K2,3)=A11*BI(J)*A22
      
1218. 100  CONTINUE
      
1219. !     WRITE(6,102)((S1(I,J),J=1,3),I=1,6)
      
1220. !102  FORMAT(25X,'S1************='//(5X,5F12.4))
      
1221.       DO 105 I=1,NX
      
1222.       DO 105 K=1,3
      
1223.       RS(I,K)=0.0
      
1224.       DO 110 J=1,6
      
1225.       RS(I,K)=RS(I,K)+RE(I,J)*S1(J,K)
      
1226. 110  CONTINUE
      
1227. 105  CONTINUE
      
1228.       LL=2*(L+1)
      
1229.       DO 115 I=1,3
      
1230.       RS(LL,I)=RS(LL,I)+A2(I,IQ)/XA(LL)
      
1231. 115  CONTINUE
      
1232. !     WRITE(6,116)((RS(I,J),J=1,3),I=1,NX)
      
1233. !116  FORMAT(25X,'RS(I,J)=*****'//(5X,5F12.4))
      
1234.       CB=ABS(CI(1))
      
1235.       IF(NB.NE.2) GOTO 55
      
1236.       CB=ABS(CI(2))
      
1237. 55   DO 216 I=1,NX
      
1238.       RC(I)=RC(I)-CB*(XA(NX+1)*RS(I,2)+RS(I,3))
      
1239. 216  CONTINUE
      
1240.       RC(NX+1)=RC(NX+1)+ABS(A2(2,IQ))*CB
      
1241.       RC(NX+2)=RC(NX+2)+CB
      
1242.       GY=GY+CB*(XA(NX+1)*ABS(A2(2,IQ))+XA(NX+2)-A2(3,IQ))
      
1243. !     WRITE (6,266) (RC(I),I=1,NV)
      
1244. !266  FORMAT(2X,'PARSURE XX'/(5X,6F11.4))
      
1245. 260  CONTINUE
      
1246.       IF(NS.NE.0) GOTO 122
      
1247.       CALL RI3(GY,RC,GA)
      
1248.       GOTO 550
      
1249. 122  DO 160 I=1,NV
      
1250. 160  RC(I)=RC(I)*SP(I)
      
1251.       DO 200 I=1,NV
      
1252.       DO 200 J=1,NV
      
1253.       AS(I,J)=0.0
      
1254.       AS(I,J)=SQRT(A(I,I))*SQ(J,I)
      
1255. 200  CONTINUE
      
1256.       DO 220 I=1,NV
      
1257.       RP(I)=0.0
      
1258. 220  CONTINUE
      
1259.       DO 300 I=1,NV
      
1260.       DO 300 J=1,NV
      
1261.       RP(I)=RP(I)+AS(I,J)*RC(J)
      
1262. 300  CONTINUE
      
1263. !     WRITE(6,262) (RP(I),I=1,NV)
      
1264. !262  FORMAT(2X,'PARSURE ZZ'/(5X,6F11.4))
      
1265.       GG=0.0
      
1266.       DO 330 I=1,NV
      
1267.       GG=GG+RP(I)*RP(I)
      
1268. 330  CONTINUE
      
1269.       GG=SQRT(GG)
      
1270.       DO 164 I=1,NV
      
1271.       GA(I)=-RP(I)/GG
      
1272.       YA(I)=SD(I)*YA(I)/SP(I)+(ED(I)-EP(I))/SP(I)
      
1273. 164  CONTINUE
      
1274.       WRITE(6,190) (YA(I),I=1,NV)
      
1275. 190  FORMAT(10X,'YYYYY'//(10X,5F13.4))
      
1276.       Y1=0.0
      
1277.       DO 165 I=1,NV
      
1278.       Y1=Y1+YA(I)*GA(I)
      
1279. 165  CONTINUE
      
1280.       Y1=Y1+GY/GG
      
1281.       DO 170 I=1,NV
      
1282.       YA(I)=GA(I)*Y1
      
1283. 170  CONTINUE
      
1284.       DO 350 I=1,NV
      
1285.       RC(I)=0.0
      
1286.       DO 350 J=1,NV
      
1287.       RC(I)=RC(I)+AS(J,I)*YA(J)
      
1288. 350  CONTINUE
      
1289.       DO 180 I=1,NV
      
1290.       XA(I)=RC(I)*SP(I)+EP(I)
      
1291. 180  CONTINUE
      
1292. 550  CONTINUE
      
1293.       RETURN
      
1294.       END
      
1295. !     *****************************************************************
      
1296.       SUBROUTINE RI3(GY,RC,GA)
      
1297.       DIMENSION RC(20),GA(20),XA(20),YA(20)
      
1298.       COMMON/CA/NP,NE,NM,NR,NI,NL,NG,ND,NC,NA,NN1,DH,NV,NS
      
1299.       COMMON/CD/EE(20),SS(20),A1(6),XA,YA,JA(20),SD(20), &
      
1300.                 ED(20),SP(20),EP(20)
      
1301.       GG=0.0
      
1302.       DO 160 I=1,NV
      
1303.       RC(I)=RC(I)*SP(I)
      
1304.       GG=GG+RC(I)**2
      
1305. 160  CONTINUE
      
1306.       GG=SQRT(GG)
      
1307.       DO 164 I=1,NV
      
1308.       GA(I)=-RC(I)/GG
      
1309.       YA(I)=SD(I)*YA(I)/SP(I)+(ED(I)-EP(I))/SP(I)
      
1310. 164  CONTINUE
      
1311.       WRITE(6,190)(YA(I),I=1,NV)
      
1312. 190  FORMAT(10X,'YYYYY'//(10X,5F13.4))
      
1313.       Y1=0.0
      
1314.       DO 165 I=1,NV
      
1315.       Y1=Y1+YA(I)*GA(I)
      
1316. 165  CONTINUE
      
1317.       Y1=Y1+GY/GG
      
1318.       DO 170 I=1,NV
      
1319.       YA(I)=GA(I)*Y1
      
1320. 170  CONTINUE
      
1321.       DO 180 I=1,NV
      
1322.       XA(I)=YA(I)*SP(I)+EP(I)
      
1323. 180  CONTINUE
      
1324.       RETURN
      
1325.       END
      
1326. !     *****************************************************************
      
1327.       SUBROUTINE COV(N1,A,SQ,EO)
      
1328.       DIMENSION A(N1,N1),SQ(N1,N1)
      
1329.       READ(5,*)((A(I,J),I=1,N1),J=1,N1)
      
1330.       WRITE(6,6)((A(I,J),I=1,N1),J=1,N1)
      
1331.    6  FORMAT(10X,'A(I,J)='/(5X,5F13.5))
      
1332.       DO 10 I=1,N1
      
1333.       DO 10 J=1,N1
      
1334.       SQ(I,J)=0.0
      
1335. 10   CONTINUE
      
1336.       DO 15 I=1,N1
      
1337.       SQ(I,I)=1.0
      
1338. 15   CONTINUE
      
1339.       G=0.0
      
1340.       DO 40 I=2,N1
      
1341.       I1=I-1
      
1342.       DO 40 J=1,I1
      
1343.       G=G+2.0*A(I,J)**2
      
1344. 40   CONTINUE
      
1345. !     WRITE(6,42) G
      
1346. !42   FORMAT(5X,'G=',F13.5)
      
1347.       S1=SQRT(G)
      
1348.       FN=FLOAT(N1)
      
1349.       S2=EO*S1/FN
      
1350. !     WRITE (6,45) S2
      
1351.       S3=S1
      
1352. !45   FORMAT (10X,'S2=',F13.5)
      
1353.       L=0
      
1354. 50   S3=S3/FN
      
1355. !     WRITE (6,55) S3
      
1356. !55   FORMAT (10X,'S3=',F13.5)
      
1357. 60   DO 130 IQ=2,N1
      
1358.       IQ1=IQ-1
      
1359.       DO 130 IP=1,IQ1
      
1360.       IF(ABS(A(IP,IQ)).LT.S3) GOTO 130
      
1361.       L=1
      
1362.       V1=A(IP,IP)
      
1363.       V2=A(IP,IQ)
      
1364.       V3=A(IQ,IQ)
      
1365.       U=0.5*(V1-V3)
      
1366. !     WRITE(6,65) V1,V2,V3,U
      
1367. !65   FORMAT(10X,'V1*V2*V3*U=',4F13.5)
      
1368.       IF(U) 70,80,90
      
1369. 70   B=-1.0
      
1370.       GOTO 100
      
1371. 90   B=1.0
      
1372. 100   G=-B*V2/SQRT(V2*V2+U*U)
      
1373.       GOTO 105
      
1374. 80   IF(V2.GT.S3) GOTO 85
      
1375.       G=1.0
      
1376.       GOTO 105
      
1377. 85   G=-1.0
      
1378. 105  ST=G/SQRT(2.0*(1.0+SQRT(1.0-G*G)))
      
1379.       CT=SQRT(1.0-ST*ST)
      
1380. !     WRITE(6,107) ST,CT
      
1381. !107  FORMAT(10X,'ST*CT=',2F13.5)
      
1382.       DO 110 I=1,N1
      
1383.       G=A(I,IP)*CT-A(I,IQ)*ST
      
1384.       A(I,IQ)=A(I,IP)*ST+A(I,IQ)*CT
      
1385.       A(I,IP)=G
      
1386.       G=SQ(I,IP)*CT-SQ(I,IQ)*ST
      
1387.       SQ(I,IQ)=SQ(I,IP)*ST+SQ(I,IQ)*CT
      
1388.       SQ(I,IP)=G
      
1389. 110  CONTINUE
      
1390.       DO 120 I=1,N1
      
1391.       A(IP,I)=A(I,IP)
      
1392.       A(IQ,I)=A(I,IQ)
      
1393. 120  CONTINUE
      
1394.       G=2.0*V2*ST*CT
      
1395.       A(IP,IP)=V1*CT*CT+V3*ST*ST-G
      
1396.       A(IQ,IQ)=V1*ST*ST+V3*CT*CT+G
      
1397.       A(IP,IQ)=(V1-V3)*ST*CT+V2*(CT*CT-ST*ST)
      
1398.       A(IQ,IP)=A(IP,IQ)
      
1399. 130  CONTINUE
      
1400.       IF(L-1) 150,140,150
      
1401. 140  L=0
      
1402.       GOTO 60
      
1403. 150  WRITE(6,145)((A(I,J),I=1,N1),J=1,N1)
      
1404.       WRITE(6,142)((SQ(I,J),I=1,N1),J=1,N1)
      
1405. 142  FORMAT(20X,'SQ(I,J)='/(5X,5F13.5))
      
1406. 145  FORMAT(20X,'A(I,J)='/(5X,5F13.5))
      
1407.       IF(S3.GT.S2) GOTO 50
      
1408.       RETURN
      
1409.       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

schultz

请问楼主：该随机有限元程序是不是 武清玺编著的那本书上的程序？<BR>

linqus

谢谢先，
有以下疑问：
（1）没有有关程序注释的readme文件。
（2）随时.90后缀的f90文件，却完全是f77标准的程序。
（3）敬请楼主说明程序的配套书籍^_^

再次感谢。

jieli

谢谢，能不能说得详细点儿。

yoyoxo

就是武清玺那本书上的程序

linqus

<P>查了一下，图书馆没有；超星也没有。<BR>不知道哪位可以提供一下电子版的，<BR>非常感谢！！</P>

haha6877

是啊，急需电子版的，哪位好心人做点好事啊

小林子

谢谢

supervb

汗！看来要靠自己写注释了

幻觉

谢谢楼主,辛苦了

coldwind042

ou的神！终于弄到武老师的这个程序啦！拜谢！

coldwind042

怎么运行时出现问题啊？哪位运行通过了的啊？

ChaChing

好些年前的程式, 个人以为不见得编译仍相同, 可能需稍为修改下吧!

wei_x

不能编译，而且是fortran77的代码，看起来太累了··
要是fortran90的就太好了

ccdong

要是有电子版的就好了