Télécharger pre122.eso

Retour à la liste

Numérotation des lignes :

  1. C PRE122 SOURCE PV 09/03/12 21:30:37 6325
  2. C PRE121 SOURCE BECC 98/08/24 21:17:00 3286
  3. SUBROUTINE PRE122(LOGTEM,
  4. & ICEN,IFACE,IFACEL,INORM,
  5. & IROC, IGRROC, IALROC,
  6. & IVITC, IGRVC, IALVC,
  7. & IPC ,IGRPC, IALPC,
  8. & IGAMC,
  9. & DELTAT,
  10. & IROF,IVITF,IPF,IGAMF,
  11. & LOGAN,LOGNEG,LOGBOR,MESERR,VALER,VAL1,VAL2)
  12. C************************************************************************
  13. C
  14. C PROJET : CASTEM 2000
  15. C
  16. C NOM : PRE121
  17. C
  18. C DESCRIPTION : Voir PRE12
  19. C
  20. C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec extensions CISI)
  21. C
  22. C AUTEUR : A. BECCANTINI, DRN/DMT/SEMT/TTMF
  23. C
  24. C************************************************************************
  25. C
  26. C
  27. C APPELES (Outils) : KRIPAD, LICHT
  28. C
  29. C APPELES (Calcul) : AUCUN
  30. C
  31. C
  32. C************************************************************************
  33. C
  34. C ENTREES
  35. C
  36. C LOGTEM : LOGICAL; si .TRUE. 2em ordre en temps
  37. C sinon 1er ordre en temps
  38. C
  39. C 1) Pointeurs de MELEMEs et de CHPOINTs de la table DOMAINE
  40. C
  41. C ICEN : MELEME de 'POI1' SPG des CENTRES
  42. C
  43. C IFACE : MELEME de 'POI1' SPG des FACES
  44. C
  45. C IFACEL : MELEME de 'SEG3' avec
  46. C CENTRE d'Elt "gauche"
  47. C CENTRE de Face
  48. C CENTRE d'Elt "droite"
  49. C
  50. C N.B. = IFACE.NUM(i,1) = IFACEL.NUM(i,2)
  51. C
  52. C INORM : CHPOINT des cosinus directeurs de normales aux faces
  53. C
  54. C 2) Pointeurs des CHPOINTs
  55. C
  56. C IROC : CHPOINT "CENTRE" contenant la masse volumique RHO
  57. C
  58. C IGRROC : CHPOINT "CENTRE" contenant le gradient de la
  59. C masse volumique RHO (2 composantes)
  60. C
  61. C IALROC : CHPOINT "CENTRE" contenant le limiteur du gradient
  62. C de la masse volumique
  63. C
  64. C IVITC : CHPOINT "CENTRE" contenant la vitesse UX, UY ;
  65. C
  66. C IGRVC : CHPOINT "CENTRE" contenant le gradient de la
  67. C vitesse (4 composantes)
  68. C
  69. C IALVC : CHPOINT "CENTRE" contenant le limiteur du gradient
  70. C de la vitesse (2 composantes)
  71. C
  72. C IPC : CHPOINT "CENTRE" contenat la pression P;
  73. C
  74. C IGRPC : CHPOINT "CENTRE" contenant le gradient de la
  75. C pression (2 composantes)
  76. C
  77. C IALPC : CHPOINT "CENTRE" contenant le limiteur du gradient
  78. C de la pression
  79. C
  80. C IGAMC : CHPOINT "CENTRE" contenat le "Gamma" du gaz
  81. C
  82. C 3)
  83. C
  84. C DELTAT : REAL*8, encrement en temps pour calculer la prediction
  85. C
  86. C
  87. C SORTIES
  88. C
  89. C
  90. C IROF : MCHAML defini sur le MELEME de pointeur IFACEL,
  91. C contenant la masse volumique RHO
  92. C
  93. C IVITF : MCHAML defini sur le MELEME de pointeur IFACEL,
  94. C contenant la vitesse UN, UT dans le repaire local
  95. C (n,t) et defini sur le MELEME de pointeur IFACE,
  96. C contenant les cosinus directeurs du repere local
  97. C
  98. C IPF : MCHAML defini sur le MELEME de pointeur IFACEL,
  99. C contenant la pression P
  100. C
  101. C IGAMF : MCHAML defini sur le MELEME de pointeur IFACEL,
  102. C contenant le "gamma" du gaz
  103. C
  104. C LOGAN : anomalie detectee
  105. C
  106. C LOGNEG : (LOGICAL): si .TRUE. une pression ou une densité
  107. C negative a été detectée -> en interactif le
  108. C programme s'arrete en GIBIANE
  109. C (erreur stocké en MESERR et VALER)
  110. C
  111. C LOGBOR : (LOGICAL): si .TRUE. un gamma a ete detecte
  112. C dehor 1 et 3 (sa valeur stockée en MESERR et VALER;
  113. C en VAL1 et en VAL2 on stocke 1.0 et 3.0)
  114. C
  115. C MESERR
  116. C VALER
  117. C VAL1,
  118. C VAL2 : pour les messages d'erreur
  119. C
  120. C************************************************************************
  121. C
  122. C HISTORIQUE (Anomalies et modifications éventuelles)
  123. C
  124. C HISTORIQUE : Créée le 11.6.98.
  125. C
  126. C************************************************************************
  127. C
  128. C
  129. C ATTENTION: Cet programme marche si le MAILLAGE est convex;
  130. C si non il faut changer l'argoritme de calcul de
  131. C l'orientation des normales aux faces.
  132. C
  133. C La positivité n'est pas controlle parce que c'est déjà fait
  134. C dans l'operateur PRIM
  135. C
  136. C
  137. C************************************************************************
  138. C
  139. C**** Variables de COOPTIO
  140. C
  141. C INTEGER IPLLB, IERPER, IERMAX, IERR, INTERR
  142. C & ,IOTER, IOLEC, IOIMP, IOCAR, IOACQ
  143. C & ,IOPER, IOSGB, IOGRA, IOSAU, IORES
  144. C & ,IECHO, IIMPI, IOSPI
  145. C & ,IDIM
  146. CC & ,MCOORD
  147. C & ,IFOMOD, NIFOUR, IFOUR, NSDPGE, IONIVE
  148. C & ,NGMAXY, IZROSF, ISOTYP, IOSCR,LTEXLU
  149. C & ,NORINC,NORVAL,NORIND,NORVAD
  150. C & ,NUCROU, IPSAUV
  151. C
  152. C**** Les variables
  153. C
  154. IMPLICIT INTEGER(I-N)
  155. INTEGER ICEN, IFACE, IFACEL, INORM,IROC, IGRROC, IALROC
  156. & , IVITC, IGRVC, IALVC
  157. & , IPC ,IGRPC, IALPC
  158. & , IGAMC
  159. & , IROF, IVITF, IPF, IGAMF
  160. & , IGEOM, NFAC, NCEN
  161. & , N1PTEL, N1EL, N2PTEL, N2EL, N2, N1, N3, L1, NLCE
  162. & , NLCF, NGCF, NGCEG, NLCEG, NGCED, NLCED, NGCF1
  163. & , IDIMP1, INDCEL
  164. REAL*8 VALER, VAL1, VAL2, XG, YG, ZG, XC, YC, ZC, XD, YD, ZD,
  165. & DELTAT
  166. & ,DXG, DYG, DZG, DXD, DYD, DZD,ORIENT
  167. & , CNX, CNY, CNZ, CTX, CTY, CTZ, CVX, CVY, CVZ
  168. & , ROG, PG, GAMG, UXG, UYG, UZG, UNG, UTG, UVG
  169. & , ROD, PD, GAMD, UXD, UYD, UZD, UND, UTD, UVD
  170. & , VALCEL, DCEL, ALCEL
  171. & , DROX, DROY, DROZ, DUXX, DUXY, DUXZ, DUYX, DUYY,
  172. & DUYZ, DUZX, DUZY, DUZZ, DPX, DPY, DPZ
  173. & , DRO, DUX, DUY, DUZ, DP, ALPHA(5)
  174. CHARACTER*(40) MESERR
  175. CHARACTER*(8) TYPE
  176. LOGICAL LOGAN,LOGNEG, LOGBOR, LOGTEM, LOGI1, LOGALP
  177. C
  178. C**** LOGALP = .TRUE. -> Prediction avec limiteur
  179. C
  180. PARAMETER(LOGALP = .TRUE.)
  181. C
  182. C**** Les Includes
  183. C
  184. -INC SMCOORD
  185. -INC CCOPTIO
  186. -INC SMCHPOI
  187. POINTEUR MPROC.MPOVAL, MPGRR.MPOVAL,
  188. & MPVITC.MPOVAL, MPGRV.MPOVAL,
  189. & MPPC.MPOVAL, MPGRP.MPOVAL,
  190. & MPGAMC.MPOVAL, MPNORM.MPOVAL,
  191. & MPROP.MPOVAL, MPPP.MPOVAL, MPVITP.MPOVAL
  192. -INC SMCHAML
  193. POINTEUR MELVNX.MELVAL, MELVNY.MELVAL, MELVNZ.MELVAL,
  194. & MELT1X.MELVAL, MELT1Y.MELVAL, MELT1Z.MELVAL,
  195. & MELT2X.MELVAL, MELT2Y.MELVAL, MELT2Z.MELVAL
  196. POINTEUR MELVUN.MELVAL, MELVUT.MELVAL, MELVUV.MELVAL
  197. POINTEUR MELRO.MELVAL, MELP.MELVAL,
  198. & MELGAM.MELVAL
  199. POINTEUR MELALR.MPOVAL,
  200. & MELALV.MPOVAL,
  201. & MELALP.MPOVAL
  202. -INC SMLENTI
  203. -INC SMELEME
  204. C
  205. C
  206. C**** Initialisation des parametres d'erreur déjà faite, i.e.
  207. C
  208. C LOGNEG = .FALSE.
  209. C LOGBOR = .FALSE.
  210. C MESERR = ' '
  211. C MOTERR(1:40) = MESERR(1:40)
  212. C VALER = 0.0D0
  213. C VAL1 = 0.0D0
  214. C VAL2 = 0.0D0
  215. C
  216. C
  217. C**** KRIPAD pour la correspondance global/local de centre
  218. C
  219. CALL KRIPAD(ICEN,MLENT1)
  220. C
  221. C**** MLENTI1 a MCORD.XCOORD(/1)/(IDIM+1) elements
  222. C
  223. C Si i est le numero global d'un noeud de ICEN,
  224. C MLENT1.LECT(i) contient sa position, i.e.
  225. C
  226. C I = numero global du noeud centre
  227. C MLENT1.LECT(i) = numero local du noeud centre
  228. C
  229. C MLENT1 déjà activé, i.e.
  230. C
  231. C SEGACT MLENT1
  232. C
  233. C**** Activation de CHPOINTs
  234. C
  235. C densité + grad + limiteur
  236. C vitesse + grad + limiteur
  237. C pression + grad + limiteur
  238. C gamma
  239. C cosinus directeurs des normales aux surface
  240. C
  241. CALL LICHT(IROC , MPROC , TYPE, IGEOM)
  242. CALL LICHT(IGRROC, MPGRR , TYPE, IGEOM)
  243. CALL LICHT(IVITC, MPVITC , TYPE, IGEOM)
  244. CALL LICHT(IGRVC, MPGRV , TYPE, IGEOM)
  245. CALL LICHT(IPC , MPPC , TYPE, IGEOM)
  246. CALL LICHT(IGRPC, MPGRP , TYPE, IGEOM)
  247. CALL LICHT(IGAMC, MPGAMC , TYPE, IGEOM)
  248. CALL LICHT(INORM, MPNORM , TYPE, IGEOM)
  249. C
  250. C**** Les MPOVALs 'Prediction'
  251. C
  252. IF(LOGTEM)THEN
  253. SEGINI, MPROP = MPROC
  254. SEGINI, MPPP = MPPC
  255. SEGINI, MPVITP = MPVITC
  256. ELSE
  257. MPROP = MPROC
  258. MPPP = MPPC
  259. MPVITP = MPVITC
  260. ENDIF
  261. C
  262. C**** Les Limiteurs
  263. C
  264. CALL LICHT(IALROC, MELALR , TYPE, IGEOM)
  265. CALL LICHT(IALVC, MELALV , TYPE, IGEOM)
  266. CALL LICHT(IALPC, MELALP , TYPE, IGEOM)
  267. C
  268. C
  269. C**** MPOVAL sont déjà activés i.e.:
  270. C
  271. C SEGACT MPROC
  272. C SEGACT MPGRR
  273. C SEGACT MPIALR
  274. C SEGACT MPVITC
  275. C SEGACT MPGRV
  276. C SEGACT MPIALV
  277. C SEGACT MPPC
  278. C SEGACT MPGRP
  279. C SEGACT MPIALP
  280. C SEGACT MPGAMC
  281. C SEGACT MPNORM
  282. C
  283. C**** Le MELEME FACEL
  284. C
  285. IPT1 = IFACEL
  286. IPT2 = IFACE
  287. SEGACT IPT1
  288. SEGACT IPT2
  289. NFAC = IPT1.NUM(/2)
  290. C
  291. C**** Creation de MCHAMLs contenant les etats gauche et droite,
  292. C
  293. C i.e.:
  294. C
  295. C vitesse + cosinus directors du repere local
  296. C densité
  297. C pression
  298. C gamma
  299. C
  300. C**** Cosinus directors du repere local et vitesse
  301. C
  302. C Les cosinus directeurs
  303. C
  304. N1 = 2
  305. N3 = 6
  306. L1 = 28
  307. SEGINI MCHEL1
  308. IVITF = MCHEL1
  309. MCHEL1.TITCHE = 'U '
  310. MCHEL1.IMACHE(1) = IFACE
  311. MCHEL1.IMACHE(2) = IFACEL
  312. MCHEL1.CONCHE(1) = '(n,t,v)in(x,y,z)'
  313. MCHEL1.CONCHE(2) = ' U in (n,t,v) '
  314. C
  315. C**** Valeurs des cosinus definies par respect au repair global, i.e.
  316. C
  317. MCHEL1.INFCHE(1,1) = 2
  318. MCHEL1.INFCHE(1,3) = NIFOUR
  319. MCHEL1.INFCHE(1,4) = 0
  320. MCHEL1.INFCHE(1,5) = 0
  321. MCHEL1.INFCHE(1,6) = 0
  322. MCHEL1.IFOCHE = IFOUR
  323. C
  324. C**** Valeurs de vitesse definies par respect au repair local, i.e.
  325. C
  326. MCHEL1.INFCHE(2,1) = 1
  327. MCHEL1.INFCHE(2,3) = NIFOUR
  328. MCHEL1.INFCHE(2,4) = 0
  329. MCHEL1.INFCHE(2,5) = 0
  330. MCHEL1.INFCHE(2,6) = 0
  331. C
  332. C**** Le cosinus directeurs
  333. C
  334. N1PTEL = 1
  335. N1EL = NFAC
  336. N2PTEL = 0
  337. N2EL = 0
  338. C
  339. C**** MCHAML a N2 composantes:
  340. C
  341. C cosinus directeurs du repere local (n,t1,t2)
  342. C
  343. C IDIM = 3 -> 9 composantes
  344. C
  345. N2 = 9
  346. SEGINI MCHAM1
  347. MCHEL1.ICHAML(1) = MCHAM1
  348. MCHAM1.NOMCHE(1) = 'NX '
  349. MCHAM1.NOMCHE(2) = 'NY '
  350. MCHAM1.NOMCHE(3) = 'NZ '
  351. MCHAM1.NOMCHE(4) = 'TX '
  352. MCHAM1.NOMCHE(5) = 'TY '
  353. MCHAM1.NOMCHE(6) = 'TZ '
  354. MCHAM1.NOMCHE(7) = 'VX '
  355. MCHAM1.NOMCHE(8) = 'VY '
  356. MCHAM1.NOMCHE(9) = 'VZ '
  357. MCHAM1.TYPCHE(1) = 'REAL*8 '
  358. MCHAM1.TYPCHE(2) = 'REAL*8 '
  359. MCHAM1.TYPCHE(3) = 'REAL*8 '
  360. MCHAM1.TYPCHE(4) = 'REAL*8 '
  361. MCHAM1.TYPCHE(5) = 'REAL*8 '
  362. MCHAM1.TYPCHE(6) = 'REAL*8 '
  363. MCHAM1.TYPCHE(7) = 'REAL*8 '
  364. MCHAM1.TYPCHE(8) = 'REAL*8 '
  365. MCHAM1.TYPCHE(9) = 'REAL*8 '
  366. SEGINI MELVNX
  367. SEGINI MELVNY
  368. SEGINI MELVNZ
  369. SEGINI MELT1X
  370. SEGINI MELT1Y
  371. SEGINI MELT1Z
  372. SEGINI MELT2X
  373. SEGINI MELT2Y
  374. SEGINI MELT2Z
  375. MCHAM1.IELVAL(1) = MELVNX
  376. MCHAM1.IELVAL(2) = MELVNY
  377. MCHAM1.IELVAL(3) = MELVNZ
  378. MCHAM1.IELVAL(4) = MELT1X
  379. MCHAM1.IELVAL(5) = MELT1Y
  380. MCHAM1.IELVAL(6) = MELT1Z
  381. MCHAM1.IELVAL(7) = MELT2X
  382. MCHAM1.IELVAL(8) = MELT2Y
  383. MCHAM1.IELVAL(9) = MELT2Z
  384. SEGDES MCHAM1
  385. C
  386. C**** Vitesse
  387. C
  388. N1EL = NFAC
  389. N1PTEL = 3
  390. N2EL = 0
  391. N2PTEL = 0
  392. C
  393. C**** MCHAML a N2 composantes:
  394. C
  395. C IDIM = 3 -> 3 composantes
  396. C
  397. N2 = 3
  398. SEGINI MCHAM1
  399. MCHEL1.ICHAML(2) = MCHAM1
  400. SEGDES MCHEL1
  401. MCHAM1.NOMCHE(1) = 'UN '
  402. MCHAM1.NOMCHE(2) = 'UT '
  403. MCHAM1.NOMCHE(3) = 'UV '
  404. MCHAM1.TYPCHE(1) = 'REAL*8 '
  405. MCHAM1.TYPCHE(2) = 'REAL*8 '
  406. MCHAM1.TYPCHE(3) = 'REAL*8 '
  407. SEGINI MELVUN
  408. SEGINI MELVUT
  409. SEGINI MELVUV
  410. MCHAM1.IELVAL(1) = MELVUN
  411. MCHAM1.IELVAL(2) = MELVUT
  412. MCHAM1.IELVAL(3) = MELVUV
  413. SEGDES MCHAM1
  414. C
  415. C**** Densite
  416. C
  417. N1 = 1
  418. N3 = 6
  419. L1 = 15
  420. SEGINI MCHEL2
  421. IROF = MCHEL2
  422. MCHEL2.IMACHE(1) = IFACEL
  423. MCHEL2.TITCHE = 'RO '
  424. MCHEL2.CONCHE(1) = ' '
  425. C
  426. C**** Valeurs independente du repére, i.e.
  427. C
  428. MCHEL2.INFCHE(1,1) = 0
  429. MCHEL2.INFCHE(1,3) = NIFOUR
  430. MCHEL2.INFCHE(1,4) = 0
  431. MCHEL2.INFCHE(1,5) = 0
  432. MCHEL2.INFCHE(1,6) = 0
  433. MCHEL2.IFOCHE = IFOUR
  434. N2 = 1
  435. SEGINI MCHAM1
  436. MCHEL2.ICHAML(1) = MCHAM1
  437. SEGDES MCHEL2
  438. MCHAM1.NOMCHE(1) = 'SCAL '
  439. MCHAM1.TYPCHE(1) = 'REAL*8 '
  440. SEGINI MELRO
  441. MCHAM1.IELVAL(1) = MELRO
  442. SEGDES MCHAM1
  443. C
  444. C**** Pression
  445. C
  446. MCHEL1 = IROF
  447. SEGINI, MCHEL2 = MCHEL1
  448. IPF = MCHEL2
  449. MCHEL2.TITCHE = 'P '
  450. C
  451. C**** MCHAM1 = MCHAML de la densite
  452. C
  453. SEGINI, MCHAM2 = MCHAM1
  454. MCHEL2.ICHAML(1) = MCHAM2
  455. SEGDES MCHEL2
  456. SEGINI MELP
  457. MCHAM2.IELVAL(1) = MELP
  458. SEGDES MCHAM2
  459. C
  460. C**** Gamma
  461. C
  462. MCHEL1 = IROF
  463. SEGINI, MCHEL2 = MCHEL1
  464. IGAMF = MCHEL2
  465. MCHEL2.TITCHE = 'GAMMA '
  466. C
  467. C**** MCHAM1 = MCHAML de la densite
  468. C
  469. SEGINI, MCHAM2 = MCHAM1
  470. MCHEL2.ICHAML(1) = MCHAM2
  471. SEGDES MCHEL2
  472. SEGINI MELGAM
  473. MCHAM2.IELVAL(1) = MELGAM
  474. SEGDES MCHAM2
  475. C
  476. C**** Recapitulatif: les MELVALs et les MPOVALs actives
  477. C
  478. C MELVNX, MELVNY, MELVNZ
  479. C MELT1X, MELT1Y, MELT1Z
  480. C MELT2X, MELT2Y, MELT2Z
  481. C
  482. C MELVUN, MELVUT -> vitesse
  483. C
  484. C MELRO -> densite
  485. C
  486. C MELP -> pression
  487. C
  488. C MELGAM -> gamma
  489. C
  490. C MPROC, MPROP -> densite
  491. C
  492. C MPVITC, MPVITP -> vitesse
  493. C
  494. C MPPC, MPPP -> pression
  495. C
  496. C MPGAMC -> gamma
  497. C
  498. C MPNORM -> normales aux faces
  499. C
  500. C
  501. C***********************************************************************
  502. C********* PREDICTION **************************************************
  503. C***********************************************************************
  504. C
  505. IF(LOGTEM)THEN
  506. C
  507. IPT3 = ICEN
  508. SEGACT IPT3
  509. NCEN = IPT3.NUM(/2)
  510. DO NLCE = 1, NCEN
  511. IF(LOGALP)THEN
  512. C
  513. C************* Gradients limités
  514. C
  515. ALPHA(1) = MELALR.VPOCHA(NLCE,1)
  516. ALPHA(2) = MELALV.VPOCHA(NLCE,1)
  517. ALPHA(3) = MELALV.VPOCHA(NLCE,2)
  518. ALPHA(4) = MELALV.VPOCHA(NLCE,3)
  519. ALPHA(5) = MELALP.VPOCHA(NLCE,1)
  520. ELSE
  521. C
  522. C************* Gradients non limités
  523. C
  524. ALPHA(1) = 1.0D0
  525. ALPHA(2) = 1.0D0
  526. ALPHA(3) = 1.0D0
  527. ALPHA(4) = 1.0D0
  528. ALPHA(5) = 1.0D0
  529. ENDIF
  530. ROG = MPROP.VPOCHA(NLCE,1)
  531. UXG = MPVITP.VPOCHA(NLCE,1)
  532. UYG = MPVITP.VPOCHA(NLCE,2)
  533. UZG = MPVITP.VPOCHA(NLCE,3)
  534. PG = MPPP.VPOCHA(NLCE,1)
  535. DROX = MPGRR.VPOCHA(NLCE,1)*ALPHA(1)
  536. DROY = MPGRR.VPOCHA(NLCE,2)*ALPHA(1)
  537. DROZ = MPGRR.VPOCHA(NLCE,3)*ALPHA(1)
  538. DUXX = MPGRV.VPOCHA(NLCE,1)*ALPHA(2)
  539. DUXY = MPGRV.VPOCHA(NLCE,2)*ALPHA(2)
  540. DUXZ = MPGRV.VPOCHA(NLCE,3)*ALPHA(2)
  541. DUYX = MPGRV.VPOCHA(NLCE,4)*ALPHA(3)
  542. DUYY = MPGRV.VPOCHA(NLCE,5)*ALPHA(3)
  543. DUYZ = MPGRV.VPOCHA(NLCE,6)*ALPHA(3)
  544. DUZX = MPGRV.VPOCHA(NLCE,7)*ALPHA(4)
  545. DUZY = MPGRV.VPOCHA(NLCE,8)*ALPHA(4)
  546. DUZZ = MPGRV.VPOCHA(NLCE,9)*ALPHA(4)
  547. DPX = MPGRP.VPOCHA(NLCE,1)*ALPHA(5)
  548. DPY = MPGRP.VPOCHA(NLCE,2)*ALPHA(5)
  549. DPZ = MPGRP.VPOCHA(NLCE,3)*ALPHA(5)
  550. GAMG = MPGAMC.VPOCHA(NLCE,1)
  551. DRO = UXG * DROX + ROG * ( DUXX + DUYY + DUZZ)
  552. & + UYG * DROY + UZG * DROZ
  553. DUX = UXG * DUXX + DPX / ROG + UYG * DUXY
  554. & + UZG * DUXZ
  555. DUY = UXG * DUYX + UYG * DUYY + DPY / ROG
  556. & + UZG * DUYZ
  557. DUZ = UXG * DUZX + UYG * DUZY + UZG * DUZZ
  558. & + DPZ / ROG
  559. DP = GAMG * PG * (DUXX + DUYY + DUZZ)
  560. & + UXG * DPX + UYG * DPY + UZG * DPZ
  561. C
  562. C************* Prediction avec/sans gradients limités
  563. C
  564. MPROP.VPOCHA(NLCE,1) = ROG - DELTAT * DRO
  565. MPVITP.VPOCHA(NLCE,1) = UXG - DELTAT * DUX
  566. MPVITP.VPOCHA(NLCE,2) = UYG - DELTAT * DUY
  567. MPVITP.VPOCHA(NLCE,3) = UZG - DELTAT * DUZ
  568. MPPP.VPOCHA(NLCE,1) = PG - DELTAT * DP
  569. ENDDO
  570. C
  571. ENDIF
  572. C
  573. C
  574. C***********************************************************************
  575. C********* CORRECTION **************************************************
  576. C***********************************************************************
  577. C
  578. C**** Boucle sur le faces
  579. C
  580. IDIMP1 = IDIM + 1
  581. DO NLCF = 1, NFAC
  582. C
  583. C******* NLCF = numero local du centre de face
  584. C NGCF = numero global du centre de face
  585. C NGCEG = numero global du centre ELT "gauche"
  586. C NLCEG = numero local du centre ELT "gauche"
  587. C NGCED = numero global du centre ELT "droite"
  588. C NLCED = numero local du centre ELT "droite"
  589. C
  590. NGCEG = IPT1.NUM(1,NLCF)
  591. NGCF = IPT1.NUM(2,NLCF)
  592. NGCED = IPT1.NUM(3,NLCF)
  593. NLCEG = MLENT1.LECT(NGCEG)
  594. NLCED = MLENT1.LECT(NGCED)
  595. C
  596. C******* TEST: IPT2.NUM(1,NLCF) = IPT1.NUM(2,NLCF)
  597. C
  598. NGCF1 = IPT2.NUM(1,NLCF)
  599. IF( NGCF1 .NE. NGCF) THEN
  600. LOGAN = .TRUE.
  601. MESERR(1:40) = 'PRET, subroutine pre122.eso '
  602. GOTO 9999
  603. ENDIF
  604. C
  605. C******* Cosinus directeurs des NORMALES aux faces
  606. C
  607. C On impose que les normales sont direct "Gauche" -> "Centre"
  608. C
  609. INDCEL = (NGCEG-1)*IDIMP1
  610. XG = XCOOR(INDCEL+1)
  611. YG = XCOOR(INDCEL+2)
  612. ZG = XCOOR(INDCEL+3)
  613. INDCEL = (NGCF-1)*IDIMP1
  614. XC = XCOOR(INDCEL + 1)
  615. YC = XCOOR(INDCEL + 2)
  616. ZC = XCOOR(INDCEL+3)
  617. INDCEL = (NGCED-1)*IDIMP1
  618. XD = XCOOR(INDCEL+1)
  619. YD = XCOOR(INDCEL+2)
  620. ZD = XCOOR(INDCEL+3)
  621. DXG = XC - XG
  622. DYG = YC - YG
  623. DZG = ZC - ZG
  624. DXD = XC - XD
  625. DYD = YC - YD
  626. DZD = ZC - ZD
  627. C
  628. C******* On calcule le sign du pruduit scalare
  629. C (Normales de Castem) * (vecteur "gauche" -> "centre")
  630. C
  631. CNX = MPNORM.VPOCHA(NLCF,7)
  632. CNY = MPNORM.VPOCHA(NLCF,8)
  633. CNZ = MPNORM.VPOCHA(NLCF,9)
  634. ORIENT = CNX * DXG + CNY * DYG + CNZ * DZG
  635. ORIENT = SIGN(1.0D0,ORIENT)
  636. IF(ORIENT .NE. 1.0D0)THEN
  637. LOGAN = .TRUE.
  638. MESERR(1:30)=
  639. & 'PRET , subroutine pre122.eso. '
  640. GOTO 9999
  641. ENDIF
  642. CNX = CNX * ORIENT
  643. CNY = CNY * ORIENT
  644. CNZ = CNZ * ORIENT
  645. C
  646. C********** Cosinus directeurs de tangente 1
  647. C
  648. CTX = MPNORM.VPOCHA(NLCF,1) * ORIENT
  649. CTY = MPNORM.VPOCHA(NLCF,2) * ORIENT
  650. CTZ = MPNORM.VPOCHA(NLCF,3) * ORIENT
  651. C
  652. C********** Cosinus directeurs de tangente 2
  653. C
  654. CVX = MPNORM.VPOCHA(NLCF,4) * ORIENT
  655. CVY = MPNORM.VPOCHA(NLCF,5) * ORIENT
  656. CVZ = MPNORM.VPOCHA(NLCF,6) * ORIENT
  657. C
  658. C
  659. C******* Les autres MELVALs
  660. C
  661. C
  662. C******* N.B.: On suppose qu'on a déjà controlle RO, P > 0
  663. C GAMMA \in (1,3)
  664. C Si non il faut le faire, en utilisant LOGBOR,
  665. C LOGNEG, VALER, VAL1, VAL2
  666. C
  667. C
  668. C
  669. C******* NGCEG = NGCED -> Mur
  670. C
  671. IF( NGCEG .EQ. NGCED)THEN
  672. C
  673. C********** Sur le mur on fait de reconstruction sur l'etat gauche
  674. C
  675. C
  676. C********** Etat gauche
  677. C
  678. VALCEL = MPROP.VPOCHA(NLCEG, 1)
  679. ALCEL = MELALR.VPOCHA(NLCEG, 1)
  680. DCEL = MPGRR.VPOCHA(NLCEG, 1)*DXG +
  681. & MPGRR.VPOCHA(NLCEG, 2)*DYG +
  682. & MPGRR.VPOCHA(NLCEG, 3)*DZG
  683. ROG = VALCEL + ALCEL * DCEL
  684. C
  685. VALCEL = MPPP.VPOCHA(NLCEG, 1)
  686. ALCEL = MELALP.VPOCHA(NLCEG, 1)
  687. DCEL = MPGRP.VPOCHA(NLCEG, 1)*DXG +
  688. & MPGRP.VPOCHA(NLCEG, 2)*DYG +
  689. & MPGRP.VPOCHA(NLCEG, 3)*DZG
  690. PG = VALCEL + ALCEL * DCEL
  691. C
  692. GAMG = MPGAMC.VPOCHA(NLCEG, 1)
  693. C
  694. VALCEL = MPVITP.VPOCHA(NLCEG, 1)
  695. ALCEL = MELALV.VPOCHA(NLCEG, 1)
  696. DCEL = MPGRV.VPOCHA(NLCEG, 1)*DXG +
  697. & MPGRV.VPOCHA(NLCEG, 2)*DYG +
  698. & MPGRV.VPOCHA(NLCEG, 3)*DZG
  699. UXG = VALCEL + ALCEL * DCEL
  700. C
  701. VALCEL = MPVITP.VPOCHA(NLCEG, 2)
  702. ALCEL = MELALV.VPOCHA(NLCEG, 2)
  703. DCEL = MPGRV.VPOCHA(NLCEG, 4)*DXG +
  704. & MPGRV.VPOCHA(NLCEG, 5)*DYG +
  705. & MPGRV.VPOCHA(NLCEG, 6)*DZG
  706. UYG = VALCEL + ALCEL * DCEL
  707. C
  708. VALCEL = MPVITP.VPOCHA(NLCEG, 3)
  709. ALCEL = MELALV.VPOCHA(NLCEG, 3)
  710. DCEL = MPGRV.VPOCHA(NLCEG, 7)*DXG +
  711. & MPGRV.VPOCHA(NLCEG, 8)*DYG +
  712. & MPGRV.VPOCHA(NLCEG, 9)*DZG
  713. UZG = VALCEL + ALCEL * DCEL
  714. C
  715. UNG = UXG * CNX + UYG * CNY + UZG * CNZ
  716. UTG = UXG * CTX + UYG * CTY + UZG * CTZ
  717. UVG = UXG * CVX + UYG * CVY + UZG * CVZ
  718. C
  719. C********** Si l'on fait pas de prediction, ce n'est pas necessaire de
  720. C controller la positivite' de la pression et de la densité; elle
  721. C est déjà garantie par la proprieté LED de limiteur; mais, vu
  722. C que le limiteur n'est pas calculé ici, mais dans un autre
  723. C operateur, on le fait
  724. C
  725. LOGI1 = (PG .LT. 0.0D0) .OR. (ROG .LT. 0.0D0)
  726. C
  727. IF(LOGI1)THEN
  728. C
  729. C************* Premier ordre en espace local
  730. C
  731. ROG = MPROC.VPOCHA(NLCEG,1)
  732. PG = MPPC.VPOCHA(NLCEG,1)
  733. UNG = MPVITC.VPOCHA(NLCEG,1)*CNX +
  734. & MPVITC.VPOCHA(NLCEG,2)*CNY +
  735. & MPVITC.VPOCHA(NLCEG,3)*CNZ
  736. UTG = MPVITC.VPOCHA(NLCEG,1)*CTX +
  737. & MPVITC.VPOCHA(NLCEG,2)*CTY +
  738. & MPVITC.VPOCHA(NLCEG,3)*CTZ
  739. UVG = MPVITC.VPOCHA(NLCEG,1)*CVX +
  740. & MPVITC.VPOCHA(NLCEG,2)*CVY +
  741. & MPVITC.VPOCHA(NLCEG,3)*CVZ
  742. ENDIF
  743. C
  744. C********** Son etat droite
  745. C
  746. ROD = ROG
  747. PD = PG
  748. GAMD = GAMG
  749. UND = -1.0D0 * UNG
  750. UTD = UTG
  751. UVD = UVG
  752. C
  753. C************* Fin cas mur
  754. C
  755. ELSE
  756. C
  757. C************* Etat gauche
  758. C
  759. VALCEL = MPROP.VPOCHA(NLCEG, 1)
  760. ALCEL = MELALR.VPOCHA(NLCEG, 1)
  761. DCEL = MPGRR.VPOCHA(NLCEG, 1)*DXG +
  762. & MPGRR.VPOCHA(NLCEG, 2)*DYG +
  763. & MPGRR.VPOCHA(NLCEG, 3)*DZG
  764. ROG = VALCEL + ALCEL * DCEL
  765. C
  766. VALCEL = MPPP.VPOCHA(NLCEG, 1)
  767. ALCEL = MELALP.VPOCHA(NLCEG, 1)
  768. DCEL = MPGRP.VPOCHA(NLCEG, 1)*DXG +
  769. & MPGRP.VPOCHA(NLCEG, 2)*DYG +
  770. & MPGRP.VPOCHA(NLCEG, 3)*DZG
  771. PG = VALCEL + ALCEL * DCEL
  772. C
  773. GAMG = MPGAMC.VPOCHA(NLCEG, 1)
  774. C
  775. VALCEL = MPVITP.VPOCHA(NLCEG, 1)
  776. ALCEL = MELALV.VPOCHA(NLCEG, 1)
  777. DCEL = MPGRV.VPOCHA(NLCEG, 1)*DXG +
  778. & MPGRV.VPOCHA(NLCEG, 2)*DYG +
  779. & MPGRV.VPOCHA(NLCEG, 3)*DZG
  780. UXG = VALCEL + ALCEL * DCEL
  781. C
  782. VALCEL = MPVITP.VPOCHA(NLCEG, 2)
  783. ALCEL = MELALV.VPOCHA(NLCEG, 2)
  784. DCEL = MPGRV.VPOCHA(NLCEG, 4)*DXG +
  785. & MPGRV.VPOCHA(NLCEG, 5)*DYG +
  786. & MPGRV.VPOCHA(NLCEG, 6)*DZG
  787. UYG = VALCEL + ALCEL * DCEL
  788. C
  789. VALCEL = MPVITP.VPOCHA(NLCEG, 3)
  790. ALCEL = MELALV.VPOCHA(NLCEG, 3)
  791. DCEL = MPGRV.VPOCHA(NLCEG, 7)*DXG +
  792. & MPGRV.VPOCHA(NLCEG, 8)*DYG +
  793. & MPGRV.VPOCHA(NLCEG, 9)*DZG
  794. UZG = VALCEL + ALCEL * DCEL
  795. C
  796. UNG = UXG * CNX + UYG * CNY + UZG * CNZ
  797. UTG = UXG * CTX + UYG * CTY + UZG * CTZ
  798. UVG = UXG * CVX + UYG * CVY + UZG * CVZ
  799. C
  800. C********** Positivite
  801. C
  802. LOGI1 = (PG .LT. 0.0D0) .OR. (ROG .LT. 0.0D0)
  803. C
  804. IF(LOGI1)THEN
  805. C
  806. C************* Premier ordre en espace local
  807. C
  808. ROG = MPROC.VPOCHA(NLCEG,1)
  809. PG = MPPC.VPOCHA(NLCEG,1)
  810. UNG = MPVITC.VPOCHA(NLCEG,1)*CNX +
  811. & MPVITC.VPOCHA(NLCEG,2)*CNY +
  812. & MPVITC.VPOCHA(NLCEG,3)*CNZ
  813. UTG = MPVITC.VPOCHA(NLCEG,1)*CTX +
  814. & MPVITC.VPOCHA(NLCEG,2)*CTY +
  815. & MPVITC.VPOCHA(NLCEG,3)*CTZ
  816. UVG = MPVITC.VPOCHA(NLCEG,1)*CVX +
  817. & MPVITC.VPOCHA(NLCEG,2)*CVY +
  818. & MPVITC.VPOCHA(NLCEG,3)*CVZ
  819. ENDIF
  820. C
  821. C********** Etat droite
  822. C
  823. VALCEL = MPROP.VPOCHA(NLCED, 1)
  824. ALCEL = MELALR.VPOCHA(NLCED, 1)
  825. DCEL = MPGRR.VPOCHA(NLCED, 1)*DXD +
  826. & MPGRR.VPOCHA(NLCED, 2)*DYD +
  827. & MPGRR.VPOCHA(NLCED, 3)*DZD
  828. ROD = VALCEL + ALCEL * DCEL
  829. C
  830. VALCEL = MPPP.VPOCHA(NLCED, 1)
  831. ALCEL = MELALP.VPOCHA(NLCED, 1)
  832. DCEL = MPGRP.VPOCHA(NLCED, 1)*DXD +
  833. & MPGRP.VPOCHA(NLCED, 2)*DYD +
  834. & MPGRP.VPOCHA(NLCED, 3)*DZD
  835. PD = VALCEL + ALCEL * DCEL
  836. C
  837. GAMD = MPGAMC.VPOCHA(NLCED, 1)
  838. C
  839. VALCEL = MPVITP.VPOCHA(NLCED, 1)
  840. ALCEL = MELALV.VPOCHA(NLCED, 1)
  841. DCEL = MPGRV.VPOCHA(NLCED, 1)*DXD +
  842. & MPGRV.VPOCHA(NLCED, 2)*DYD +
  843. & MPGRV.VPOCHA(NLCED, 3)*DZD
  844. UXD = VALCEL + ALCEL * DCEL
  845. C
  846. VALCEL = MPVITP.VPOCHA(NLCED, 2)
  847. ALCEL = MELALV.VPOCHA(NLCED, 2)
  848. DCEL = MPGRV.VPOCHA(NLCED, 4)*DXD +
  849. & MPGRV.VPOCHA(NLCED, 5)*DYD +
  850. & MPGRV.VPOCHA(NLCED, 6)*DZD
  851. UYD = VALCEL + ALCEL * DCEL
  852. C
  853. VALCEL = MPVITP.VPOCHA(NLCED, 3)
  854. ALCEL = MELALV.VPOCHA(NLCED, 3)
  855. DCEL = MPGRV.VPOCHA(NLCED, 7)*DXD +
  856. & MPGRV.VPOCHA(NLCED, 8)*DYD +
  857. & MPGRV.VPOCHA(NLCED, 9)*DZD
  858. UZD = VALCEL + ALCEL * DCEL
  859. C
  860. UND = UXD * CNX + UYD * CNY + UZD * CNZ
  861. UTD = UXD * CTX + UYD * CTY + UZD * CTZ
  862. UVD = UXD * CVX + UYD * CVY + UZD * CVZ
  863. C
  864. C********** Positivite
  865. C
  866. LOGI1 = (PD .LT. 0.0D0) .OR. (ROD .LT. 0.0D0)
  867. C
  868. IF(LOGI1)THEN
  869. C
  870. C************* Premier ordre en espace local
  871. C
  872. ROD = MPROC.VPOCHA(NLCED,1)
  873. PD = MPPC.VPOCHA(NLCED,1)
  874. UND = MPVITC.VPOCHA(NLCED,1)*CNX +
  875. & MPVITC.VPOCHA(NLCED,2)*CNY +
  876. & MPVITC.VPOCHA(NLCED,3)*CNZ
  877. UTD = MPVITC.VPOCHA(NLCED,1)*CTX +
  878. & MPVITC.VPOCHA(NLCED,2)*CTY +
  879. & MPVITC.VPOCHA(NLCED,3)*CTZ
  880. UVD = MPVITC.VPOCHA(NLCED,1)*CVX +
  881. & MPVITC.VPOCHA(NLCED,2)*CVY +
  882. & MPVITC.VPOCHA(NLCED,3)*CVZ
  883. ENDIF
  884. ENDIF
  885. C
  886. C******** Les MELVALs
  887. C
  888. C
  889. MELRO.VELCHE(1,NLCF) = ROG
  890. MELRO.VELCHE(3,NLCF) = ROD
  891. MELP.VELCHE(1,NLCF) = PG
  892. MELP.VELCHE(3,NLCF) = PD
  893. MELGAM.VELCHE(1,NLCF) = GAMG
  894. MELGAM.VELCHE(3,NLCF) = GAMD
  895. MELVUN.VELCHE(1,NLCF) = UNG
  896. MELVUN.VELCHE(3,NLCF) = UND
  897. MELVUT.VELCHE(1,NLCF) = UTG
  898. MELVUT.VELCHE(3,NLCF) = UTD
  899. MELVUV.VELCHE(1,NLCF) = UVG
  900. MELVUV.VELCHE(3,NLCF) = UVD
  901. MELVNX.VELCHE(1,NLCF) = CNX
  902. MELVNY.VELCHE(1,NLCF) = CNY
  903. MELVNZ.VELCHE(1,NLCF) = CNZ
  904. MELT1X.VELCHE(1,NLCF) = CTX
  905. MELT1Y.VELCHE(1,NLCF) = CTY
  906. MELT1Z.VELCHE(1,NLCF) = CTZ
  907. MELT2X.VELCHE(1,NLCF) = CVX
  908. MELT2Y.VELCHE(1,NLCF) = CVY
  909. MELT2Z.VELCHE(1,NLCF) = CVZ
  910. ENDDO
  911. C
  912. C**** Desactivation des SEGMENTs
  913. C
  914. SEGDES IPT1
  915. SEGDES IPT2
  916. C
  917. C**** Le MPOVALs 'Prediction' sont detruits (si existentes)
  918. C
  919. IF(LOGTEM)THEN
  920. SEGSUP MPROP
  921. SEGSUP MPVITP
  922. SEGSUP MPPP
  923. ENDIF
  924. C
  925. SEGDES MPROC
  926. SEGDES MPGRR
  927. SEGDES MELALR
  928. SEGDES MPVITC
  929. SEGDES MPGRV
  930. SEGDES MELALV
  931. SEGDES MPPC
  932. SEGDES MPGRP
  933. SEGDES MELALP
  934. SEGDES MPGAMC
  935. SEGDES MPNORM
  936. C
  937. SEGDES MELRO
  938. SEGDES MELP
  939. SEGDES MELGAM
  940. SEGDES MELVUN
  941. SEGDES MELVUT
  942. SEGDES MELVUV
  943. SEGDES MELVNX
  944. SEGDES MELVNY
  945. SEGDES MELVNZ
  946. SEGDES MELT1X
  947. SEGDES MELT1Y
  948. SEGDES MELT1Z
  949. SEGDES MELT2X
  950. SEGDES MELT2Y
  951. SEGDES MELT2Z
  952. C
  953. C
  954. C**** Destruction du MELNTI correspondance local/global
  955. C
  956. SEGSUP MLENT1
  957. C
  958. 9999 CONTINUE
  959. C
  960. RETURN
  961. END
  962.  
  963.  
  964.  
  965.  
  966.  
  967.  
  968.  
  969.  
  970.  
  971.  

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