Télécharger act3.eso

Retour à la liste

Numérotation des lignes :

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

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