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.  
  186. -INC PPARAM
  187. -INC CCOPTIO
  188. -INC SMCHPOI
  189. POINTEUR MPROC.MPOVAL, MPGRR.MPOVAL,
  190. & MPVITC.MPOVAL, MPGRV.MPOVAL,
  191. & MPPC.MPOVAL, MPGRP.MPOVAL,
  192. & MPGAMC.MPOVAL, MPNORM.MPOVAL,
  193. & MPROP.MPOVAL, MPPP.MPOVAL, MPVITP.MPOVAL
  194. -INC SMCHAML
  195. POINTEUR MELVNX.MELVAL, MELVNY.MELVAL, MELVNZ.MELVAL,
  196. & MELT1X.MELVAL, MELT1Y.MELVAL, MELT1Z.MELVAL,
  197. & MELT2X.MELVAL, MELT2Y.MELVAL, MELT2Z.MELVAL
  198. POINTEUR MELVUN.MELVAL, MELVUT.MELVAL, MELVUV.MELVAL
  199. POINTEUR MELRO.MELVAL, MELP.MELVAL,
  200. & MELGAM.MELVAL
  201. POINTEUR MELALR.MPOVAL,
  202. & MELALV.MPOVAL,
  203. & MELALP.MPOVAL
  204. -INC SMLENTI
  205. -INC SMELEME
  206. C
  207. C
  208. C**** Initialisation des parametres d'erreur déjà faite, i.e.
  209. C
  210. C LOGNEG = .FALSE.
  211. C LOGBOR = .FALSE.
  212. C MESERR = ' '
  213. C MOTERR(1:40) = MESERR(1:40)
  214. C VALER = 0.0D0
  215. C VAL1 = 0.0D0
  216. C VAL2 = 0.0D0
  217. C
  218. C
  219. C**** KRIPAD pour la correspondance global/local de centre
  220. C
  221. CALL KRIPAD(ICEN,MLENT1)
  222. C
  223. C**** MLENTI1 a MCORD.XCOORD(/1)/(IDIM+1) elements
  224. C
  225. C Si i est le numero global d'un noeud de ICEN,
  226. C MLENT1.LECT(i) contient sa position, i.e.
  227. C
  228. C I = numero global du noeud centre
  229. C MLENT1.LECT(i) = numero local du noeud centre
  230. C
  231. C MLENT1 déjà activé, i.e.
  232. C
  233. C SEGACT MLENT1
  234. C
  235. C**** Activation de CHPOINTs
  236. C
  237. C densité + grad + limiteur
  238. C vitesse + grad + limiteur
  239. C pression + grad + limiteur
  240. C gamma
  241. C cosinus directeurs des normales aux surface
  242. C
  243. CALL LICHT(IROC , MPROC , TYPE, IGEOM)
  244. CALL LICHT(IGRROC, MPGRR , TYPE, IGEOM)
  245. CALL LICHT(IVITC, MPVITC , TYPE, IGEOM)
  246. CALL LICHT(IGRVC, MPGRV , TYPE, IGEOM)
  247. CALL LICHT(IPC , MPPC , TYPE, IGEOM)
  248. CALL LICHT(IGRPC, MPGRP , TYPE, IGEOM)
  249. CALL LICHT(IGAMC, MPGAMC , TYPE, IGEOM)
  250. CALL LICHT(INORM, MPNORM , TYPE, IGEOM)
  251. C
  252. C**** Les MPOVALs 'Prediction'
  253. C
  254. IF(LOGTEM)THEN
  255. SEGINI, MPROP = MPROC
  256. SEGINI, MPPP = MPPC
  257. SEGINI, MPVITP = MPVITC
  258. ELSE
  259. MPROP = MPROC
  260. MPPP = MPPC
  261. MPVITP = MPVITC
  262. ENDIF
  263. C
  264. C**** Les Limiteurs
  265. C
  266. CALL LICHT(IALROC, MELALR , TYPE, IGEOM)
  267. CALL LICHT(IALVC, MELALV , TYPE, IGEOM)
  268. CALL LICHT(IALPC, MELALP , TYPE, IGEOM)
  269. C
  270. C
  271. C**** MPOVAL sont déjà activés i.e.:
  272. C
  273. C SEGACT MPROC
  274. C SEGACT MPGRR
  275. C SEGACT MPIALR
  276. C SEGACT MPVITC
  277. C SEGACT MPGRV
  278. C SEGACT MPIALV
  279. C SEGACT MPPC
  280. C SEGACT MPGRP
  281. C SEGACT MPIALP
  282. C SEGACT MPGAMC
  283. C SEGACT MPNORM
  284. C
  285. C**** Le MELEME FACEL
  286. C
  287. IPT1 = IFACEL
  288. IPT2 = IFACE
  289. SEGACT IPT1
  290. SEGACT IPT2
  291. NFAC = IPT1.NUM(/2)
  292. C
  293. C**** Creation de MCHAMLs contenant les etats gauche et droite,
  294. C
  295. C i.e.:
  296. C
  297. C vitesse + cosinus directors du repere local
  298. C densité
  299. C pression
  300. C gamma
  301. C
  302. C**** Cosinus directors du repere local et vitesse
  303. C
  304. C Les cosinus directeurs
  305. C
  306. N1 = 2
  307. N3 = 6
  308. L1 = 28
  309. SEGINI MCHEL1
  310. IVITF = MCHEL1
  311. MCHEL1.TITCHE = 'U '
  312. MCHEL1.IMACHE(1) = IFACE
  313. MCHEL1.IMACHE(2) = IFACEL
  314. MCHEL1.CONCHE(1) = '(n,t,v)in(x,y,z)'
  315. MCHEL1.CONCHE(2) = ' U in (n,t,v) '
  316. C
  317. C**** Valeurs des cosinus definies par respect au repair global, i.e.
  318. C
  319. MCHEL1.INFCHE(1,1) = 2
  320. MCHEL1.INFCHE(1,3) = NIFOUR
  321. MCHEL1.INFCHE(1,4) = 0
  322. MCHEL1.INFCHE(1,5) = 0
  323. MCHEL1.INFCHE(1,6) = 0
  324. MCHEL1.IFOCHE = IFOUR
  325. C
  326. C**** Valeurs de vitesse definies par respect au repair local, i.e.
  327. C
  328. MCHEL1.INFCHE(2,1) = 1
  329. MCHEL1.INFCHE(2,3) = NIFOUR
  330. MCHEL1.INFCHE(2,4) = 0
  331. MCHEL1.INFCHE(2,5) = 0
  332. MCHEL1.INFCHE(2,6) = 0
  333. C
  334. C**** Le cosinus directeurs
  335. C
  336. N1PTEL = 1
  337. N1EL = NFAC
  338. N2PTEL = 0
  339. N2EL = 0
  340. C
  341. C**** MCHAML a N2 composantes:
  342. C
  343. C cosinus directeurs du repere local (n,t1,t2)
  344. C
  345. C IDIM = 3 -> 9 composantes
  346. C
  347. N2 = 9
  348. SEGINI MCHAM1
  349. MCHEL1.ICHAML(1) = MCHAM1
  350. MCHAM1.NOMCHE(1) = 'NX '
  351. MCHAM1.NOMCHE(2) = 'NY '
  352. MCHAM1.NOMCHE(3) = 'NZ '
  353. MCHAM1.NOMCHE(4) = 'TX '
  354. MCHAM1.NOMCHE(5) = 'TY '
  355. MCHAM1.NOMCHE(6) = 'TZ '
  356. MCHAM1.NOMCHE(7) = 'VX '
  357. MCHAM1.NOMCHE(8) = 'VY '
  358. MCHAM1.NOMCHE(9) = 'VZ '
  359. MCHAM1.TYPCHE(1) = 'REAL*8 '
  360. MCHAM1.TYPCHE(2) = 'REAL*8 '
  361. MCHAM1.TYPCHE(3) = 'REAL*8 '
  362. MCHAM1.TYPCHE(4) = 'REAL*8 '
  363. MCHAM1.TYPCHE(5) = 'REAL*8 '
  364. MCHAM1.TYPCHE(6) = 'REAL*8 '
  365. MCHAM1.TYPCHE(7) = 'REAL*8 '
  366. MCHAM1.TYPCHE(8) = 'REAL*8 '
  367. MCHAM1.TYPCHE(9) = 'REAL*8 '
  368. SEGINI MELVNX
  369. SEGINI MELVNY
  370. SEGINI MELVNZ
  371. SEGINI MELT1X
  372. SEGINI MELT1Y
  373. SEGINI MELT1Z
  374. SEGINI MELT2X
  375. SEGINI MELT2Y
  376. SEGINI MELT2Z
  377. MCHAM1.IELVAL(1) = MELVNX
  378. MCHAM1.IELVAL(2) = MELVNY
  379. MCHAM1.IELVAL(3) = MELVNZ
  380. MCHAM1.IELVAL(4) = MELT1X
  381. MCHAM1.IELVAL(5) = MELT1Y
  382. MCHAM1.IELVAL(6) = MELT1Z
  383. MCHAM1.IELVAL(7) = MELT2X
  384. MCHAM1.IELVAL(8) = MELT2Y
  385. MCHAM1.IELVAL(9) = MELT2Z
  386. SEGDES MCHAM1
  387. C
  388. C**** Vitesse
  389. C
  390. N1EL = NFAC
  391. N1PTEL = 3
  392. N2EL = 0
  393. N2PTEL = 0
  394. C
  395. C**** MCHAML a N2 composantes:
  396. C
  397. C IDIM = 3 -> 3 composantes
  398. C
  399. N2 = 3
  400. SEGINI MCHAM1
  401. MCHEL1.ICHAML(2) = MCHAM1
  402. SEGDES MCHEL1
  403. MCHAM1.NOMCHE(1) = 'UN '
  404. MCHAM1.NOMCHE(2) = 'UT '
  405. MCHAM1.NOMCHE(3) = 'UV '
  406. MCHAM1.TYPCHE(1) = 'REAL*8 '
  407. MCHAM1.TYPCHE(2) = 'REAL*8 '
  408. MCHAM1.TYPCHE(3) = 'REAL*8 '
  409. SEGINI MELVUN
  410. SEGINI MELVUT
  411. SEGINI MELVUV
  412. MCHAM1.IELVAL(1) = MELVUN
  413. MCHAM1.IELVAL(2) = MELVUT
  414. MCHAM1.IELVAL(3) = MELVUV
  415. SEGDES MCHAM1
  416. C
  417. C**** Densite
  418. C
  419. N1 = 1
  420. N3 = 6
  421. L1 = 15
  422. SEGINI MCHEL2
  423. IROF = MCHEL2
  424. MCHEL2.IMACHE(1) = IFACEL
  425. MCHEL2.TITCHE = 'RO '
  426. MCHEL2.CONCHE(1) = ' '
  427. C
  428. C**** Valeurs independente du repére, i.e.
  429. C
  430. MCHEL2.INFCHE(1,1) = 0
  431. MCHEL2.INFCHE(1,3) = NIFOUR
  432. MCHEL2.INFCHE(1,4) = 0
  433. MCHEL2.INFCHE(1,5) = 0
  434. MCHEL2.INFCHE(1,6) = 0
  435. MCHEL2.IFOCHE = IFOUR
  436. N2 = 1
  437. SEGINI MCHAM1
  438. MCHEL2.ICHAML(1) = MCHAM1
  439. SEGDES MCHEL2
  440. MCHAM1.NOMCHE(1) = 'SCAL '
  441. MCHAM1.TYPCHE(1) = 'REAL*8 '
  442. SEGINI MELRO
  443. MCHAM1.IELVAL(1) = MELRO
  444. SEGDES MCHAM1
  445. C
  446. C**** Pression
  447. C
  448. MCHEL1 = IROF
  449. SEGINI, MCHEL2 = MCHEL1
  450. IPF = MCHEL2
  451. MCHEL2.TITCHE = 'P '
  452. C
  453. C**** MCHAM1 = MCHAML de la densite
  454. C
  455. SEGINI, MCHAM2 = MCHAM1
  456. MCHEL2.ICHAML(1) = MCHAM2
  457. SEGDES MCHEL2
  458. SEGINI MELP
  459. MCHAM2.IELVAL(1) = MELP
  460. SEGDES MCHAM2
  461. C
  462. C**** Gamma
  463. C
  464. MCHEL1 = IROF
  465. SEGINI, MCHEL2 = MCHEL1
  466. IGAMF = MCHEL2
  467. MCHEL2.TITCHE = 'GAMMA '
  468. C
  469. C**** MCHAM1 = MCHAML de la densite
  470. C
  471. SEGINI, MCHAM2 = MCHAM1
  472. MCHEL2.ICHAML(1) = MCHAM2
  473. SEGDES MCHEL2
  474. SEGINI MELGAM
  475. MCHAM2.IELVAL(1) = MELGAM
  476. SEGDES MCHAM2
  477. C
  478. C**** Recapitulatif: les MELVALs et les MPOVALs actives
  479. C
  480. C MELVNX, MELVNY, MELVNZ
  481. C MELT1X, MELT1Y, MELT1Z
  482. C MELT2X, MELT2Y, MELT2Z
  483. C
  484. C MELVUN, MELVUT -> vitesse
  485. C
  486. C MELRO -> densite
  487. C
  488. C MELP -> pression
  489. C
  490. C MELGAM -> gamma
  491. C
  492. C MPROC, MPROP -> densite
  493. C
  494. C MPVITC, MPVITP -> vitesse
  495. C
  496. C MPPC, MPPP -> pression
  497. C
  498. C MPGAMC -> gamma
  499. C
  500. C MPNORM -> normales aux faces
  501. C
  502. C
  503. C***********************************************************************
  504. C********* PREDICTION **************************************************
  505. C***********************************************************************
  506. C
  507. IF(LOGTEM)THEN
  508. C
  509. IPT3 = ICEN
  510. SEGACT IPT3
  511. NCEN = IPT3.NUM(/2)
  512. DO NLCE = 1, NCEN
  513. IF(LOGALP)THEN
  514. C
  515. C************* Gradients limités
  516. C
  517. ALPHA(1) = MELALR.VPOCHA(NLCE,1)
  518. ALPHA(2) = MELALV.VPOCHA(NLCE,1)
  519. ALPHA(3) = MELALV.VPOCHA(NLCE,2)
  520. ALPHA(4) = MELALV.VPOCHA(NLCE,3)
  521. ALPHA(5) = MELALP.VPOCHA(NLCE,1)
  522. ELSE
  523. C
  524. C************* Gradients non limités
  525. C
  526. ALPHA(1) = 1.0D0
  527. ALPHA(2) = 1.0D0
  528. ALPHA(3) = 1.0D0
  529. ALPHA(4) = 1.0D0
  530. ALPHA(5) = 1.0D0
  531. ENDIF
  532. ROG = MPROP.VPOCHA(NLCE,1)
  533. UXG = MPVITP.VPOCHA(NLCE,1)
  534. UYG = MPVITP.VPOCHA(NLCE,2)
  535. UZG = MPVITP.VPOCHA(NLCE,3)
  536. PG = MPPP.VPOCHA(NLCE,1)
  537. DROX = MPGRR.VPOCHA(NLCE,1)*ALPHA(1)
  538. DROY = MPGRR.VPOCHA(NLCE,2)*ALPHA(1)
  539. DROZ = MPGRR.VPOCHA(NLCE,3)*ALPHA(1)
  540. DUXX = MPGRV.VPOCHA(NLCE,1)*ALPHA(2)
  541. DUXY = MPGRV.VPOCHA(NLCE,2)*ALPHA(2)
  542. DUXZ = MPGRV.VPOCHA(NLCE,3)*ALPHA(2)
  543. DUYX = MPGRV.VPOCHA(NLCE,4)*ALPHA(3)
  544. DUYY = MPGRV.VPOCHA(NLCE,5)*ALPHA(3)
  545. DUYZ = MPGRV.VPOCHA(NLCE,6)*ALPHA(3)
  546. DUZX = MPGRV.VPOCHA(NLCE,7)*ALPHA(4)
  547. DUZY = MPGRV.VPOCHA(NLCE,8)*ALPHA(4)
  548. DUZZ = MPGRV.VPOCHA(NLCE,9)*ALPHA(4)
  549. DPX = MPGRP.VPOCHA(NLCE,1)*ALPHA(5)
  550. DPY = MPGRP.VPOCHA(NLCE,2)*ALPHA(5)
  551. DPZ = MPGRP.VPOCHA(NLCE,3)*ALPHA(5)
  552. GAMG = MPGAMC.VPOCHA(NLCE,1)
  553. DRO = UXG * DROX + ROG * ( DUXX + DUYY + DUZZ)
  554. & + UYG * DROY + UZG * DROZ
  555. DUX = UXG * DUXX + DPX / ROG + UYG * DUXY
  556. & + UZG * DUXZ
  557. DUY = UXG * DUYX + UYG * DUYY + DPY / ROG
  558. & + UZG * DUYZ
  559. DUZ = UXG * DUZX + UYG * DUZY + UZG * DUZZ
  560. & + DPZ / ROG
  561. DP = GAMG * PG * (DUXX + DUYY + DUZZ)
  562. & + UXG * DPX + UYG * DPY + UZG * DPZ
  563. C
  564. C************* Prediction avec/sans gradients limités
  565. C
  566. MPROP.VPOCHA(NLCE,1) = ROG - DELTAT * DRO
  567. MPVITP.VPOCHA(NLCE,1) = UXG - DELTAT * DUX
  568. MPVITP.VPOCHA(NLCE,2) = UYG - DELTAT * DUY
  569. MPVITP.VPOCHA(NLCE,3) = UZG - DELTAT * DUZ
  570. MPPP.VPOCHA(NLCE,1) = PG - DELTAT * DP
  571. ENDDO
  572. C
  573. ENDIF
  574. C
  575. C
  576. C***********************************************************************
  577. C********* CORRECTION **************************************************
  578. C***********************************************************************
  579. C
  580. C**** Boucle sur le faces
  581. C
  582. IDIMP1 = IDIM + 1
  583. DO NLCF = 1, NFAC
  584. C
  585. C******* NLCF = numero local du centre de face
  586. C NGCF = numero global du centre de face
  587. C NGCEG = numero global du centre ELT "gauche"
  588. C NLCEG = numero local du centre ELT "gauche"
  589. C NGCED = numero global du centre ELT "droite"
  590. C NLCED = numero local du centre ELT "droite"
  591. C
  592. NGCEG = IPT1.NUM(1,NLCF)
  593. NGCF = IPT1.NUM(2,NLCF)
  594. NGCED = IPT1.NUM(3,NLCF)
  595. NLCEG = MLENT1.LECT(NGCEG)
  596. NLCED = MLENT1.LECT(NGCED)
  597. C
  598. C******* TEST: IPT2.NUM(1,NLCF) = IPT1.NUM(2,NLCF)
  599. C
  600. NGCF1 = IPT2.NUM(1,NLCF)
  601. IF( NGCF1 .NE. NGCF) THEN
  602. LOGAN = .TRUE.
  603. MESERR(1:40) = 'PRET, subroutine pre122.eso '
  604. GOTO 9999
  605. ENDIF
  606. C
  607. C******* Cosinus directeurs des NORMALES aux faces
  608. C
  609. C On impose que les normales sont direct "Gauche" -> "Centre"
  610. C
  611. INDCEL = (NGCEG-1)*IDIMP1
  612. XG = XCOOR(INDCEL+1)
  613. YG = XCOOR(INDCEL+2)
  614. ZG = XCOOR(INDCEL+3)
  615. INDCEL = (NGCF-1)*IDIMP1
  616. XC = XCOOR(INDCEL + 1)
  617. YC = XCOOR(INDCEL + 2)
  618. ZC = XCOOR(INDCEL+3)
  619. INDCEL = (NGCED-1)*IDIMP1
  620. XD = XCOOR(INDCEL+1)
  621. YD = XCOOR(INDCEL+2)
  622. ZD = XCOOR(INDCEL+3)
  623. DXG = XC - XG
  624. DYG = YC - YG
  625. DZG = ZC - ZG
  626. DXD = XC - XD
  627. DYD = YC - YD
  628. DZD = ZC - ZD
  629. C
  630. C******* On calcule le sign du pruduit scalare
  631. C (Normales de Castem) * (vecteur "gauche" -> "centre")
  632. C
  633. CNX = MPNORM.VPOCHA(NLCF,7)
  634. CNY = MPNORM.VPOCHA(NLCF,8)
  635. CNZ = MPNORM.VPOCHA(NLCF,9)
  636. ORIENT = CNX * DXG + CNY * DYG + CNZ * DZG
  637. ORIENT = SIGN(1.0D0,ORIENT)
  638. IF(ORIENT .NE. 1.0D0)THEN
  639. LOGAN = .TRUE.
  640. MESERR(1:30)=
  641. & 'PRET , subroutine pre122.eso. '
  642. GOTO 9999
  643. ENDIF
  644. CNX = CNX * ORIENT
  645. CNY = CNY * ORIENT
  646. CNZ = CNZ * ORIENT
  647. C
  648. C********** Cosinus directeurs de tangente 1
  649. C
  650. CTX = MPNORM.VPOCHA(NLCF,1) * ORIENT
  651. CTY = MPNORM.VPOCHA(NLCF,2) * ORIENT
  652. CTZ = MPNORM.VPOCHA(NLCF,3) * ORIENT
  653. C
  654. C********** Cosinus directeurs de tangente 2
  655. C
  656. CVX = MPNORM.VPOCHA(NLCF,4) * ORIENT
  657. CVY = MPNORM.VPOCHA(NLCF,5) * ORIENT
  658. CVZ = MPNORM.VPOCHA(NLCF,6) * ORIENT
  659. C
  660. C
  661. C******* Les autres MELVALs
  662. C
  663. C
  664. C******* N.B.: On suppose qu'on a déjà controlle RO, P > 0
  665. C GAMMA \in (1,3)
  666. C Si non il faut le faire, en utilisant LOGBOR,
  667. C LOGNEG, VALER, VAL1, VAL2
  668. C
  669. C
  670. C
  671. C******* NGCEG = NGCED -> Mur
  672. C
  673. IF( NGCEG .EQ. NGCED)THEN
  674. C
  675. C********** Sur le mur on fait de reconstruction sur l'etat gauche
  676. C
  677. C
  678. C********** Etat gauche
  679. C
  680. VALCEL = MPROP.VPOCHA(NLCEG, 1)
  681. ALCEL = MELALR.VPOCHA(NLCEG, 1)
  682. DCEL = MPGRR.VPOCHA(NLCEG, 1)*DXG +
  683. & MPGRR.VPOCHA(NLCEG, 2)*DYG +
  684. & MPGRR.VPOCHA(NLCEG, 3)*DZG
  685. ROG = VALCEL + ALCEL * DCEL
  686. C
  687. VALCEL = MPPP.VPOCHA(NLCEG, 1)
  688. ALCEL = MELALP.VPOCHA(NLCEG, 1)
  689. DCEL = MPGRP.VPOCHA(NLCEG, 1)*DXG +
  690. & MPGRP.VPOCHA(NLCEG, 2)*DYG +
  691. & MPGRP.VPOCHA(NLCEG, 3)*DZG
  692. PG = VALCEL + ALCEL * DCEL
  693. C
  694. GAMG = MPGAMC.VPOCHA(NLCEG, 1)
  695. C
  696. VALCEL = MPVITP.VPOCHA(NLCEG, 1)
  697. ALCEL = MELALV.VPOCHA(NLCEG, 1)
  698. DCEL = MPGRV.VPOCHA(NLCEG, 1)*DXG +
  699. & MPGRV.VPOCHA(NLCEG, 2)*DYG +
  700. & MPGRV.VPOCHA(NLCEG, 3)*DZG
  701. UXG = VALCEL + ALCEL * DCEL
  702. C
  703. VALCEL = MPVITP.VPOCHA(NLCEG, 2)
  704. ALCEL = MELALV.VPOCHA(NLCEG, 2)
  705. DCEL = MPGRV.VPOCHA(NLCEG, 4)*DXG +
  706. & MPGRV.VPOCHA(NLCEG, 5)*DYG +
  707. & MPGRV.VPOCHA(NLCEG, 6)*DZG
  708. UYG = VALCEL + ALCEL * DCEL
  709. C
  710. VALCEL = MPVITP.VPOCHA(NLCEG, 3)
  711. ALCEL = MELALV.VPOCHA(NLCEG, 3)
  712. DCEL = MPGRV.VPOCHA(NLCEG, 7)*DXG +
  713. & MPGRV.VPOCHA(NLCEG, 8)*DYG +
  714. & MPGRV.VPOCHA(NLCEG, 9)*DZG
  715. UZG = VALCEL + ALCEL * DCEL
  716. C
  717. UNG = UXG * CNX + UYG * CNY + UZG * CNZ
  718. UTG = UXG * CTX + UYG * CTY + UZG * CTZ
  719. UVG = UXG * CVX + UYG * CVY + UZG * CVZ
  720. C
  721. C********** Si l'on fait pas de prediction, ce n'est pas necessaire de
  722. C controller la positivite' de la pression et de la densité; elle
  723. C est déjà garantie par la proprieté LED de limiteur; mais, vu
  724. C que le limiteur n'est pas calculé ici, mais dans un autre
  725. C operateur, on le fait
  726. C
  727. LOGI1 = (PG .LT. 0.0D0) .OR. (ROG .LT. 0.0D0)
  728. C
  729. IF(LOGI1)THEN
  730. C
  731. C************* Premier ordre en espace local
  732. C
  733. ROG = MPROC.VPOCHA(NLCEG,1)
  734. PG = MPPC.VPOCHA(NLCEG,1)
  735. UNG = MPVITC.VPOCHA(NLCEG,1)*CNX +
  736. & MPVITC.VPOCHA(NLCEG,2)*CNY +
  737. & MPVITC.VPOCHA(NLCEG,3)*CNZ
  738. UTG = MPVITC.VPOCHA(NLCEG,1)*CTX +
  739. & MPVITC.VPOCHA(NLCEG,2)*CTY +
  740. & MPVITC.VPOCHA(NLCEG,3)*CTZ
  741. UVG = MPVITC.VPOCHA(NLCEG,1)*CVX +
  742. & MPVITC.VPOCHA(NLCEG,2)*CVY +
  743. & MPVITC.VPOCHA(NLCEG,3)*CVZ
  744. ENDIF
  745. C
  746. C********** Son etat droite
  747. C
  748. ROD = ROG
  749. PD = PG
  750. GAMD = GAMG
  751. UND = -1.0D0 * UNG
  752. UTD = UTG
  753. UVD = UVG
  754. C
  755. C************* Fin cas mur
  756. C
  757. ELSE
  758. C
  759. C************* Etat gauche
  760. C
  761. VALCEL = MPROP.VPOCHA(NLCEG, 1)
  762. ALCEL = MELALR.VPOCHA(NLCEG, 1)
  763. DCEL = MPGRR.VPOCHA(NLCEG, 1)*DXG +
  764. & MPGRR.VPOCHA(NLCEG, 2)*DYG +
  765. & MPGRR.VPOCHA(NLCEG, 3)*DZG
  766. ROG = VALCEL + ALCEL * DCEL
  767. C
  768. VALCEL = MPPP.VPOCHA(NLCEG, 1)
  769. ALCEL = MELALP.VPOCHA(NLCEG, 1)
  770. DCEL = MPGRP.VPOCHA(NLCEG, 1)*DXG +
  771. & MPGRP.VPOCHA(NLCEG, 2)*DYG +
  772. & MPGRP.VPOCHA(NLCEG, 3)*DZG
  773. PG = VALCEL + ALCEL * DCEL
  774. C
  775. GAMG = MPGAMC.VPOCHA(NLCEG, 1)
  776. C
  777. VALCEL = MPVITP.VPOCHA(NLCEG, 1)
  778. ALCEL = MELALV.VPOCHA(NLCEG, 1)
  779. DCEL = MPGRV.VPOCHA(NLCEG, 1)*DXG +
  780. & MPGRV.VPOCHA(NLCEG, 2)*DYG +
  781. & MPGRV.VPOCHA(NLCEG, 3)*DZG
  782. UXG = VALCEL + ALCEL * DCEL
  783. C
  784. VALCEL = MPVITP.VPOCHA(NLCEG, 2)
  785. ALCEL = MELALV.VPOCHA(NLCEG, 2)
  786. DCEL = MPGRV.VPOCHA(NLCEG, 4)*DXG +
  787. & MPGRV.VPOCHA(NLCEG, 5)*DYG +
  788. & MPGRV.VPOCHA(NLCEG, 6)*DZG
  789. UYG = VALCEL + ALCEL * DCEL
  790. C
  791. VALCEL = MPVITP.VPOCHA(NLCEG, 3)
  792. ALCEL = MELALV.VPOCHA(NLCEG, 3)
  793. DCEL = MPGRV.VPOCHA(NLCEG, 7)*DXG +
  794. & MPGRV.VPOCHA(NLCEG, 8)*DYG +
  795. & MPGRV.VPOCHA(NLCEG, 9)*DZG
  796. UZG = VALCEL + ALCEL * DCEL
  797. C
  798. UNG = UXG * CNX + UYG * CNY + UZG * CNZ
  799. UTG = UXG * CTX + UYG * CTY + UZG * CTZ
  800. UVG = UXG * CVX + UYG * CVY + UZG * CVZ
  801. C
  802. C********** Positivite
  803. C
  804. LOGI1 = (PG .LT. 0.0D0) .OR. (ROG .LT. 0.0D0)
  805. C
  806. IF(LOGI1)THEN
  807. C
  808. C************* Premier ordre en espace local
  809. C
  810. ROG = MPROC.VPOCHA(NLCEG,1)
  811. PG = MPPC.VPOCHA(NLCEG,1)
  812. UNG = MPVITC.VPOCHA(NLCEG,1)*CNX +
  813. & MPVITC.VPOCHA(NLCEG,2)*CNY +
  814. & MPVITC.VPOCHA(NLCEG,3)*CNZ
  815. UTG = MPVITC.VPOCHA(NLCEG,1)*CTX +
  816. & MPVITC.VPOCHA(NLCEG,2)*CTY +
  817. & MPVITC.VPOCHA(NLCEG,3)*CTZ
  818. UVG = MPVITC.VPOCHA(NLCEG,1)*CVX +
  819. & MPVITC.VPOCHA(NLCEG,2)*CVY +
  820. & MPVITC.VPOCHA(NLCEG,3)*CVZ
  821. ENDIF
  822. C
  823. C********** Etat droite
  824. C
  825. VALCEL = MPROP.VPOCHA(NLCED, 1)
  826. ALCEL = MELALR.VPOCHA(NLCED, 1)
  827. DCEL = MPGRR.VPOCHA(NLCED, 1)*DXD +
  828. & MPGRR.VPOCHA(NLCED, 2)*DYD +
  829. & MPGRR.VPOCHA(NLCED, 3)*DZD
  830. ROD = VALCEL + ALCEL * DCEL
  831. C
  832. VALCEL = MPPP.VPOCHA(NLCED, 1)
  833. ALCEL = MELALP.VPOCHA(NLCED, 1)
  834. DCEL = MPGRP.VPOCHA(NLCED, 1)*DXD +
  835. & MPGRP.VPOCHA(NLCED, 2)*DYD +
  836. & MPGRP.VPOCHA(NLCED, 3)*DZD
  837. PD = VALCEL + ALCEL * DCEL
  838. C
  839. GAMD = MPGAMC.VPOCHA(NLCED, 1)
  840. C
  841. VALCEL = MPVITP.VPOCHA(NLCED, 1)
  842. ALCEL = MELALV.VPOCHA(NLCED, 1)
  843. DCEL = MPGRV.VPOCHA(NLCED, 1)*DXD +
  844. & MPGRV.VPOCHA(NLCED, 2)*DYD +
  845. & MPGRV.VPOCHA(NLCED, 3)*DZD
  846. UXD = VALCEL + ALCEL * DCEL
  847. C
  848. VALCEL = MPVITP.VPOCHA(NLCED, 2)
  849. ALCEL = MELALV.VPOCHA(NLCED, 2)
  850. DCEL = MPGRV.VPOCHA(NLCED, 4)*DXD +
  851. & MPGRV.VPOCHA(NLCED, 5)*DYD +
  852. & MPGRV.VPOCHA(NLCED, 6)*DZD
  853. UYD = VALCEL + ALCEL * DCEL
  854. C
  855. VALCEL = MPVITP.VPOCHA(NLCED, 3)
  856. ALCEL = MELALV.VPOCHA(NLCED, 3)
  857. DCEL = MPGRV.VPOCHA(NLCED, 7)*DXD +
  858. & MPGRV.VPOCHA(NLCED, 8)*DYD +
  859. & MPGRV.VPOCHA(NLCED, 9)*DZD
  860. UZD = VALCEL + ALCEL * DCEL
  861. C
  862. UND = UXD * CNX + UYD * CNY + UZD * CNZ
  863. UTD = UXD * CTX + UYD * CTY + UZD * CTZ
  864. UVD = UXD * CVX + UYD * CVY + UZD * CVZ
  865. C
  866. C********** Positivite
  867. C
  868. LOGI1 = (PD .LT. 0.0D0) .OR. (ROD .LT. 0.0D0)
  869. C
  870. IF(LOGI1)THEN
  871. C
  872. C************* Premier ordre en espace local
  873. C
  874. ROD = MPROC.VPOCHA(NLCED,1)
  875. PD = MPPC.VPOCHA(NLCED,1)
  876. UND = MPVITC.VPOCHA(NLCED,1)*CNX +
  877. & MPVITC.VPOCHA(NLCED,2)*CNY +
  878. & MPVITC.VPOCHA(NLCED,3)*CNZ
  879. UTD = MPVITC.VPOCHA(NLCED,1)*CTX +
  880. & MPVITC.VPOCHA(NLCED,2)*CTY +
  881. & MPVITC.VPOCHA(NLCED,3)*CTZ
  882. UVD = MPVITC.VPOCHA(NLCED,1)*CVX +
  883. & MPVITC.VPOCHA(NLCED,2)*CVY +
  884. & MPVITC.VPOCHA(NLCED,3)*CVZ
  885. ENDIF
  886. ENDIF
  887. C
  888. C******** Les MELVALs
  889. C
  890. C
  891. MELRO.VELCHE(1,NLCF) = ROG
  892. MELRO.VELCHE(3,NLCF) = ROD
  893. MELP.VELCHE(1,NLCF) = PG
  894. MELP.VELCHE(3,NLCF) = PD
  895. MELGAM.VELCHE(1,NLCF) = GAMG
  896. MELGAM.VELCHE(3,NLCF) = GAMD
  897. MELVUN.VELCHE(1,NLCF) = UNG
  898. MELVUN.VELCHE(3,NLCF) = UND
  899. MELVUT.VELCHE(1,NLCF) = UTG
  900. MELVUT.VELCHE(3,NLCF) = UTD
  901. MELVUV.VELCHE(1,NLCF) = UVG
  902. MELVUV.VELCHE(3,NLCF) = UVD
  903. MELVNX.VELCHE(1,NLCF) = CNX
  904. MELVNY.VELCHE(1,NLCF) = CNY
  905. MELVNZ.VELCHE(1,NLCF) = CNZ
  906. MELT1X.VELCHE(1,NLCF) = CTX
  907. MELT1Y.VELCHE(1,NLCF) = CTY
  908. MELT1Z.VELCHE(1,NLCF) = CTZ
  909. MELT2X.VELCHE(1,NLCF) = CVX
  910. MELT2Y.VELCHE(1,NLCF) = CVY
  911. MELT2Z.VELCHE(1,NLCF) = CVZ
  912. ENDDO
  913. C
  914. C**** Desactivation des SEGMENTs
  915. C
  916. SEGDES IPT1
  917. SEGDES IPT2
  918. C
  919. C**** Le MPOVALs 'Prediction' sont detruits (si existentes)
  920. C
  921. IF(LOGTEM)THEN
  922. SEGSUP MPROP
  923. SEGSUP MPVITP
  924. SEGSUP MPPP
  925. ENDIF
  926. C
  927. SEGDES MPROC
  928. SEGDES MPGRR
  929. SEGDES MELALR
  930. SEGDES MPVITC
  931. SEGDES MPGRV
  932. SEGDES MELALV
  933. SEGDES MPPC
  934. SEGDES MPGRP
  935. SEGDES MELALP
  936. SEGDES MPGAMC
  937. SEGDES MPNORM
  938. C
  939. SEGDES MELRO
  940. SEGDES MELP
  941. SEGDES MELGAM
  942. SEGDES MELVUN
  943. SEGDES MELVUT
  944. SEGDES MELVUV
  945. SEGDES MELVNX
  946. SEGDES MELVNY
  947. SEGDES MELVNZ
  948. SEGDES MELT1X
  949. SEGDES MELT1Y
  950. SEGDES MELT1Z
  951. SEGDES MELT2X
  952. SEGDES MELT2Y
  953. SEGDES MELT2Z
  954. C
  955. C
  956. C**** Destruction du MELNTI correspondance local/global
  957. C
  958. SEGSUP MLENT1
  959. C
  960. 9999 CONTINUE
  961. C
  962. RETURN
  963. END
  964.  
  965.  
  966.  
  967.  
  968.  
  969.  
  970.  
  971.  
  972.  
  973.  

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