Télécharger act3.eso

Retour à la liste

Numérotation des lignes :

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

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