Télécharger act3.eso

Retour à la liste

Numérotation des lignes :

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

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