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