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=VONMIS(E,ITYP,ALFAH,COVNMS)
  180. C
  181. EPST=EPSINT+HDEP
  182. WEP=1.D00
  183. CALL CHAFON( EPST ,AMTRI,AMTRI(7),AMTRI(13),DYK4,
  184. . SI,C1,C2,ITYP,ICENT2,IDIAM,G,R,IBOU,ELT,HDEP,
  185. . R0,RM,B,ICAS,WEP,EINT,W1INT,W2INT,PSI,OME,ecou)
  186. C
  187. C ESTIMATION FINALE DE LA SOLUTION
  188. C
  189. FAC1=1.+PHMAX*XM3+2.*PHMAX*XM2
  190. FAC2=XL1+XM3+0.5*XM2*PHMAX
  191. FAC3=XM2*(1.+0.5*PHMAX)
  192. DO 4105 IB=1,IBOU
  193. E(IB)=EINT(IB)*FAC1+DYK1(IB)*FAC2+DYK2(IB)*FAC3
  194. . +DYK3(IB)*XM2+DYK4(IB)*XM1+YK5(IB)*XM1*PHMAX
  195. SIGT(IB)=E(IB)
  196. W1(IB)=W1INT(IB)*FAC1+DYK1(6+IB)*FAC2
  197. . +DYK2(6+IB)*FAC3 +DYK3(6+IB)*XM2 +DYK4(6+IB)*XM1
  198. . +YK5(6+IB)*XM1*PHMAX
  199. IF(ICENT2.EQ.0) GO TO 4105
  200. W2(IB)= W2INT(IB)*FAC1 +DYK1(12+IB)*FAC2
  201. . +DYK2(12+IB)*FAC3 +DYK3(12+IB)*XM2 +DYK4(12+IB)*XM1
  202. . +YK5(12+IB)*XM1*PHMAX
  203. 4105 CONTINUE
  204. SI=VONMIS(SIGT,ITYP,ALFAH,COVNMS)
  205. C
  206. C CALCUL DE L'ERREUR
  207. C
  208. ERMAX=0.
  209. DO 4106 IB=1,IBOU
  210. ERM=ABS(SIGT(IB)-YK5(IB))/MAX(PRESIG,
  211. . ABS(SIGT(IB)-EINT(IB)))
  212. IF(ERMAX.LT.ERM) ERMAX=ERM
  213. ERM=ABS(W1(IB)-YK5(6+IB))/MAX(PRESIG,
  214. . ABS(W1(IB)-W1INT(IB)))
  215. IF(ERMAX.LT.ERM) ERMAX=ERM
  216. IF(ICENT2.EQ.0) GO TO 4106
  217. ERM=ABS(W2(IB)-YK5(12+IB))/MAX(PRESIG,
  218. . ABS(W2(IB)- W2INT(IB)))
  219. IF(ERMAX.LT.ERM) ERMAX=ERM
  220. 4106 CONTINUE
  221. C
  222. IF(ERMAX.LE.PRERUN) GO TO 1240
  223. C-----------------------------------------------------------------------
  224. C POUR ITER >1 ET IDECO=1 , ON REPREND LE HDEP0 SI ON N'A PAS REUSSI
  225. C A INTEGRER TOUT LE DEPSI EN UN SEUL SOUS-PAS
  226. C-----------------------------------------------------------------------
  227. HDEP=HDEP*0.5D00
  228. IF(JESSAY.GT.1) GO TO 424
  229. IF(IDECO.EQ.1.AND.ITER.GT.1.AND.ABS(HDEP).GT.ABS(HDEP0))
  230. . HDEP=HDEP0
  231. GO TO 424
  232. 1240 DEDEP=DEDEP+HDEP
  233. C-----------------------------------------------------------------------
  234. C SI ON A FINI , ON SORT EN 514 . SINON ON CONTINUE
  235. C-----------------------------------------------------------------------
  236. IF(ABS(DEDEP).GE.ABS(DEPSI)) GO TO 514
  237. IF(ERMAX.LE.PREMIN) HDEP=2.D00*HDEP
  238. EPSINT=EPST
  239. SINT=SI
  240. DO 502 IB=1,IBOU
  241. EINT(IB)=SIGT(IB)
  242. W1INT(IB)=W1(IB)
  243. IF(ICENT2.EQ.0) GO TO 502
  244. W2INT(IB)=W2(IB)
  245. 502 CONTINUE
  246. C.......................................................................
  247. C FIN DE LA BOUCLE IDECO
  248. C.......................................................................
  249. 416 CONTINUE
  250. C
  251. C ON CONTINUE AVEC UN DEPS DIMINUE
  252. C
  253. DEPS=DEPS-DEPSI+DEDEP
  254. C-----------------------------------------------------------------------
  255. C FIN DE TREANOR
  256. C-----------------------------------------------------------------------
  257. 514 CONTINUE
  258. HDEP0=HDEP
  259. IF(ABS(SI-R)/R.LT.TEST) GO TO 57
  260. IF(ITER.GE.ITMAX) GO TO 56
  261. C
  262. C CALCUL DE LA PENTE ET DU DEPSI
  263. C
  264. ICAS=3
  265. WEP=0.0D00
  266. CALL CHAFON(EPST,AMTRI,AMTRI(7),AMTRI(13),DYK4,
  267. . SI,C1,C2,ITYP,ICENT2,IDIAM,G,R,IBOU,ELT,DEPS,
  268. . R0,RM,B,ICAS,WEP,E,W1,W2,PSI,OME,ecou)
  269. DEPSI=(SI-R)/WEP
  270. DEPS=DEPS+DEPSI
  271. IF(DEPS.LT.0.) GO TO 54
  272. GO TO 555
  273. C
  274. 54 DEPS=0.5D00*DEPS0
  275. GO TO 670
  276. C
  277. 56 KERRE=2
  278. 57 JX=JA
  279. DO 301 I=1,IBOU
  280. SIGEL(I)=SIGT(I)
  281. DALPHA(I)=W1(I)-ALPHA1(JX)
  282. IF(ICENT2.EQ.0) GO TO 301
  283. DALPHA(I)=DALPHA(I)+W2(I)
  284. 301 JX=JX+1
  285. C
  286. C LES DEUX CENTRES SONT CUMULES DANS ALPHA1
  287. C
  288. SN=R
  289. C
  290. C MISE A JOUR DES CENTRES DES SPHERES
  291. C
  292. JX=JA
  293. DO 303 I=1,IBOU
  294. ALPHA1(JX)=W1(I)
  295. IF(ICENT2.EQ.0) GO TO 303
  296. ALPHA1(JX)=ALPHA1(JX)+W2(I)
  297. ALPHA2(JX)=W2(I)
  298. 303 JX=JX+1
  299. RETURN
  300. 999 WRITE(6,7999)
  301. 7999 FORMAT('0 CHABOK - CAS NON IMPLEMENTE '/)
  302. RETURN
  303. END
  304.  
  305.  
  306.  
  307.  
  308.  
  309.  
  310.  
  311.  

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