Télécharger chabok.eso

Retour à la liste

Numérotation des lignes :

  1. C CHABOK SOURCE CHAT 05/01/12 21:53:52 5004
  2. SUBROUTINE CHABOK(YUNG,XNU,IA,EI,
  3. 1 XMAT,ALPHA1,JA,IBOU,SI,DEPS,EPST,EPSTAR,AMTRI,KERRE,
  4. 2 SN,ALPHA2,NUMCHA,ecou,necou)
  5. C***********************************************************************
  6. C INTEGRATION MODELE DE CHABOCHE
  7. C***********************************************************************
  8. IMPLICIT INTEGER(I-N)
  9. IMPLICIT REAL *8(A-H,O-Z)
  10. DIMENSION XMAT(*),EI(*),ALPHA1(*),ALPHA2(*)
  11. DIMENSION AMTRI(*),YK5(18),DYK1(18),DYK2(18),DYK3(18),
  12. . DYK4(18),EINT(6),W1INT(6),W2INT(6)
  13. C
  14. * Commun ECOU: sert de fourre-tout pour les tableaux
  15. *
  16. SEGMENT ECOU
  17. *** COMMON/ECOU/TEST,ALFAH,
  18. REAL*8 TEST, ALFAH,
  19. 1 HPAS,TEMPS,COVNMS(6),VECPRO(9),VALPRO(6),
  20. 2 CVNMSD(12),STOT(6),SIGEL(6),DSIGP(6),SIGT(6),W1(6),W2(6),
  21. 1 DALPHA(6),DSIGO(6),E(12),XINV(3),
  22. 2 SIPLAD(6),DSIGP0(6),TET,TETI
  23. ENDSEGMENT
  24.  
  25. C COMMON/ECOU/TEST,ALFAH,
  26. C 1 HPAS,TEMPS,COVNMS(6),VECPRO(9),VALPRO(6),
  27. C 2 CVNMSD(12),STOT(6),SIGEL(6),DSIGP(6),SIGT(6),W1(6),W2(6),
  28. C 1 DALPHA(6),DSIGO(6),E(12),XINV(3),
  29. C 2 SIPLAD(6),DSIGP0(6),TET,TETI
  30. C
  31. *
  32. * Commun NECOU utilisé dans ECOINC
  33. *
  34. SEGMENT NECOU
  35. * COMMON/NECOU/NCOURB,IPLAST,IT,IMAPLA,ISOTRO,
  36. INTEGER NCOURB,IPLAST,IT,IMAPLA,ISOTRO,
  37. . ITYP,IFOUR,IFLUAG,
  38. . ICINE,ITHER,IFLUPL,ICYCL,IBI,
  39. . JFLUAG,KFLUAG,LFLUAG,
  40. . IRELAX,JNTRIN,MFLUAG,JSOUFL,JGRDEF
  41. ENDSEGMENT
  42. C COMMON/NECOU/NCOURB,IPLAST,IT,IMAPLA,ISOTRO,
  43. C . ITYP,IFOUR,IFLUAG,
  44. C . ICINE,ITHER,IFLUPL,ICYCL,IBI,
  45. C . JFLUAG,KFLUAG,LFLUAG,
  46. C . IRELAX,JNTRIN,MFLUAG,JSOUFL,JGRDEF
  47. C
  48. C EN SORTIE
  49. C KERRE = 2 SI ON ATTEINT LE NOMBRE MAX D'ITERATIONS INTERNES
  50. C
  51. DATA ITMAX/15/
  52. DATA IDECOU/20/
  53. DATA RAC2/0.707106781D00/
  54. PREDEF=1.D-6
  55. PRESIG=YUNG*PREDEF
  56. PREPH=1.D-3
  57. PRERUN=TEST
  58. PREMIN=3.D-1*PRERUN
  59. C
  60. C INITIALISATION A 0. DES VECTEURS POUR TREANOR
  61. C
  62. NVY=18
  63. DO 3341 IB=1,NVY
  64. YK5(IB)=0.
  65. DYK1(IB)=0.
  66. DYK2(IB)=0.
  67. DYK3(IB)=0.
  68. DYK4(IB)=0.
  69. 3341 CONTINUE
  70. A2=0.
  71. C2=0.
  72. B=0.
  73. RM=0.
  74. PHI=1.D00
  75. ICOD=0
  76. CALL CHALIM(EPSTAR,R,XMAT,TET,ICOD,A1,C1,A2,C2,R0,RM,B,
  77. . PHI,PSI,OME,ICENT2,IDIAM,NUMCHA)
  78. ELT=YUNG/(1.+XNU)
  79. IF(ITHER.NE.0) ELT=ELT*EI(IA)/YUNG
  80. G=ELT*0.5
  81. SAC1=A1*C1
  82. AC1=SAC1*0.66666667
  83. SAC2=A2*C2
  84. AC2=SAC2*0.66666667
  85. SAC12=SAC1+SAC2
  86. AC12=AC1+AC2
  87. GO TO (101,102,102,103,105,102,103,103,109,999,
  88. . 999,999,102,101),ITYP
  89. 101 E(7)=ELT*(1.-XNU)/(1.-2.*XNU)
  90. E(8)=ELT*XNU/(1.-2.*XNU)
  91. E(9)=AC12
  92. E(10)=AC1
  93. E(11)=AC2
  94. GO TO 200
  95. 102 E(7)=ELT/(1.-XNU)
  96. E(8)=E(7)*XNU
  97. E(9)=SAC12
  98. E(10)=SAC1
  99. E(11)=SAC2
  100. GO TO 200
  101. 103 E(7)=ELT*(1.+XNU)
  102. E(8)=0.
  103. E(9)=SAC12
  104. E(10)=SAC1
  105. E(11)=SAC2
  106. GO TO 200
  107. 105 CONTINUE
  108. 109 CONTINUE
  109. C
  110. C INITIALISATIONS
  111. C
  112. 200 ITER=0
  113. CALL CHAREM(IDIAM,R0,RM,B,R,EPSTAR,EPST,DEPS,STOT,
  114. .W1,W2,E,IBOU,ICENT2,JA,ALPHA1,ALPHA2,PHI,PSI,OME)
  115. SI=VONMIS(E,ITYP,ALFAH,COVNMS)
  116. C
  117. C PREMIER DEPS
  118. C
  119. ICAS=1
  120. CALL CHAFON(EPST,AMTRI,AMTRI(7),AMTRI(13),DYK1,
  121. . SI,C1,C2,ITYP,ICENT2,IDIAM,G,R,IBOU,ELT,
  122. . DEPS,R0,RM,B,ICAS,WEP,EINT,W1INT,W2INT,PSI,OME,ecou)
  123. 670 DEPS0=DEPS
  124. DEPSI=DEPS
  125. IF(ITER.EQ.0) GO TO 666
  126. CALL CHAREM(IDIAM,R0,RM,B,R,EPSTAR,EPST,DBID,STOT,
  127. . W1,W2,E,IBOU,ICENT2,JA,ALPHA1,ALPHA2,PHI,PSI,OME)
  128. SI=VONMIS(E,ITYP,ALFAH,COVNMS)
  129. 666 CONTINUE
  130. C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  131. C ITERATIONS INTERNES
  132. C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  133. 555 ITER=ITER+1
  134. C
  135. C CALCUL PAR LA METHODE DE TREANOR
  136. C
  137. DO 413 IB=1,IBOU
  138. EINT(IB)=E(IB)
  139. W1INT(IB)=W1(IB)
  140. IF(ICENT2.EQ.0) GO TO 413
  141. W2INT(IB)=W2(IB)
  142. 413 CONTINUE
  143. SINT=SI
  144. DEDEP=0.
  145. HDEP=DEPSI
  146. IF(ITER.EQ.1) HDEP0=HDEP
  147. EPSINT=EPST
  148. C-----------------------------------------------------------------------
  149. C BOUCLE SUR LES SOUS-PAS
  150. C-----------------------------------------------------------------------
  151. DO 416 IDECO=1,IDECOU
  152. JESSAY=0
  153. DEPRES=DEPSI-DEDEP
  154. IF(ABS(HDEP).GT.ABS(DEPRES)) HDEP=DEPRES
  155. IF(IDECO.EQ.1) GO TO 417
  156. C
  157. 424 CONTINUE
  158. JESSAY=JESSAY+1
  159. DO 418 IB=1,IBOU
  160. E(IB)=EINT(IB)
  161. W1(IB)=W1INT(IB)
  162. IF(ICENT2.EQ.0) GO TO 418
  163. W2(IB)=W2INT(IB)
  164. 418 CONTINUE
  165. SI=SINT
  166. 417 ICAS=2
  167. WEP=0.5D00
  168. CALL CHAFON( EPSINT ,AMTRI,AMTRI(7),AMTRI(13),DYK1,
  169. . SI,C1,C2,ITYP,ICENT2,IDIAM,G,R,IBOU,ELT,HDEP,
  170. . R0,RM,B,ICAS,WEP,EINT,W1INT,W2INT,PSI,OME,ecou)
  171. C
  172. EPST=EPSINT+0.5D00*HDEP
  173. CALL CHAFON( EPST ,AMTRI,AMTRI(7),AMTRI(13),DYK2,
  174. . SI,C1,C2,ITYP,ICENT2,IDIAM,G,R,IBOU,ELT,HDEP,
  175. . R0,RM,B,ICAS,WEP,EINT,W1INT,W2INT,PSI,OME,ecou)
  176. C
  177. WEP=1.D00
  178. CALL CHAFON( EPST ,AMTRI,AMTRI(7),AMTRI(13),DYK3,
  179. . SI,C1,C2,ITYP,ICENT2,IDIAM,G,R,IBOU,ELT,HDEP,
  180. . R0,RM,B,ICAS,WEP,EINT,W1INT,W2INT,PSI,OME,ecou)
  181. C
  182. PHMAX=0.
  183. DO 4101 IB=1,NVY
  184. IF(ABS(DYK2(IB)-DYK1(IB)).LT.PRESIG) GO TO 4101
  185. PH=2.*(DYK2(IB)-DYK3(IB))/(DYK2(IB)-DYK1(IB))
  186. IF(PHMAX.LT.PH) PHMAX=PH
  187. 4101 CONTINUE
  188. IF(PHMAX.LE.0.) GO TO 4102
  189. IF(PHMAX.GE.PREPH) GO TO 4123
  190. XL1=1.-PHMAX*0.5
  191. XL2=0.5-PHMAX*0.166666666666666667
  192. XL3=0.166666666666666667-PHMAX*0.04166666666666666667
  193. XL4=-PHMAX*0.166666666666666667
  194. GO TO 4103
  195. 4123 CONTINUE
  196. XL1=(1.-EXP(-PHMAX))/PHMAX
  197. XL2=(1.-XL1)/PHMAX
  198. XL3=(0.5-XL2)/PHMAX
  199. XL4=XL1-2.*XL2
  200. GO TO 4103
  201. 4102 PHMAX=0.
  202. XL1=1.
  203. XL2=0.5
  204. XL3=0.166666666666666667
  205. XL4=0.
  206. 4103 CONTINUE
  207. XM1=-XL2+4.*XL3
  208. XM2=2.*(XL2-2.*XL3)
  209. XM3=4.*XL3-3.*XL2
  210. C
  211. C NOUVELLE ESTIMATION DE LA SOLUTION
  212. C
  213. DO 4104 IB=1,IBOU
  214. E(IB)=EINT(IB)+2.*XL2*DYK3(IB)+XL2*PHMAX*DYK2(IB)+XL4*DYK1(IB)
  215. YK5(IB)=E(IB)
  216. W1(IB)=W1INT(IB)+2.*XL2*DYK3(6+IB)
  217. . +XL2*PHMAX*DYK2(6+IB)+XL4*DYK1(6+IB)
  218. YK5(6+IB)=W1(IB)
  219. IF(ICENT2.EQ.0) GO TO 4104
  220. W2(IB)=W2INT(IB)+2.*XL2*DYK3(12+IB)
  221. . +XL2*PHMAX*DYK2(12+IB)+XL4*DYK1(12+IB)
  222. YK5(12+IB)=W2(IB)
  223. 4104 CONTINUE
  224. SI=VONMIS(E,ITYP,ALFAH,COVNMS)
  225. C
  226. EPST=EPSINT+HDEP
  227. WEP=1.D00
  228. CALL CHAFON( EPST ,AMTRI,AMTRI(7),AMTRI(13),DYK4,
  229. . SI,C1,C2,ITYP,ICENT2,IDIAM,G,R,IBOU,ELT,HDEP,
  230. . R0,RM,B,ICAS,WEP,EINT,W1INT,W2INT,PSI,OME,ecou)
  231. C
  232. C ESTIMATION FINALE DE LA SOLUTION
  233. C
  234. FAC1=1.+PHMAX*XM3+2.*PHMAX*XM2
  235. FAC2=XL1+XM3+0.5*XM2*PHMAX
  236. FAC3=XM2*(1.+0.5*PHMAX)
  237. DO 4105 IB=1,IBOU
  238. E(IB)=EINT(IB)*FAC1+DYK1(IB)*FAC2+DYK2(IB)*FAC3
  239. . +DYK3(IB)*XM2+DYK4(IB)*XM1+YK5(IB)*XM1*PHMAX
  240. SIGT(IB)=E(IB)
  241. W1(IB)=W1INT(IB)*FAC1+DYK1(6+IB)*FAC2
  242. . +DYK2(6+IB)*FAC3 +DYK3(6+IB)*XM2 +DYK4(6+IB)*XM1
  243. . +YK5(6+IB)*XM1*PHMAX
  244. IF(ICENT2.EQ.0) GO TO 4105
  245. W2(IB)= W2INT(IB)*FAC1 +DYK1(12+IB)*FAC2
  246. . +DYK2(12+IB)*FAC3 +DYK3(12+IB)*XM2 +DYK4(12+IB)*XM1
  247. . +YK5(12+IB)*XM1*PHMAX
  248. 4105 CONTINUE
  249. SI=VONMIS(SIGT,ITYP,ALFAH,COVNMS)
  250. C
  251. C CALCUL DE L'ERREUR
  252. C
  253. ERMAX=0.
  254. DO 4106 IB=1,IBOU
  255. ERM=ABS(SIGT(IB)-YK5(IB))/MAX(PRESIG,
  256. . ABS(SIGT(IB)-EINT(IB)))
  257. IF(ERMAX.LT.ERM) ERMAX=ERM
  258. ERM=ABS(W1(IB)-YK5(6+IB))/MAX(PRESIG,
  259. . ABS(W1(IB)-W1INT(IB)))
  260. IF(ERMAX.LT.ERM) ERMAX=ERM
  261. IF(ICENT2.EQ.0) GO TO 4106
  262. ERM=ABS(W2(IB)-YK5(12+IB))/MAX(PRESIG,
  263. . ABS(W2(IB)- W2INT(IB)))
  264. IF(ERMAX.LT.ERM) ERMAX=ERM
  265. 4106 CONTINUE
  266. C
  267. IF(ERMAX.LE.PRERUN) GO TO 1240
  268. C-----------------------------------------------------------------------
  269. C POUR ITER >1 ET IDECO=1 , ON REPREND LE HDEP0 SI ON N'A PAS REUSSI
  270. C A INTEGRER TOUT LE DEPSI EN UN SEUL SOUS-PAS
  271. C-----------------------------------------------------------------------
  272. HDEP=HDEP*0.5D00
  273. IF(JESSAY.GT.1) GO TO 424
  274. IF(IDECO.EQ.1.AND.ITER.GT.1.AND.ABS(HDEP).GT.ABS(HDEP0))
  275. . HDEP=HDEP0
  276. GO TO 424
  277. 1240 DEDEP=DEDEP+HDEP
  278. C-----------------------------------------------------------------------
  279. C SI ON A FINI , ON SORT EN 514 . SINON ON CONTINUE
  280. C-----------------------------------------------------------------------
  281. IF(ABS(DEDEP).GE.ABS(DEPSI)) GO TO 514
  282. IF(ERMAX.LE.PREMIN) HDEP=2.D00*HDEP
  283. EPSINT=EPST
  284. SINT=SI
  285. DO 502 IB=1,IBOU
  286. EINT(IB)=SIGT(IB)
  287. W1INT(IB)=W1(IB)
  288. IF(ICENT2.EQ.0) GO TO 502
  289. W2INT(IB)=W2(IB)
  290. 502 CONTINUE
  291. C.......................................................................
  292. C FIN DE LA BOUCLE IDECO
  293. C.......................................................................
  294. 416 CONTINUE
  295. C
  296. C ON CONTINUE AVEC UN DEPS DIMINUE
  297. C
  298. DEPS=DEPS-DEPSI+DEDEP
  299. C-----------------------------------------------------------------------
  300. C FIN DE TREANOR
  301. C-----------------------------------------------------------------------
  302. 514 CONTINUE
  303. HDEP0=HDEP
  304. IF(ABS(SI-R)/R.LT.TEST) GO TO 57
  305. IF(ITER.GE.ITMAX) GO TO 56
  306. C
  307. C CALCUL DE LA PENTE ET DU DEPSI
  308. C
  309. ICAS=3
  310. WEP=0.0D00
  311. CALL CHAFON(EPST,AMTRI,AMTRI(7),AMTRI(13),DYK4,
  312. . SI,C1,C2,ITYP,ICENT2,IDIAM,G,R,IBOU,ELT,DEPS,
  313. . R0,RM,B,ICAS,WEP,E,W1,W2,PSI,OME,ecou)
  314. DEPSI=(SI-R)/WEP
  315. DEPS=DEPS+DEPSI
  316. IF(DEPS.LT.0.) GO TO 54
  317. GO TO 555
  318. C
  319. 54 DEPS=0.5D00*DEPS0
  320. GO TO 670
  321. C
  322. 56 KERRE=2
  323. 57 JX=JA
  324. DO 301 I=1,IBOU
  325. SIGEL(I)=SIGT(I)
  326. DALPHA(I)=W1(I)-ALPHA1(JX)
  327. IF(ICENT2.EQ.0) GO TO 301
  328. DALPHA(I)=DALPHA(I)+W2(I)
  329. 301 JX=JX+1
  330. C
  331. C LES DEUX CENTRES SONT CUMULES DANS ALPHA1
  332. C
  333. SN=R
  334. C
  335. C MISE A JOUR DES CENTRES DES SPHERES
  336. C
  337. JX=JA
  338. DO 303 I=1,IBOU
  339. ALPHA1(JX)=W1(I)
  340. IF(ICENT2.EQ.0) GO TO 303
  341. ALPHA1(JX)=ALPHA1(JX)+W2(I)
  342. ALPHA2(JX)=W2(I)
  343. 303 JX=JX+1
  344. RETURN
  345. 999 WRITE(6,7999)
  346. 7999 FORMAT('0 CHABOK - CAS NON IMPLEMENTE '/)
  347. RETURN
  348. END
  349.  
  350.  
  351.  
  352.  
  353.  
  354.  
  355.  
  356.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales