Télécharger act3.eso

Retour à la liste

Numérotation des lignes :

  1. C ACT3 SOURCE PV 18/04/03 21:15:01 9791
  2. SUBROUTINE ACT3
  3. C--------------------------------------------------------------------
  4. C ACCELERATION POUR DES CHAMPOINTS
  5. C
  6. C APPELE DANS INCREME
  7. C
  8. C LIT DANS L'ORDRE 3 CHAMPS PRIMAUX
  9. C ET 4 CHAMPS DE DUAUX
  10. C
  11. C REND UN CHAMP PRIMAL EXTRAPOLE QUI MINIMISE LE CHAMP DUAL.
  12. C METHODE : ON ESTIME L'APPLICATION TANGENTE ET ON L'UTILISE
  13. C POUR RESOUDRE LE PROBLEME
  14. C
  15. C--------------------------------------------------------------------
  16. C
  17. C AOUT 97 : RAJOUTE L'ACCELERATION AVEC 2 OU 3 CHAMPS DE DEPLACEMENTS
  18. C
  19. C CHAT 05/10/2007 : suppression de la limitation à certaines composantes
  20. C
  21. C PM 19/07/2007 : cas à 2 points mettait IDCH0 à 0 au lieu de IDCH1
  22. C
  23. C--------------------------------------------------------------------
  24. IMPLICIT INTEGER(I-N)
  25. IMPLICIT REAL*8(A-H,O-Z)
  26. -INC CCREEL
  27. -INC SMCHPOI
  28. -INC SMELEME
  29. -INC SMLMOTS
  30. -INC CCOPTIO
  31. xref=1d-6
  32. C
  33. C FABRICATION DES LISTMOTS
  34. C
  35. JGN=4
  36. JGM=100
  37. SEGINI MLMOTS
  38. JGM=0
  39.  
  40. MLMOTF=MLMOTS
  41. *
  42. CALL LIROBJ('CHPOINT',ICHD2,1,IRETOU)
  43. CALL LIROBJ('CHPOINT',ICHD1,1,IRETOU)
  44. CALL LIROBJ('CHPOINT',ICHD0,1,IRETOU)
  45. CALL LIROBJ('CHPOINT',IFO3 ,1,IRETOU)
  46. CALL LIROBJ('CHPOINT',IFO2 ,1,IRETOU)
  47. CALL LIROBJ('CHPOINT',IFO1 ,1,IRETOU)
  48. CALL LIROBJ('CHPOINT',IFO0 ,1,IRETOU)
  49. IF (IERR.NE.0) RETURN
  50.  
  51. * Acceleration avec 3 points
  52. IF (ICHD2.EQ.ICHD1) THEN
  53. write(6,*) ' acceleration sur 3 points '
  54. ICHD2 = 0
  55. IFO3 = 0
  56. ENDIF
  57.  
  58. * Acceleration avec 2 points
  59. IF (ICHD1.EQ.ICHD0) THEN
  60. write(6,*) ' acceleration sur 2 points '
  61. ICHD2 = 0
  62. *PM ICHD0 = 0
  63. ICHD1 = 0
  64. IFO3 = 0
  65. IFO2 = 0
  66. ENDIF
  67.  
  68. *
  69. * Recherche des noms de composante des forces
  70. *
  71. MCHPOI=IFO3
  72. IPA=1
  73. 1 CONTINUE
  74. IF (MCHPOI.NE.0) THEN
  75. SEGACT MCHPOI
  76. DO IA=1,IPCHP(/1)
  77. MSOUPO=IPCHP(IA)
  78. SEGACT MSOUPO
  79. DO IB=1,NOCOMP(/2)
  80. DO IC=1,JGM
  81. IF (MOTS(IC).EQ.NOCOMP(IB)) GO TO 2
  82. ENDDO
  83. ** write (6,*) ' act3 ncomp ',nocomp(ib)
  84. if (nocomp(ib).eq.'FLX ') goto 2
  85.  
  86.  
  87. * on ajoute la nouvelle composante à la liste
  88. JGM=JGM+1
  89. IF (JGM.GT.MOTS(/2)) THEN
  90. SEGADJ MLMOTS
  91. ENDIF
  92. MOTS(JGM)=NOCOMP(IB)
  93. 2 CONTINUE
  94. ENDDO
  95. ENDDO
  96. ENDIF
  97. * On fait de même pour les trois autres champs en entrée
  98. IF (IPA.EQ.1) THEN
  99. IPA=2
  100. MCHPOI=IFO2
  101. GO TO 1
  102. ELSEIF(IPA.EQ.2) THEN
  103. IPA=3
  104. MCHPOI=IFO1
  105. GO TO 1
  106. ELSEIF(IPA.EQ.3) THEN
  107. IPA=4
  108. MCHPOI=IFO0
  109. GO TO 1
  110. ENDIF
  111. SEGADJ MLMOTS
  112. *
  113. * On utilise un procédé d'orthogonalisation de Schmidt
  114. * pour obtenir une base de l'espace engendrée par les quatre forces
  115. * puis on effectue la projection (de coordonnées XL, XM et XN) du
  116. * champ de force nulle sur cet espace
  117. *
  118. IF (IFO3.NE.0) THEN
  119.  
  120. * CAS A 4 POINTS :
  121. * -----------------
  122. * famille des 3 vecteurs variations du champ dual
  123. CALL ADCHPO(IFO1,IFO0,IFO10,1D0,-1D0)
  124. CALL ADCHPO(IFO2,IFO0,IFO20,1D0,-1D0)
  125. CALL ADCHPO(IFO3,IFO0,IFO30,1D0,-1D0)
  126. *
  127. CALL XTX1(IFO10,F12)
  128. CALL XTX1(IFO20,F22)
  129. CALL XTX1(IFO30,F32)
  130. if(f12.LT.XPETIT*1d2) f12=XPETIT*1d2
  131. if(f22.LT.XPETIT*1d2) f22=XPETIT*1d2
  132. if(f32.LT.XPETIT*1d2) f32=XPETIT*1d2
  133. * write (6,*) ' dans act3 f12/f22 ',f12/f22,f12/f32
  134. *
  135. * Détection de flip-flop ?
  136. * (1 et 3 sont trop proches, malgré 2 différent de 1)
  137. if (f12/f22.gt.1.2d0 .and. f12/f32.gt.0.707d0 .and.
  138. > f12/f32.lt.1.414d0) then
  139. * act3 : plan B
  140. CALL ERREUR(858)
  141. XL = (sqrt(f12)+sqrt(f22))/(2*sqrt(f12))
  142. XM = 0.
  143. XN = 0.
  144. GOTO 1010
  145. endif
  146. *
  147. * Construction d'une famille libre orthogonale
  148. CALL XTY1(IFO10,IFO20,MLMOTF,MLMOTF,F1F2)
  149. CALL XTY1(IFO10,IFO30,MLMOTF,MLMOTF,F1F3)
  150. CALL XTY1(IFO20,IFO30,MLMOTF,MLMOTF,F2F3)
  151. CALL XTY1(IFO0 ,IFO10,MLMOTF,MLMOTF,F0F1)
  152. CALL XTY1(IFO0 ,IFO20,MLMOTF,MLMOTF,F0F2)
  153. CALL XTY1(IFO0 ,IFO30,MLMOTF,MLMOTF,F0F3)
  154. *
  155. * normalisation pour ameliorer la precision
  156. * les iteres successives pouvant diminer pas mal si on converge
  157. * A1 = F12
  158. * B1 = F1F2
  159. * C1 = F1F3
  160. * XJ1 = -F0F1
  161. * D1 = F1F2
  162. * E1 = F22
  163. * F1 = F2F3
  164. * XK1 = -F0F2
  165. * G1 = F1F3
  166. * H1 = F2F3
  167. * XI1 = F32
  168. * XL1 = -F0F3
  169. A1 = 1.d0
  170. B1 = F1F2/f12
  171. C1 = F1F3/f12
  172. XJ1 = -F0F1/f12
  173. D1 = F1F2/f22
  174. E1 = 1.D0
  175. F1 = F2F3/f22
  176. XK1 = -F0F2/f22
  177. G1 = F1F3/f32
  178. H1 = F2F3/f32
  179. XI1 = 1.d0
  180. XL1 = -F0F3/f32
  181.  
  182.  
  183.  
  184. * write (6,*) ' f12 f22 f32 ',f12,f22,f32
  185.  
  186. * Calcul des coordonnées de IFO0 dans la famille libre
  187. * et élimination des vecteurs liés
  188. IF (ABS(XI1) .LT. xpetit) GOTO 2000
  189. *
  190. A2 = A1 - ( G1*C1/XI1)
  191. B2 = B1 - ( H1*C1/XI1)
  192. XJ2 = XJ1 - (XL1*C1/XI1)
  193. D2 = D1 - ( G1*F1/XI1)
  194. E2 = E1 - ( H1*F1/XI1)
  195. XK2 = XK1 - (XL1*F1/XI1)
  196. *
  197. IF (E2 .LE. e1*xref) THEN
  198. * write (6,*) ' act3 e2 e1 ',e2,e1
  199. GOTO 2000
  200. ENDIF
  201. *
  202. A3 = A2 -( D2*B2/E2)
  203. XJ3 = XJ2-(XK2*B2/E2)
  204. *
  205. IF (A3 .LE. a1*xref) THEN
  206. * write (6,*) ' act3 a3 a1 ',a3,a1
  207. GOTO 2000
  208. ENDIF
  209. IF (ABS(A2) .LE. A1*1d-6) THEN
  210. * write (6,*) ' act3 a2 a1 ',a2,a1
  211. GOTO 2000
  212. ENDIF
  213.  
  214. XL=XJ3/A3
  215. *
  216. B4 = B2
  217. XJ4 = XJ2 - (A2*XL)
  218.  
  219.  
  220. *
  221. IF (ABS(B4) .LE. ABS(B1)*1D-6) THEN
  222. * write (6,*) ' act3 b4 b1 ',b4,b1
  223. GOTO 2000
  224. ENDIF
  225. *
  226. XM = XJ4/B4
  227. *
  228. C5 = C1
  229. XJ5 = XJ1 - (A1*XL) - (B1*XM)
  230. *
  231. IF (ABS(C5) .LT. xref) GOTO 2000
  232. *
  233. XN=XJ5/C5
  234. *
  235. ref=sqrt(a1**2+b1**2+c1**2)
  236. * write (6,*) ' act3 ',ref,XI1,A3,b4,c5
  237. * write (6,*) ' act3 a3 e2 xi1 ',a3,e2,xi1
  238. if (a3.lt.a1*xref.or.e2.lt.e1*xref) then
  239. write (6,*) ' act3 operateur secant non positif'
  240. goto 2000
  241. endif
  242. GOTO 1010
  243.  
  244. 2000 CONTINUE
  245. * le systeme est singulier on essaye avec un vecteur de moins
  246. *
  247. * write (6,*) ' act3 ',xi1,e2,e1,a3,a1,b4,b1,c5
  248. IF (ABS(E1) .LT. xref) GOTO 2100
  249. *
  250. A2 = A1 - (D1 *B1/E1)
  251. XJ2 = XJ1 - (XK1*B1/E1)
  252. *
  253. IF (A2 .LE. A1*xref) THEN
  254. * write (6,*) ' act3 a2 a1 ',a2,a1
  255. GOTO 2100
  256. ENDIF
  257. *
  258. XL = XJ2/A2
  259. *
  260. B3 = B1
  261. XJ3 = XJ1 - (A1*XL)
  262. *
  263. IF (ABS(B3) .LT. xref) GOTO 2100
  264. *
  265. XM = XJ3/B3
  266. *
  267. XN = 0
  268.  
  269. if (a2.lt.a1*xref) then
  270. * write (6,*) ' act3 a2 a1 ',a2,e1
  271. write (6,*) ' act3 operateur secant non positif'
  272. goto 2100
  273. endif
  274. * act3 : reduction a %i1 dimensions
  275. INTERR(1)=2
  276. CALL ERREUR(859)
  277. GOTO 1010
  278.  
  279. 2100 CONTINUE
  280. * le systeme est singulier on essaye avec un vecteur de moins
  281. *
  282. IF (ABS(A1) .LT. xref ) GOTO 1000
  283. *
  284. XL = XJ1/A1
  285. *
  286. XM = 0
  287. XN = 0
  288.  
  289. * act3 : reduction a %i1 dimensions
  290. * write (6,*) ' act3 a1 ',a1
  291. INTERR(1)=1
  292.  
  293. CALL ERREUR(859)
  294. GOTO 1010
  295.  
  296. 1000 CONTINUE
  297. * act3 : pas d accélération
  298. call erreur(860)
  299. XL = 0
  300. XM = 0
  301. XN = 0
  302.  
  303. 1010 CONTINUE
  304. *
  305. * WRITE (6,*) ' f12*f22/(f1f2*2) ',f12*f22/(f1f2**2)
  306. * WRITE (6,*) ' f12*f32/(f1f3*2) ',f12*f32/(f1f3**2)
  307. * WRITE (6,*) ' f22*f32/(f2f3*2) ',f22*f32/(f2f3**2)
  308. * write (6,*) ' xi1,e2,a3,b4,c5 ',xi1,e2,a3,b4,c5
  309. * write (6,*) ' coefficients acc ',-xl,-xm,-xn
  310. * xl=sign(min(abs(xl),1d1),xl)
  311. * xm=sign(min(abs(xm),1d1),xm)
  312. * xn=sign(min(abs(xn),1d1),xn)
  313. * write (6,*) ' coefficients cor ',-xl,-xm,-xn
  314.  
  315. * Vérification qu'on ne retombe pas sur l'un des duals entrés
  316. * Sinon, on réessaye en tenant compte de moins de vecteurs
  317. if (abs(xl-1).lt.1e-1 .and. abs(xm).lt.1e-1
  318. > .and. abs(xn).lt.1e-1) goto 1000
  319. if (abs(xm-1).lt.1e-1 .and. abs(xn).lt.1e-1
  320. > .and. abs(xl).lt.1e-1) goto 2100
  321. if (abs(xn-1).lt.1e-1 .and. abs(xl).lt.1e-1
  322. > .and. abs(xm).lt.1e-1) goto 2000
  323. *
  324. CALL ADCHPO(ICHD0,ICHD1,ITEMP4,XL+XM+XN,XM+XN)
  325. CALL ADCHPO(ICHD2,ITEMP4,ITEMP5,XN,1D0)
  326. *
  327. * ajustement de la nature du champ
  328. * on prend la meme nature que ICHDO
  329. MCHPOI = ICHD0
  330. MCHPO1 = ITEMP5
  331. SEGACT MCHPOI,MCHPO1
  332. NAT = JATTRI(/1)
  333. NSOUPO = MCHPO1.IPCHP(/1)
  334. SEGADJ MCHPO1
  335. DO 1020 I=1,NAT
  336. MCHPO1.JATTRI(I)=JATTRI(I)
  337. 1020 CONTINUE
  338. SEGDES MCHPOI,MCHPO1
  339. *
  340. * Affectation
  341. *
  342. CALL ECROBJ('CHPOINT',ITEMP5)
  343. *
  344. * Quelques destructions
  345. *
  346. MLMOTS=MLMOTF
  347. SEGSUP MLMOTS
  348. CALL DTCHPO(IFO10)
  349. CALL DTCHPO(IFO20)
  350. CALL DTCHPO(IFO30)
  351. CALL DTCHPO(ITEMP4)
  352. RETURN
  353. ENDIF
  354.  
  355. * CAS A 3 POINTS :
  356. * -----------------
  357. * famille des 2 vecteurs variations du champ dual
  358. IF (IFO2.NE.0.AND.IFO3.EQ.0) THEN
  359. CALL ADCHPO(IFO1,IFO0,IFO10,1D0,-1D0)
  360. CALL ADCHPO(IFO2,IFO0,IFO20,1D0,-1D0)
  361. *
  362. CALL XTX1(IFO10,F12)
  363. CALL XTX1(IFO20,F22)
  364. *
  365. * Construction d'une famille libre orthogonale
  366. CALL XTY1(IFO10,IFO20,MLMOTF,MLMOTF,F1F2)
  367. CALL XTY1(IFO0 ,IFO10,MLMOTF,MLMOTF,F0F1)
  368. CALL XTY1(IFO0 ,IFO20,MLMOTF,MLMOTF,F0F2)
  369. *
  370. * A1 = F12
  371. * B1 = F1F2
  372. * XJ1 = -F0F1
  373. * D1 = F1F2
  374. * E1 = F22
  375. * XK1 = -F0F2
  376. A1 = 1.d0
  377. B1 = F1F2/F12
  378. XJ1 = -F0F1/F12
  379. D1 = F1F2/F22
  380. E1 = 1.D0
  381. XK1 = -F0F2/F22
  382. *
  383. * Calcul des coordonnées de IFO0 dans la famille libre
  384. * et élimination des vecteurs liés
  385. IF (ABS(E1) .LT. xref ) GOTO 2200
  386. *
  387. A2 = A1 - (D1 *B1/E1)
  388. XJ2 = XJ1 - (XK1*B1/E1)
  389. *
  390. IF (A2 .LE. A1*xref) THEN
  391. * write (6,*) ' act3 a2 a1 ',a2,a1
  392. GOTO 2200
  393. ENDIF
  394. *
  395. XL = XJ2/A2
  396. *
  397. B3 = B1
  398. XJ3 = XJ1 - (A1*XL)
  399. *
  400. IF (ABS(B3) .LT. xref) GOTO 2200
  401. *
  402. XM = XJ3/B3
  403. *
  404. XN = 0
  405. *
  406. * write (6,*) ' act3 xl xm ',xl,xm
  407. GOTO 1210
  408.  
  409. 2200 CONTINUE
  410. * le systeme est singulier on essaye avec un vecteur de moins
  411. *
  412. IF (A1 .LT. xref ) GOTO 2300
  413. *
  414. XL = XJ1/A1
  415. *
  416. XM = 0
  417. XN = 0
  418. *
  419. * act3 : reduction a %i1 dimensions
  420. INTERR(1)=1
  421. * CALL ERREUR(859)
  422. GOTO 1210
  423.  
  424. 2300 CONTINUE
  425. * act3 : pas d accélération
  426. call erreur(860)
  427. XL = 0
  428. XM = 0
  429. XN = 0
  430.  
  431. 1210 CONTINUE
  432. * Vérification qu'on ne retombe pas sur l'un des duals entrés
  433. * Sinon, on réessaye en tenant compte de moins de vecteurs
  434. if (abs(xl-1).lt.1e-1.and.abs(xm).lt.1e-1.and.
  435. > abs(xn).lt.1e-1) goto 2300
  436. if (abs(xm-1).lt.1e-1.and.abs(xn).lt.1e-1.and.
  437. > abs(xl).lt.1e-1) goto 2200
  438. * if (abs(xn-1).lt.1e-1.and.abs(xl).lt.1e-1.and.
  439. * > abs(xm).lt.1e-1) goto 2300
  440. *
  441. CALL ADCHPO(ICHD0,ICHD1,ITEMP5,XL+XM,XM)
  442. *
  443. * ajustement de la nature du champ
  444. * on prend la meme nature que ICHDO
  445. *
  446. MCHPOI = ICHD0
  447. MCHPO1 = ITEMP5
  448. SEGACT MCHPOI,MCHPO1
  449. NAT = JATTRI(/1)
  450. NSOUPO = MCHPO1.IPCHP(/1)
  451. SEGADJ MCHPO1
  452. DO 1220 I=1,NAT
  453. MCHPO1.JATTRI(I)=JATTRI(I)
  454. 1220 CONTINUE
  455. SEGDES MCHPOI,MCHPO1
  456. *
  457. * Affectation
  458. *
  459. CALL ECROBJ('CHPOINT',ITEMP5)
  460. *
  461. * Quelques destructions
  462. *
  463. MLMOTS=MLMOTF
  464. SEGSUP MLMOTS
  465. CALL DTCHPO(IFO10)
  466. CALL DTCHPO(IFO20)
  467. RETURN
  468. ENDIF
  469.  
  470. * CAS A 2 POINTS :
  471. * -----------------
  472. IF (IFO2.EQ.0) THEN
  473. CALL ADCHPO(IFO1,IFO0,IFO10,1D0,-1D0)
  474. *
  475. CALL XTX1(IFO10,F12)
  476. *
  477. * Un seul vecteur = famille libre d'emblée
  478. CALL XTY1(IFO0,IFO10,MLMOTF,MLMOTF,F0F1)
  479. *
  480. * A1 = F12
  481. * XJ1 = -F0F1
  482. A1 = 1.D0
  483. XJ1 = -F0F1/F12
  484. *
  485. IF (ABS(A1) .LT. xref ) GOTO 1100
  486. *
  487. XL = XJ1/A1
  488. *
  489. XM = 0
  490. XN = 0
  491. GOTO 1110
  492.  
  493. 1100 CONTINUE
  494. * act3 : pas d accélération
  495. call erreur(860)
  496. XL = 0
  497. XM = 0
  498. XN = 0
  499.  
  500. 1110 CONTINUE
  501. if (abs(xl-1).lt.1e-1) goto 1100
  502. CALL MUCHPO(ICHD0,XL,ITEMP5,1)
  503. * write (6,*) ' act3 xl vaut ',xl
  504. *
  505. * ajustement de la nature du champ
  506. * on prend la meme nature que ICHDO
  507. *
  508. MCHPOI = ICHD0
  509. MCHPO1 = ITEMP5
  510. SEGACT MCHPOI,MCHPO1
  511. NAT = JATTRI(/1)
  512. NSOUPO = MCHPO1.IPCHP(/1)
  513. SEGADJ MCHPO1
  514. DO 1120 I=1,NAT
  515. MCHPO1.JATTRI(I)=JATTRI(I)
  516. 1120 CONTINUE
  517. SEGDES MCHPOI,MCHPO1
  518. *
  519. * Affectation
  520. *
  521. CALL ECROBJ('CHPOINT',ITEMP5)
  522. *
  523. * QUELQUES DESTRUCTIONS
  524. *
  525. MLMOTS=MLMOTF
  526. SEGSUP MLMOTS
  527. CALL DTCHPO(IFO10)
  528. RETURN
  529. ENDIF
  530. END
  531.  
  532.  
  533.  
  534.  
  535.  
  536.  
  537.  
  538.  
  539.  

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