Télécharger fatig3.eso

Retour à la liste

Numérotation des lignes :

  1. C FATIG3 SOURCE KICH 19/10/25 21:15:15 10351
  2. SUBROUTINE FATIG3(sig,p,raytau,nbrobl,nombre,icri,
  3. &alpha,beta,ycrit,s,ib,igau)
  4. IMPLICIT INTEGER(I-N)
  5. IMPLICIT REAL*8 (A-H,O-Z)
  6. REAL*8 centre(5),stab(3,3),d(3),e(3),v(3),xnorm(3,3),u(5)
  7. REAL*8 sig(nbrobl,nombre),p(nombre),raytau(nombre)
  8. REAL*8 s(nbrobl-1,nombre)
  9.  
  10. SQ2 = DSQRT(2.d0)
  11. SQ6 = DSQRT(6.d0)
  12. C
  13. c *********************************************************
  14. c appel des critere
  15. c *********************************************************
  16. c
  17. if (icri.eq.3) then
  18. C PAPADOPOULOS
  19. c
  20. * call dvpoly(s,centre,elmin,nbrobl,nombre,ib,igau)
  21. call dvboul(s,centre,ray,nbrobl,nombre,iel,igau)
  22. c
  23. pmax = -1.0d10
  24. pmin = 1.0d10
  25. c
  26. stab(1,2) = 0.d0
  27. stab(2,1) = stab(1,2)
  28. stab(1,3) = 0.d0
  29. stab(3,1) = stab(1,3)
  30. do i=1,nombre
  31. if (p(i).gt.pmax) pmax = p(i)
  32. if (p(i).lt.pmin) pmin = p(i)
  33. stab(1,1) = (s(1,i) - centre(1))/SQ2
  34. . +(s(2,i) - centre(2))/SQ6
  35. stab(2,2) =-(s(1,i) - centre(1))/SQ2
  36. . +(s(2,i) - centre(2))/SQ6
  37. stab(3,3) = - stab(1,1) - stab(2,2)
  38. if(nbrobl.gt.4) then
  39. stab(1,2) = s(5,i) - centre(5)
  40. stab(2,1) = s(5,i) - centre(5)
  41. stab(1,3) = s(4,i) - centre(4)
  42. stab(3,1) = s(4,i) - centre(4)
  43. endif
  44. stab(2,3) = s(3,i) - centre(3)
  45. stab(3,2) = s(3,i) - centre(3)
  46. raytau(i) = dsqrt(stab(1,1)**2+stab(2,2)**2+stab(3,3)**2+
  47. . 2.*(stab(1,2)**2+stab(2,3)**2+stab(1,3)**2))/1.4142135
  48. enddo
  49. ycrit = (ray + alpha*pmax - beta)/(beta-alpha*pmax)
  50. c
  51. return
  52. endif
  53. c
  54. if(icri.eq.2) then
  55. C DANG VAN
  56. c
  57. * pmin = 1.d10
  58. * pmax = -1.d10
  59. critere = -1.d10
  60. c
  61. call dvboul(s,centre,ray,nbrobl,nombre,iel,igau)
  62. * if (ib.eq.1.and.igau.eq.1) then
  63. * write(6,*) 'ct',((centre(jj)/1.e6),jj= 1,5)
  64. * write(6,*) (ray/1.e6),nbrobl,nombre
  65. * endif
  66. c
  67. stab(1,2) = 0.d0
  68. stab(2,1) = stab(1,2)
  69. stab(1,3) = 0.d0
  70. stab(3,1) = stab(1,3)
  71. do i=1,nombre
  72. stab(1,1) = (s(1,i) - centre(1))/SQ2
  73. . +(s(2,i) - centre(2))/SQ6
  74. stab(2,2) =-(s(1,i) - centre(1))/SQ2
  75. . +(s(2,i) - centre(2))/SQ6
  76. stab(3,3) = - stab(1,1) - stab(2,2)
  77. if(nbrobl.gt.4) then
  78. stab(1,2) = s(5,i) - centre(5)
  79. stab(2,1) = s(5,i) - centre(5)
  80. stab(1,3) = s(4,i) - centre(4)
  81. stab(3,1) = s(4,i) - centre(4)
  82. endif
  83. stab(2,3) = s(3,i) - centre(3)
  84. stab(3,2) = s(3,i) - centre(3)
  85. call dvdiag(stab,d)
  86. c
  87. c Classement des valeurs propres
  88. c
  89. d3 = dabs(d(1)-d(2))
  90. d2 = dabs(d(1)-d(3))
  91. d1 = dabs(d(3)-d(2))
  92. if (d1.lt.1.d-04) d(2) = d(3)
  93. if (d2.lt.1.d-04) d(1) = d(3)
  94. if (d3.lt.1.d-04) d(2) = d(1)
  95. c
  96. n1 = 1
  97. n2 = 2
  98. n3 = 3
  99. if (d(2).gt.d(1).and.d(3).gt.d(2)) then
  100. n1 = 3
  101. n2 = 2
  102. n3 = 1
  103. endif
  104. if (d(2).gt.d(3).and.d(3).gt.d(1)) then
  105. n1 = 2
  106. n2 = 3
  107. n3 = 1
  108. endif
  109. if (d(3).gt.d(1).and.d(1).gt.d(2)) then
  110. n1 = 3
  111. n2 = 1
  112. n3 = 2
  113. endif
  114. if (d(1).gt.d(3).and.d(3).gt.d(2)) then
  115. n1 = 1
  116. n2 = 3
  117. n3 = 2
  118. endif
  119. if (d(2).gt.d(1).and.d(1).gt.d(3)) then
  120. n1 = 2
  121. n2 = 1
  122. n3 = 3
  123. endif
  124. c
  125. raytau(i) = dabs(d(n1)-d(n3))/2.
  126. crit = raytau(i) + alpha*p(i) - beta
  127. crit = crit / beta
  128. c
  129. * if(p(i).gt.pmax) pmax = p(i)
  130. * if(p(i).lt.pmin) pmin = p(i)
  131. c
  132. if (crit.gt.critere) then
  133. critere = crit
  134. * critere = crit / (beta-alpha*p(i))
  135. v(1) = d(n1)
  136. v(2) = d(n2)
  137. v(3) = d(n3)
  138. xnorm(1,1) = stab(1,n1)
  139. xnorm(2,1) = stab(2,n1)
  140. xnorm(3,1) = stab(3,n1)
  141. xnorm(1,2) = stab(1,n2)
  142. xnorm(2,2) = stab(2,n2)
  143. xnorm(3,2) = stab(3,n2)
  144. xnorm(1,3) = stab(1,n3)
  145. xnorm(2,3) = stab(2,n3)
  146. xnorm(3,3) = stab(3,n3)
  147. endif
  148. enddo
  149. c
  150. ycrit = critere
  151. c
  152. return
  153. endif
  154.  
  155. if(icri.eq.4) then
  156. C SINES
  157. c
  158. dmax = 0.d0
  159. do i=1,nombre-1
  160. do j=i+1,nombre
  161. dist = (s(1,i)-s(1,j))**2+(s(2,i)-s(2,j))**2
  162. . +2.d0*((s(3,i)-s(3,j))**2)
  163. if (nbrobl.gt.4) then
  164. dist = dist + 2.d0*((s(4,i)-s(4,j))**2
  165. . +(s(5,i)-s(5,j))**2)
  166. endif
  167. dist = sqrt(dist)
  168. if (dist.gt.dmax) then
  169. dmax=dist
  170. endif
  171. enddo
  172. enddo
  173. dmax = dmax/2.d0/SQ2
  174. c
  175. pmoyenne = 0.d0
  176. do i=1,nombre
  177. pmoyenne = pmoyenne+p(i)
  178. enddo
  179. pmoyenne = pmoyenne/float(nombre)
  180. ycrit = (dmax+alpha*pmoyenne-beta)/(beta-alpha*pmoyenne)
  181. return
  182. endif
  183.  
  184. if(icri.eq.5) then
  185. C CROSS
  186. c
  187. dmax = 0.d0
  188. do i=1,nombre-1
  189. do j=i+1,nombre
  190. dist = (s(1,i)-s(1,j))**2+(s(2,i)-s(2,j))**2
  191. . +2.d0*((s(3,i)-s(3,j))**2)
  192. if (nbrobl.gt.4) then
  193. dist = dist + 2.d0*((s(4,i)-s(4,j))**2
  194. . +(s(5,i)-s(5,j))**2)
  195. endif
  196. dist = dsqrt(dist)
  197. if (dist.gt.dmax) then
  198. dmax=dist
  199. endif
  200. enddo
  201. enddo
  202. dmax = dmax/2.d0/SQ2
  203. c
  204. pmax = -1.0d10
  205. do i=1,nombre
  206. if (p(i).gt.pmax) pmax=p(i)
  207. enddo
  208. ycrit = (dmax + alpha*pmax - beta)/(beta - alpha*pmax)
  209. return
  210. endif
  211.  
  212. if(icri.eq.6) then
  213. C DC
  214. c
  215. dmax = 0.d0
  216. dmaxp = 0.d0
  217. c
  218. c ***********************
  219. c plus grand diametre ...
  220. c ***********************
  221. do i=1,nombre-1
  222. do j=i+1,nombre
  223. dist = (s(1,i)-s(1,j))**2+(s(2,i)-s(2,j))**2
  224. . +2.d0*((s(3,i)-s(3,j))**2)
  225. if (nbrobl.gt.4) then
  226. dist = dist + 2.d0*((s(4,i)-s(4,j))**2
  227. . +(s(5,i)-s(5,j))**2)
  228. endif
  229. dist = sqrt(dist)
  230. if (dist.gt.dmax) then
  231. dmax=dist
  232. ic = i
  233. jc = j
  234. endif
  235. enddo
  236. enddo
  237. if(dmax.eq.0.) then
  238. dmaxp = 0.
  239. else
  240. c ***********************
  241. c deuxieme diametre ...
  242. c ***********************
  243. do k=1,nbrobl-1
  244. u(k) = (s(k,ic)-s(k,jc))/dmax
  245. enddo
  246. do k=1,nombre
  247. prod = u(1)*s(1,k)+u(2)*s(2,k)
  248. . +2.d0*u(3)*s(3,k)
  249. if (nbrobl.gt.4) then
  250. prod = prod + 2.d0*(u(4)*s(4,k)+u(5)*s(5,k))
  251. endif
  252. do kj=1,nbrobl-1
  253. s(kj,k) = s(kj,k)-prod*u(kj)
  254. enddo
  255. enddo
  256. do i=1,nombre-1
  257. do j=i+1,nombre
  258. dist = (s(1,i)-s(1,j))**2+(s(2,i)-s(2,j))**2
  259. . +2.*((s(3,i)-s(3,j))**2)
  260. if (nbrobl.gt.4) then
  261. dist = dist + 2.d0*((s(4,i)-s(4,j))**2
  262. . +(s(5,i)-s(5,j))**2)
  263. endif
  264. dist = sqrt(dist)
  265. if (dist.gt.dmaxp) dmaxp = dist
  266. enddo
  267. enddo
  268. endif
  269. c
  270. pmax = -1.0d10
  271. do i=1,nombre
  272. if (p(i).gt.pmax) pmax=p(i)
  273. enddo
  274. a = sqrt((dmax/2.)**2+(dmaxp/2.)**2)/SQ2
  275. ycrit = (a + alpha*pmax - beta)/(beta - alpha*pmax)
  276. return
  277. endif
  278.  
  279. if(ICRI.eq.7) then
  280. C VON MISES
  281. c
  282. vonmax = 0.d0
  283. do i=1,nombre
  284. if (nbrobl.gt.4) then
  285. von = sqrt(s(1,i)**2+s(2,i)**2
  286. . +2.*(s(3,i)**2+s(4,i)**2+s(5,i)**2))/SQ2
  287. else
  288. von = sqrt(s(1,i)**2+s(2,i)**2
  289. . +2.*(s(3,i)**2))/SQ2
  290. endif
  291. if (von.gt.vonmax) then
  292. vonmax = von
  293. endif
  294. enddo
  295. c
  296. ycrit = vonmax
  297. return
  298. endif
  299.  
  300. c
  301. 900 format(9x,'1 = dvk',/9x,'2 = papado',/9X,'3 = cross',
  302. . /9X,'4 = sines',/9X,'5 = dc',/9X,'6 = vm'/)
  303. c
  304. c
  305. c
  306. end
  307.  
  308.  
  309.  
  310.  

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