Télécharger act3.eso

Retour à la liste

Numérotation des lignes :

  1. C ACT3 SOURCE PV 19/01/31 21:15:00 10091
  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. if (xl .gt.1.D0) then
  216. ** write(6,*) ' act3 xl > 0 ==> pas d acceleration'
  217. goto 2000
  218. endif
  219. *
  220. B4 = B2
  221. XJ4 = XJ2 - (A2*XL)
  222.  
  223.  
  224. *
  225. IF (ABS(B4) .LE. ABS(B1)*1D-6) THEN
  226. * write (6,*) ' act3 b4 b1 ',b4,b1
  227. GOTO 2000
  228. ENDIF
  229. *
  230. XM = XJ4/B4
  231. *
  232. C5 = C1
  233. XJ5 = XJ1 - (A1*XL) - (B1*XM)
  234. *
  235. IF (ABS(C5) .LT. xref) GOTO 2000
  236. *
  237. XN=XJ5/C5
  238. *
  239. ref=sqrt(a1**2+b1**2+c1**2)
  240. * write (6,*) ' act3 ',ref,XI1,A3,b4,c5
  241. * write (6,*) ' act3 a3 e2 xi1 ',a3,e2,xi1
  242. if (a3.lt.a1*xref.or.e2.lt.e1*xref) then
  243. write (6,*) ' act3 operateur secant non positif'
  244. goto 2000
  245. endif
  246. GOTO 1010
  247.  
  248. 2000 CONTINUE
  249. * le systeme est singulier on essaye avec un vecteur de moins
  250. *
  251. * write (6,*) ' act3 ',xi1,e2,e1,a3,a1,b4,b1,c5
  252. IF (ABS(E1) .LT. xref) GOTO 2100
  253. *
  254. A2 = A1 - (D1 *B1/E1)
  255. XJ2 = XJ1 - (XK1*B1/E1)
  256. *
  257. IF (A2 .LE. A1*xref) THEN
  258. * write (6,*) ' act3 a2 a1 ',a2,a1
  259. GOTO 2100
  260. ENDIF
  261. *
  262. XL = XJ2/A2
  263. if (xl .gt.1.D0) then
  264. ** write(6,*) ' act3 xl > 0 ==> pas d acceleration'
  265. goto 2100
  266. endif
  267. *
  268. B3 = B1
  269. XJ3 = XJ1 - (A1*XL)
  270. *
  271. IF (ABS(B3) .LT. xref) GOTO 2100
  272. *
  273. XM = XJ3/B3
  274. *
  275. XN = 0
  276.  
  277. if (a2.lt.a1*xref) then
  278. * write (6,*) ' act3 a2 a1 ',a2,e1
  279. write (6,*) ' act3 operateur secant non positif'
  280. goto 2100
  281. endif
  282. * act3 : reduction a %i1 dimensions
  283. INTERR(1)=2
  284. CALL ERREUR(859)
  285. GOTO 1010
  286.  
  287. 2100 CONTINUE
  288. * le systeme est singulier on essaye avec un vecteur de moins
  289. *
  290. IF (ABS(A1) .LT. xref ) GOTO 1000
  291. *
  292. XL = XJ1/A1
  293. if (xl .gt.1.D0) then
  294. ** write(6,*) ' act3 xl > 0 ==> pas d acceleration'
  295. goto 1000
  296. endif
  297. *
  298. XM = 0
  299. XN = 0
  300.  
  301. * act3 : reduction a %i1 dimensions
  302. * write (6,*) ' act3 a1 ',a1
  303. INTERR(1)=1
  304.  
  305. CALL ERREUR(859)
  306. GOTO 1010
  307.  
  308. 1000 CONTINUE
  309. * act3 : pas d accélération
  310. call erreur(860)
  311. XL = 0
  312. XM = 0
  313. XN = 0
  314.  
  315. 1010 CONTINUE
  316. *
  317. * WRITE (6,*) ' f12*f22/(f1f2*2) ',f12*f22/(f1f2**2)
  318. * WRITE (6,*) ' f12*f32/(f1f3*2) ',f12*f32/(f1f3**2)
  319. * WRITE (6,*) ' f22*f32/(f2f3*2) ',f22*f32/(f2f3**2)
  320. * write (6,*) ' xi1,e2,a3,b4,c5 ',xi1,e2,a3,b4,c5
  321. * write (6,*) ' coefficients acc ',-xl,-xm,-xn
  322. * xl=sign(min(abs(xl),1d1),xl)
  323. * xm=sign(min(abs(xm),1d1),xm)
  324. * xn=sign(min(abs(xn),1d1),xn)
  325. * write (6,*) ' coefficients cor ',-xl,-xm,-xn
  326.  
  327. * Vérification qu'on ne retombe pas sur l'un des duals entrés
  328. * Sinon, on réessaye en tenant compte de moins de vecteurs
  329. if (abs(xl-1).lt.1e-1 .and. abs(xm).lt.1e-1
  330. > .and. abs(xn).lt.1e-1) goto 1000
  331. if (abs(xm-1).lt.1e-1 .and. abs(xn).lt.1e-1
  332. > .and. abs(xl).lt.1e-1) goto 2100
  333. if (abs(xn-1).lt.1e-1 .and. abs(xl).lt.1e-1
  334. > .and. abs(xm).lt.1e-1) goto 2000
  335. *
  336. CALL ADCHPO(ICHD0,ICHD1,ITEMP4,XL+XM+XN,XM+XN)
  337. CALL ADCHPO(ICHD2,ITEMP4,ITEMP5,XN,1D0)
  338. *
  339. * ajustement de la nature du champ
  340. * on prend la meme nature que ICHDO
  341. MCHPOI = ICHD0
  342. MCHPO1 = ITEMP5
  343. SEGACT MCHPOI,MCHPO1
  344. NAT = JATTRI(/1)
  345. NSOUPO = MCHPO1.IPCHP(/1)
  346. SEGADJ MCHPO1
  347. DO 1020 I=1,NAT
  348. MCHPO1.JATTRI(I)=JATTRI(I)
  349. 1020 CONTINUE
  350. SEGDES MCHPOI,MCHPO1
  351. *
  352. * Affectation
  353. *
  354. CALL ECROBJ('CHPOINT',ITEMP5)
  355. *
  356. * Quelques destructions
  357. *
  358. MLMOTS=MLMOTF
  359. SEGSUP MLMOTS
  360. CALL DTCHPO(IFO10)
  361. CALL DTCHPO(IFO20)
  362. CALL DTCHPO(IFO30)
  363. CALL DTCHPO(ITEMP4)
  364. RETURN
  365. ENDIF
  366.  
  367. * CAS A 3 POINTS :
  368. * -----------------
  369. * famille des 2 vecteurs variations du champ dual
  370. IF (IFO2.NE.0.AND.IFO3.EQ.0) THEN
  371. CALL ADCHPO(IFO1,IFO0,IFO10,1D0,-1D0)
  372. CALL ADCHPO(IFO2,IFO0,IFO20,1D0,-1D0)
  373. *
  374. CALL XTX1(IFO10,F12)
  375. CALL XTX1(IFO20,F22)
  376. *
  377. * Construction d'une famille libre orthogonale
  378. CALL XTY1(IFO10,IFO20,MLMOTF,MLMOTF,F1F2)
  379. CALL XTY1(IFO0 ,IFO10,MLMOTF,MLMOTF,F0F1)
  380. CALL XTY1(IFO0 ,IFO20,MLMOTF,MLMOTF,F0F2)
  381. *
  382. * A1 = F12
  383. * B1 = F1F2
  384. * XJ1 = -F0F1
  385. * D1 = F1F2
  386. * E1 = F22
  387. * XK1 = -F0F2
  388. A1 = 1.d0
  389. B1 = F1F2/F12
  390. XJ1 = -F0F1/F12
  391. D1 = F1F2/F22
  392. E1 = 1.D0
  393. XK1 = -F0F2/F22
  394. *
  395. * Calcul des coordonnées de IFO0 dans la famille libre
  396. * et élimination des vecteurs liés
  397. IF (ABS(E1) .LT. xref ) GOTO 2200
  398. *
  399. A2 = A1 - (D1 *B1/E1)
  400. XJ2 = XJ1 - (XK1*B1/E1)
  401. *
  402. IF (A2 .LE. A1*xref) THEN
  403. * write (6,*) ' act3 a2 a1 ',a2,a1
  404. GOTO 2200
  405. ENDIF
  406. *
  407. XL = XJ2/A2
  408. *
  409. B3 = B1
  410. XJ3 = XJ1 - (A1*XL)
  411. *
  412. IF (ABS(B3) .LT. xref) GOTO 2200
  413. *
  414. XM = XJ3/B3
  415. *
  416. XN = 0
  417. *
  418. * write (6,*) ' act3 xl xm ',xl,xm
  419. GOTO 1210
  420.  
  421. 2200 CONTINUE
  422. * le systeme est singulier on essaye avec un vecteur de moins
  423. *
  424. IF (A1 .LT. xref ) GOTO 2300
  425. *
  426. XL = XJ1/A1
  427. *
  428. XM = 0
  429. XN = 0
  430. *
  431. * act3 : reduction a %i1 dimensions
  432. INTERR(1)=1
  433. * CALL ERREUR(859)
  434. GOTO 1210
  435.  
  436. 2300 CONTINUE
  437. * act3 : pas d accélération
  438. call erreur(860)
  439. XL = 0
  440. XM = 0
  441. XN = 0
  442.  
  443. 1210 CONTINUE
  444. * Vérification qu'on ne retombe pas sur l'un des duals entrés
  445. * Sinon, on réessaye en tenant compte de moins de vecteurs
  446. if (abs(xl-1).lt.1e-1.and.abs(xm).lt.1e-1.and.
  447. > abs(xn).lt.1e-1) goto 2300
  448. if (abs(xm-1).lt.1e-1.and.abs(xn).lt.1e-1.and.
  449. > abs(xl).lt.1e-1) goto 2200
  450. * if (abs(xn-1).lt.1e-1.and.abs(xl).lt.1e-1.and.
  451. * > abs(xm).lt.1e-1) goto 2300
  452. *
  453. CALL ADCHPO(ICHD0,ICHD1,ITEMP5,XL+XM,XM)
  454. *
  455. * ajustement de la nature du champ
  456. * on prend la meme nature que ICHDO
  457. *
  458. MCHPOI = ICHD0
  459. MCHPO1 = ITEMP5
  460. SEGACT MCHPOI,MCHPO1
  461. NAT = JATTRI(/1)
  462. NSOUPO = MCHPO1.IPCHP(/1)
  463. SEGADJ MCHPO1
  464. DO 1220 I=1,NAT
  465. MCHPO1.JATTRI(I)=JATTRI(I)
  466. 1220 CONTINUE
  467. SEGDES MCHPOI,MCHPO1
  468. *
  469. * Affectation
  470. *
  471. CALL ECROBJ('CHPOINT',ITEMP5)
  472. *
  473. * Quelques destructions
  474. *
  475. MLMOTS=MLMOTF
  476. SEGSUP MLMOTS
  477. CALL DTCHPO(IFO10)
  478. CALL DTCHPO(IFO20)
  479. RETURN
  480. ENDIF
  481.  
  482. * CAS A 2 POINTS :
  483. * -----------------
  484. IF (IFO2.EQ.0) THEN
  485. CALL ADCHPO(IFO1,IFO0,IFO10,1D0,-1D0)
  486. *
  487. CALL XTX1(IFO10,F12)
  488. *
  489. * Un seul vecteur = famille libre d'emblée
  490. CALL XTY1(IFO0,IFO10,MLMOTF,MLMOTF,F0F1)
  491. *
  492. * A1 = F12
  493. * XJ1 = -F0F1
  494. A1 = 1.D0
  495. XJ1 = -F0F1/F12
  496. *
  497. IF (ABS(A1) .LT. xref ) GOTO 1100
  498. *
  499. XL = XJ1/A1
  500. *
  501. XM = 0
  502. XN = 0
  503. GOTO 1110
  504.  
  505. 1100 CONTINUE
  506. * act3 : pas d accélération
  507. call erreur(860)
  508. XL = 0
  509. XM = 0
  510. XN = 0
  511.  
  512. 1110 CONTINUE
  513. if (abs(xl-1).lt.1e-1) goto 1100
  514. CALL MUCHPO(ICHD0,XL,ITEMP5,1)
  515. * write (6,*) ' act3 xl vaut ',xl
  516. *
  517. * ajustement de la nature du champ
  518. * on prend la meme nature que ICHDO
  519. *
  520. MCHPOI = ICHD0
  521. MCHPO1 = ITEMP5
  522. SEGACT MCHPOI,MCHPO1
  523. NAT = JATTRI(/1)
  524. NSOUPO = MCHPO1.IPCHP(/1)
  525. SEGADJ MCHPO1
  526. DO 1120 I=1,NAT
  527. MCHPO1.JATTRI(I)=JATTRI(I)
  528. 1120 CONTINUE
  529. SEGDES MCHPOI,MCHPO1
  530. *
  531. * Affectation
  532. *
  533. CALL ECROBJ('CHPOINT',ITEMP5)
  534. *
  535. * QUELQUES DESTRUCTIONS
  536. *
  537. MLMOTS=MLMOTF
  538. SEGSUP MLMOTS
  539. CALL DTCHPO(IFO10)
  540. RETURN
  541. ENDIF
  542. END
  543.  
  544.  
  545.  
  546.  
  547.  
  548.  
  549.  
  550.  
  551.  
  552.  
  553.  

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