1 : C PRE711 SOURCE BECC 11/05/18 21:15:27 6974 2 : SUBROUTINE PRE711(NESP,MLMESP, 3 : & ICEN,IFACE,IFACEL,INORM, 4 : & IPHI, IGRPHI, ILIPHI, 5 : & IRN1, IGRRN1, ILIRN1, 6 : & IVN1, IGRVN1, ILIVN1, 7 : & IPN1, IGRPN1, ILIPN1, 8 : & IYMA, IGRYMA, ILIYMA, 9 : & IALC, IGRALC, ILIALC, 10 : & IPHIF, IRN1F, IVN1F, IPN1F, IYF, IALF) 11 : C************************************************************************ 12 : C 13 : C PROJET : CASTEM 2000 14 : C 15 : C NOM : PRE711 16 : C 17 : C DESCRIPTION : Voir PRE71 18 : C 19 : C LANGAGE : FORTRAN 77 + ESOPE 2000 (avec extensions CISI) 20 : C 21 : C AUTEUR : A. BECCANTINI, DEN/DM2S/SFME/LTMF 22 : C 23 : C************************************************************************ 24 : C 25 : C ENTREES 26 : C 27 : C NESP : number of the species involved in the EULER equations. 28 : C 29 : C MLMESP : MLMOTS object; names of the species involved in the 30 : C EULER equations. 31 : C 32 : C 1) Pointeurs de MELEMEs et de CHPOINTs de la table DOMAINE 33 : C 34 : C ICEN : MELEME de 'POI1' SPG des CENTRES 35 : C 36 : C IFACE : MELEME de 'POI1' SPG des FACES 37 : C 38 : C IFACEL : MELEME de 'SEG3' avec 39 : C CENTRE d'Elt "gauche" 40 : C CENTRE de Face 41 : C CENTRE d'Elt "droite" 42 : C 43 : C N.B. = IFACE.NUM(i,1) = IFACEL.NUM(i,2) 44 : C 45 : C INORM : CHPOINT des cosinus directeurs de normales aux faces 46 : C 47 : C 2) Autres pointeurs 48 : C 49 : C IPHI, IGRPHI, ILIPHI, 50 : C CHPOINT "CENTRE" de phi, gradient et limiteur 51 : C 52 : C IRN1, IGRRN1, ILIRN1, 53 : C CHPOINT "CENTRE" de densité, gradient et limiteur 54 : C 55 : C IVN1, IGRVN1, ILIVN1, 56 : C CHPOINT "CENTRE" de vitesse, gradient et limiteur 57 : C 58 : C IPN1, IGRPN1, ILIPN1. 59 : C CHPOINT "CENTRE" de pression, gradients et limiteurs 60 : C 61 : C IYMA, IGRYMA, ILIYMA, 62 : C CHPOINT "CENTRE" de Y, gradient et limiteur 63 : C 64 : C IALC, IGRALC, ILIALC, 65 : C CHPOINT "CENTRE" de ALPHA, gradient et limiteur 66 : C 67 : C SORTIES 68 : C 69 : C IPHIF, IRN1F, IVN1F, IPN1F, IYF, IALF 70 : C MCHAMLs definis sur le MELEME de pointeur IFACEL 71 : C 72 : C************************************************************************ 73 : C 74 : C HISTORIQUE (Anomalies et modifications éventuelles) 75 : C 76 : C HISTORIQUE : Crée le 03.12.10. 77 : C 78 : C************************************************************************ 79 : C 80 : C 81 : C ATTENTION: Cet programme marche si le MAILLAGE est convex; 82 : C si non il faut changer l'algoritme de calcul de 83 : C l'orientation des normales aux faces. 84 : C 85 : C La positivité n'est pas controlle parce que c'est déjà fait 86 : C dans l'operateur PRIM 87 : C 88 : C 89 : C************************************************************************ 90 : C 91 : IMPLICIT INTEGER(I-N) 92 : C 93 : C**** Les variables 94 : C 95 : INTEGER NESP, I1 96 : & , ICEN, IFACE, IFACEL, INORM 97 : & , IPHI, IGRPHI, ILIPHI 98 : & , IRN1, IGRRN1, ILIRN1 99 : & , IVN1, IGRVN1, ILIVN1 100 : & , IPN1, IGRPN1, ILIPN1 101 : & , IYMA, IGRYMA, ILIYMA 102 : & , IALC, IGRALC, ILIALC 103 : & , IPHIF, IRN1F 104 : & , IVN1F, IPN1F 105 : & , IYF, IALF 106 : & , IGEOM, NFAC 107 : & , N1PTEL, N1EL, N2PTEL, N2EL, N2, N1, N3, L1 108 : & , NLCF, NGCF, NGCEG, NLCEG, NGCED, NLCED, NGCF1 109 : & , IDIMP1, INDCEL 110 : 111 : REAL*8 XG, YG, XC, YC, XD, YD 112 : & ,DXG, DYG, DXD, DYD 113 : & , CNX, CNY, CTX, CTY, ORIENT 114 : & , PHIG, RN1G, PN1G 115 : & , UX1G, UY1G 116 : & , UN1G, UT1G 117 : & , PHID, RN1D, PN1D 118 : & , UX1D, UY1D 119 : & , UN1D, UT1D 120 : & , VALCEL, DCEL, LIMCEL 121 : CHARACTER*(40) MESERR 122 : CHARACTER*(8) TYPE 123 : CHARACTER*(4) CARCEL 124 : LOGICAL LOGI1 125 : C 126 : C**** Segments des fractions massiques gauche et droit 127 : C 128 : INTEGER NS 129 : SEGMENT FRAMAS 130 : REAL*8 FRAMG(NS), FRAMD(NS) 131 : ENDSEGMENT 132 : POINTEUR ALPHA.FRAMAS 133 : C 134 : C**** Les Includes 135 : C 136 : -INC SMCOORD 137 : -INC CCOPTIO 138 : -INC SMCHPOI 139 : POINTEUR MPPHI.MPOVAL, MPGPHI.MPOVAL, MPLPHI.MPOVAL 140 : POINTEUR MPRN1.MPOVAL, MPGRN1.MPOVAL, MPLRN1.MPOVAL 141 : POINTEUR MPVN1.MPOVAL, MPGVN1.MPOVAL, MPLVN1.MPOVAL 142 : POINTEUR MPPN1.MPOVAL, MPGPN1.MPOVAL, MPLPN1.MPOVAL 143 : POINTEUR MPYMA.MPOVAL, MPGYMA.MPOVAL, MPLYMA.MPOVAL 144 : POINTEUR MPALC.MPOVAL, MPGALC.MPOVAL, MPLALC.MPOVAL 145 : POINTEUR MPNORM.MPOVAL 146 : -INC SMCHAML 147 : C Melval des cosinus directeurs 148 : POINTEUR MELVNX.MELVAL, MELVNY.MELVAL, 149 : & MELT1X.MELVAL, MELT1Y.MELVAL 150 : C Melval des vitesses 151 : POINTEUR MEVUN1.MELVAL, MEVUT1.MELVAL 152 : C Melval des densités, pressions, alphas 153 : POINTEUR MELRN1.MELVAL 154 : POINTEUR MELPN1.MELVAL 155 : POINTEUR MELPHI.MELVAL 156 : POINTEUR MCHAMY.MCHAML 157 : POINTEUR MCHAMA.MCHAML 158 : -INC SMLENTI 159 : -INC SMELEME 160 : -INC SMLMOTS 161 : POINTEUR MLMESP.MLMOTS 162 : C 163 : LOGI1 = .FALSE. 164 : C 165 : C**** KRIPAD pour la correspondance global/local de centre 166 : C 166 : CALL KRIPAD(ICEN,MLENT1) 167 : C 168 : C**** MLENTI1 a MCORD.XCOORD(/1)/(IDIM+1) elements 169 : C 170 : C Si i est le numero global d'un noeud de ICEN, 171 : C MLENT1.LECT(i) contient sa position, i.e. 172 : C 173 : C I = numero global du noeud centre 174 : C MLENT1.LECT(i) = numero local du noeud centre 175 : C 176 : C MLENT1 déjà activé, i.e. 177 : C 178 : C SEGACT MLENT1 179 : C 180 : C**** Activation de CHPOINTs 181 : C 182 : C phi + grad + limiteur 183 : C densité + grad + limiteur 184 : C vitesse + grad + limiteur 185 : C pression + grad + limiteur 186 : C ymass + grad + limiteur 187 : C alpha + grad + limiteur 188 : C cosinus directeurs des normales aux surface 189 : C 189 : CALL LICHT(IPHI , MPPHI , TYPE, IGEOM) 190 : C SEGACT MPPHI 190 : CALL LICHT(IGRPHI , MPGPHI , TYPE, IGEOM) 191 : C SEGACT MPGPHI 191 : CALL LICHT(ILIPHI , MPLPHI , TYPE, IGEOM) 192 : C SEGACT MPLPHI 192 : CALL LICHT(IRN1 , MPRN1 , TYPE, IGEOM) 193 : C SEGACT MPRN1 193 : CALL LICHT(IGRRN1 , MPGRN1 , TYPE, IGEOM) 194 : C SEGACT MPGRN1 194 : CALL LICHT(ILIRN1 , MPLRN1 , TYPE, IGEOM) 195 : C SEGACT MPLRN1 195 : CALL LICHT(IVN1 , MPVN1 , TYPE, IGEOM) 196 : C SEGACT MPVN1 196 : CALL LICHT(IGRVN1 , MPGVN1 , TYPE, IGEOM) 197 : C SEGACT MPGVN1 197 : CALL LICHT(ILIVN1 , MPLVN1 , TYPE, IGEOM) 198 : C SEGACT MPLVN1 198 : CALL LICHT(IPN1 , MPPN1 , TYPE, IGEOM) 199 : C SEGACT MPPN1 199 : CALL LICHT(IGRPN1 , MPGPN1 , TYPE, IGEOM) 200 : C SEGACT MPGPN1 200 : CALL LICHT(ILIPN1 , MPLPN1 , TYPE, IGEOM) 201 : C SEGACT MPLPN1 202 : IF (NESP .GE. 1)THEN 202 : CALL LICHT(IYMA , MPYMA , TYPE, IGEOM) 203 : C SEGACT MPYMA 203 : CALL LICHT(IGRYMA , MPGYMA , TYPE, IGEOM) 204 : C SEGACT MPGYMA 204 : CALL LICHT(ILIYMA , MPLYMA , TYPE, IGEOM) 205 : C SEGACT MPLYMA 205 : CALL LICHT(IALC , MPALC , TYPE, IGEOM) 206 : C SEGACT MPALC 206 : CALL LICHT(IGRALC , MPGALC , TYPE, IGEOM) 207 : C SEGACT MPGALC 207 : CALL LICHT(ILIALC , MPLALC , TYPE, IGEOM) 208 : C SEGACT MPLALC 209 : ENDIF 210 : C 211 : C**** Le cosinus directeurs 212 : C 212 : CALL LICHT(INORM , MPNORM , TYPE, IGEOM) 213 : C SEGACT MPNORM 214 : C 215 : C**** Les MPOVAL sont déjà activés i.e.: 216 : C 217 : C 218 : C**** Le MELEME FACEL 219 : C 220 : IPT1 = IFACEL 221 : IPT2 = IFACE 222 : SEGACT IPT1 223 : SEGACT IPT2 224 : NFAC = IPT1.NUM(/2) 225 : C 226 : C**** Creation de MCHAMLs contenant les etats gauche et droite, 227 : C 228 : C i.e.: 229 : C 230 : C vitesse + cosinus directors du repere local 231 : C alpha 232 : C densité 233 : C pression 234 : C 235 : C********************************************************** 236 : C**** Cosinus directors du repere local et vitesse ******** 237 : C********************************************************** 238 : C 239 : C Les cosinus directeurs 240 : C 241 : N1 = 2 242 : N3 = 6 243 : L1 = 28 244 : SEGINI MCHEL1 245 : IVN1F = MCHEL1 246 : MCHEL1.TITCHE = 'U ' 247 : MCHEL1.IMACHE(1) = IFACE 248 : MCHEL1.IMACHE(2) = IFACEL 249 : MCHEL1.CONCHE(1) = ' (n,t) in (x,y) ' 250 : MCHEL1.CONCHE(2) = ' U in (n,t) ' 251 : * MCHEL1.NOPHAS(1) = ' ' 252 : * MCHEL1.NOPHAS(2) = ' ' 253 : C 254 : C**** Valeurs des cosinus definies par respect au repair global, i.e. 255 : C 256 : MCHEL1.INFCHE(1,1) = 2 257 : MCHEL1.INFCHE(1,3) = NIFOUR 258 : MCHEL1.INFCHE(1,4) = 0 259 : MCHEL1.INFCHE(1,5) = 0 260 : MCHEL1.INFCHE(1,6) = 0 261 : MCHEL1.IFOCHE = IFOUR 262 : C 263 : C**** Valeurs de vitesse definies par respect au repair local, i.e. 264 : C 265 : MCHEL1.INFCHE(2,1) = 1 266 : MCHEL1.INFCHE(2,3) = NIFOUR 267 : MCHEL1.INFCHE(2,4) = 0 268 : MCHEL1.INFCHE(2,5) = 0 269 : MCHEL1.INFCHE(2,6) = 0 270 : C 271 : C**** Le cosinus directeurs 272 : C 273 : N1PTEL = 1 274 : N1EL = NFAC 275 : N2PTEL = 0 276 : N2EL = 0 277 : C 278 : C**** MCHAML a N2 composantes: 279 : C 280 : C cosinus directeurs du repere local (n,t1) 281 : C 282 : C IDIM = 2 -> 4 composantes 283 : C 284 : N2 = 4 285 : SEGINI MCHAM1 286 : MCHEL1.ICHAML(1) = MCHAM1 287 : MCHAM1.NOMCHE(1) = 'NX ' 288 : MCHAM1.NOMCHE(2) = 'NY ' 289 : MCHAM1.NOMCHE(3) = 'TX ' 290 : MCHAM1.NOMCHE(4) = 'TY ' 291 : MCHAM1.TYPCHE(1) = 'REAL*8 ' 292 : MCHAM1.TYPCHE(2) = 'REAL*8 ' 293 : MCHAM1.TYPCHE(3) = 'REAL*8 ' 294 : MCHAM1.TYPCHE(4) = 'REAL*8 ' 295 : SEGINI MELVNX 296 : SEGINI MELVNY 297 : SEGINI MELT1X 298 : SEGINI MELT1Y 299 : MCHAM1.IELVAL(1) = MELVNX 300 : MCHAM1.IELVAL(2) = MELVNY 301 : MCHAM1.IELVAL(3) = MELT1X 302 : MCHAM1.IELVAL(4) = MELT1Y 303 : SEGDES MCHAM1 304 : C 305 : C**** Vitesse 306 : C 307 : N1EL = NFAC 308 : N1PTEL = 3 309 : N2EL = 0 310 : N2PTEL = 0 311 : C 312 : C**** MCHAML a N2 composantes: 313 : C 314 : C IDIM = 2 -> 2 composantes 315 : C 316 : N2 = 2 317 : SEGINI MCHAM1 318 : MCHEL1.ICHAML(2) = MCHAM1 319 : MCHAM1.NOMCHE(1) = 'UN ' 320 : MCHAM1.NOMCHE(2) = 'UT ' 321 : MCHAM1.TYPCHE(1) = 'REAL*8 ' 322 : MCHAM1.TYPCHE(2) = 'REAL*8 ' 323 : SEGINI MEVUN1 324 : SEGINI MEVUT1 325 : MCHAM1.IELVAL(1) = MEVUN1 326 : MCHAM1.IELVAL(2) = MEVUT1 327 : SEGDES MCHAM1 328 : C 329 : C********************************************************** 330 : C**** PHI1 ******** 331 : C********************************************************** 332 : C 333 : C**** PHI1 334 : C 335 : N1 = 1 336 : N3 = 6 337 : L1 = 15 338 : SEGINI MCHEL2 339 : IPHIF = MCHEL2 340 : MCHEL2.IMACHE(1) = IFACEL 341 : MCHEL2.TITCHE = 'PHI ' 342 : MCHEL2.CONCHE(1) = ' ' 343 : * MCHEL2.NOPHAS(1) = ' ' 344 : C 345 : C**** Valeurs independente du repére, i.e. 346 : C 347 : MCHEL2.INFCHE(1,1) = 0 348 : MCHEL2.INFCHE(1,3) = NIFOUR 349 : MCHEL2.INFCHE(1,4) = 0 350 : MCHEL2.INFCHE(1,5) = 0 351 : MCHEL2.INFCHE(1,6) = 0 352 : MCHEL2.IFOCHE = IFOUR 353 : N2 = 1 354 : SEGINI MCHAM1 355 : MCHEL2.ICHAML(1) = MCHAM1 356 : C We cannot deseactivate MCHEL2 = IPHIF now since we 357 : C use it after... 358 : MCHAM1.NOMCHE(1) = 'SCAL ' 359 : MCHAM1.TYPCHE(1) = 'REAL*8 ' 360 : SEGINI MELPHI 361 : MCHAM1.IELVAL(1) = MELPHI 362 : SEGDES MCHAM1 363 : C 364 : C********************************************************** 365 : C**** IRN1F and IRN2F ******** 366 : C********************************************************** 367 : C 368 : MCHEL1 = IPHIF 369 : SEGINI, MCHEL2 = MCHEL1 370 : IRN1F = MCHEL2 371 : MCHEL2.TITCHE = 'RHO1 ' 372 : MCHAM1 = MCHEL1.ICHAML(1) 373 : SEGINI, MCHAM2 = MCHAM1 374 : MCHEL2.ICHAML(1) = MCHAM2 375 : SEGDES MCHEL2 376 : SEGINI MELRN1 377 : MCHAM2.IELVAL(1) = MELRN1 378 : SEGDES MCHAM2 379 : C 380 : C 381 : C 382 : C********************************************************** 383 : C**** IPN1F 384 : C********************************************************** 385 : C 386 : MCHEL1 = IPHIF 387 : SEGINI, MCHEL2 = MCHEL1 388 : IPN1F = MCHEL2 389 : MCHEL2.TITCHE = 'P1 ' 390 : MCHAM1 = MCHEL1.ICHAML(1) 391 : SEGINI, MCHAM2 = MCHAM1 392 : MCHEL2.ICHAML(1) = MCHAM2 393 : SEGDES MCHEL2 394 : SEGDES MCHEL1 395 : C We desactivate MCHEL1 = IPHIF now ! 396 : SEGINI MELPN1 397 : MCHAM2.IELVAL(1) = MELPN1 398 : SEGDES MCHAM2 399 : C 400 : C write(*,*) 'Qui ci arrivo 1...' 401 : IF (NESP .GE. 1) THEN 402 : C 403 : SEGACT MLMESP 404 : C 405 : C******* YF 406 : C 407 : NS = NESP 408 : SEGINI FRAMAS 409 : MCHEL1 = IRN1F 410 : SEGINI, MCHEL2 = MCHEL1 411 : IYF = MCHEL2 412 : MCHEL2.TITCHE = 'Y ' 413 : N2 = NESP 414 : SEGINI MCHAMY 415 : MCHEL2.ICHAML(1) = MCHAMY 416 : SEGDES MCHEL2 417 : N1EL = NFAC 418 : N1PTEL = 3 419 : N2EL = 0 420 : N2PTEL = 0 421 : DO I1 = 1, NESP 422 : SEGINI MELVA1 423 : MCHAMY.IELVAL(I1) = MELVA1 424 : C AB Error. It is corrected by the following line 425 : CARCEL = MLMESP.MOTS(I1) 426 : TYPE = ' ' 427 : TYPE(1:4) = CARCEL(1:4) 428 : MCHAMY.NOMCHE(I1) = TYPE 429 : MCHAMY.TYPCHE(I1) = 'REAL*8 ' 430 : ENDDO 431 : C 432 : C******* IALF 433 : C 434 : NS = NESP 435 : SEGINI ALPHA 436 : MCHEL1 = IRN1F 437 : SEGINI, MCHEL2 = MCHEL1 438 : IALF = MCHEL2 439 : MCHEL2.TITCHE = 'ALPHA ' 440 : N2 = NESP 441 : SEGINI MCHAMA 442 : MCHEL2.ICHAML(1) = MCHAMA 443 : SEGDES MCHEL2 444 : N1EL = NFAC 445 : N1PTEL = 3 446 : N2EL = 0 447 : N2PTEL = 0 448 : DO I1 = 1, NESP 449 : SEGINI MELVA1 450 : MCHAMA.IELVAL(I1) = MELVA1 451 : C AB Error. It is corrected by the following line 452 : CARCEL = MLMESP.MOTS(I1) 453 : TYPE = ' ' 454 : TYPE(1:4) = CARCEL(1:4) 455 : MCHAMA.NOMCHE(I1) = TYPE 456 : MCHAMA.TYPCHE(I1) = 'REAL*8 ' 457 : ENDDO 458 : C 459 : SEGDES MLMESP 460 : ENDIF 461 : C 462 : C write(*,*) 'Qui ci arrivo 2...' 463 : C 464 : C 465 : C********************************************************** 466 : C**** Boucle sur le faces ********* 467 : C********************************************************** 468 : C 469 : IDIMP1 = IDIM + 1 470 : DO NLCF = 1, NFAC 471 : C 472 : C******* NLCF = numero local du centre de face 473 : C NGCF = numero global du centre de face 474 : C NGCEG = numero global du centre ELT "gauche" 475 : C NLCEG = numero local du centre ELT "gauche" 476 : C NGCED = numero global du centre ELT "droite" 477 : C NLCED = numero local du centre ELT "droite" 478 : C 479 : NGCEG = IPT1.NUM(1,NLCF) 480 : NGCF = IPT1.NUM(2,NLCF) 481 : NGCED = IPT1.NUM(3,NLCF) 482 : NLCEG = MLENT1.LECT(NGCEG) 483 : NLCED = MLENT1.LECT(NGCED) 484 : C 485 : C******* TEST: IPT2.NUM(1,NLCF) = IPT1.NUM(2,NLCF) 486 : C 487 : NGCF1 = IPT2.NUM(1,NLCF) 488 : IF( NGCF1 .NE. NGCF) THEN 489 : MESERR(1:40) = 'PRET, subroutine pre611.eso ' 490 : WRITE(IOIMP,*) MESERR 490 : CALL ERREUR(5) 491 : GOTO 9999 492 : ENDIF 493 : C 494 : C******* Cosinus directeurs des NORMALES aux faces 495 : C 496 : C On impose que les normales sont direct "Gauche" -> "Centre" 497 : C 498 : INDCEL = (NGCEG-1)*IDIMP1 499 : XG = XCOOR(INDCEL+1) 500 : YG = XCOOR(INDCEL+2) 501 : INDCEL = (NGCF-1)*IDIMP1 502 : XC = XCOOR(INDCEL + 1) 503 : YC = XCOOR(INDCEL + 2) 504 : INDCEL = (NGCED-1)*IDIMP1 505 : XD = XCOOR(INDCEL+1) 506 : YD = XCOOR(INDCEL+2) 507 : DXG = XC - XG 508 : DYG = YC - YG 509 : DXD = XC - XD 510 : DYD = YC - YD 511 : C 512 : C******* On calcule le sign du pruduit scalare 513 : C (Normales de Castem) * (vecteur "gauche" -> "centre") 514 : C 515 : CNX = MPNORM.VPOCHA(NLCF,1) 516 : CNY = MPNORM.VPOCHA(NLCF,2) 517 : ORIENT = CNX * DXG + CNY * DYG 518 : ORIENT = SIGN(1.0D0,ORIENT) 519 : IF(ORIENT .NE. 1.0D0)THEN 520 : MESERR(1:30)= 521 : & 'PRET , subroutine pre611.eso. ' 522 : GOTO 9999 523 : ENDIF 524 : CNX = CNX * ORIENT 525 : CNY = CNY * ORIENT 526 : C 527 : C******* Cosinus directeurs de tangent 2D 528 : C 529 : CTX = -1.0D0 * CNY 530 : CTY = CNX 531 : C 532 : C******* Les autres MELVALs 533 : C 534 : C 535 : C******* N.B.: On suppose qu'on a déjà controlle RO, P > 0... 536 : C 537 : C 538 : C******* Etat gauche 539 : C 540 : C PHI 541 : C 542 : VALCEL = MPPHI.VPOCHA(NLCEG, 1) 543 : LIMCEL = MPLPHI.VPOCHA(NLCEG, 1) 544 : DCEL = (MPGPHI.VPOCHA(NLCEG, 1) * DXG) + 545 : & (MPGPHI.VPOCHA(NLCEG, 2) * DYG) 546 : PHIG = VALCEL + (LIMCEL * DCEL) 547 : C write(*,*) valcel, limcel, dcel 548 : C 549 : C 550 : C RN 551 : C 552 : VALCEL = MPRN1.VPOCHA(NLCEG, 1) 553 : LIMCEL = MPLRN1.VPOCHA(NLCEG, 1) 554 : DCEL = (MPGRN1.VPOCHA(NLCEG, 1) * DXG) + 555 : & (MPGRN1.VPOCHA(NLCEG, 2) * DYG) 556 : RN1G = VALCEL + (LIMCEL * DCEL) 557 : C 558 : C 559 : C PN 560 : C 561 : VALCEL = MPPN1.VPOCHA(NLCEG, 1) 562 : LIMCEL = MPLPN1.VPOCHA(NLCEG, 1) 563 : DCEL = (MPGPN1.VPOCHA(NLCEG, 1) * DXG) + 564 : & (MPGPN1.VPOCHA(NLCEG, 2) * DYG) 565 : PN1G = VALCEL + (LIMCEL * DCEL) 566 : C 567 : C VN 568 : C 569 : VALCEL = MPVN1.VPOCHA(NLCEG, 1) 570 : LIMCEL = MPLVN1.VPOCHA(NLCEG, 1) 571 : DCEL = MPGVN1.VPOCHA(NLCEG, 1)*DXG + 572 : & MPGVN1.VPOCHA(NLCEG, 2)*DYG 573 : UX1G = VALCEL + (LIMCEL * DCEL) 574 : C 575 : VALCEL = MPVN1.VPOCHA(NLCEG, 2) 576 : LIMCEL = MPLVN1.VPOCHA(NLCEG, 2) 577 : DCEL = MPGVN1.VPOCHA(NLCEG, 3)*DXG + 578 : & MPGVN1.VPOCHA(NLCEG, 4)*DYG 579 : UY1G = VALCEL + (LIMCEL * DCEL) 580 : C 581 : C YN 582 : C 583 : DO I1 = 1, NESP 584 : INDCEL = 2 * I1 - 1 585 : VALCEL = MPYMA.VPOCHA(NLCEG,I1) 586 : DCEL = MPGYMA.VPOCHA(NLCEG, INDCEL)*DXG + 587 : & MPGYMA.VPOCHA(NLCEG,INDCEL + 1 )*DYG 588 : LIMCEL = MPLYMA.VPOCHA(NLCEG,I1) 589 : FRAMAS.FRAMG(I1) = VALCEL + (LIMCEL * DCEL) 590 : ENDDO 591 : C 592 : C ALPHAN 593 : C 594 : DO I1 = 1, NESP 595 : INDCEL = 2 * I1 - 1 596 : VALCEL = MPALC.VPOCHA(NLCEG,I1) 597 : DCEL = MPGALC.VPOCHA(NLCEG, INDCEL)*DXG + 598 : & MPGALC.VPOCHA(NLCEG,INDCEL + 1 )*DYG 599 : LIMCEL = MPLALC.VPOCHA(NLCEG,I1) 600 : ALPHA.FRAMG(I1) = VALCEL + (LIMCEL * DCEL) 601 : ENDDO 602 : C 603 : C 604 : C******* Si l'on fait pas de prediction, ce n'est pas necessaire de 605 : C controller la positivite' de la pression et de la densité; elle 606 : C est déjà garantie par la proprieté LED de limiteur. 607 : C If we want to check it, just uncomment LOGI1. 608 : C 609 : C LOGI1 = (RN1G .LT. 0.0D0) .OR. 610 : C & (PN1G .LT. 0.0D0) 611 : C 612 : IF ( NGCEG .EQ. NGCED) THEN 613 : C 614 : C********** Cas mur 615 : C 616 : IF(LOGI1)THEN 617 : C 618 : C********** Premier ordre en espace local 619 : C 620 : PHIG = MPPHI.VPOCHA(NLCEG,1) 621 : RN1G = MPRN1.VPOCHA(NLCEG,1) 622 : PN1G = MPPN1.VPOCHA(NLCEG,1) 623 : UX1G = MPVN1.VPOCHA(NLCEG,1) 624 : UY1G = MPVN1.VPOCHA(NLCEG,2) 625 : DO I1 = 1, NESP 626 : FRAMAS.FRAMG(I1) = MPYMA.VPOCHA(NLCEG,I1) 627 : ENDDO 628 : DO I1 = 1, NESP 629 : ALPHA.FRAMG(I1) = MPALC.VPOCHA(NLCEG,I1) 630 : ENDDO 631 : ENDIF 632 : C 633 : UN1G = UX1G * CNX + UY1G * CNY 634 : UT1G = UX1G * CTX + UY1G * CTY 635 : C 636 : C********** Son etat droite 637 : C 638 : PHID = PHIG 639 : RN1D = RN1G 640 : PN1D = PN1G 641 : UN1D = -1.0D0 * UN1G 642 : UT1D = UT1G 643 : DO I1 = 1, NESP 644 : FRAMAS.FRAMD(I1) = FRAMAS.FRAMG(I1) 645 : ENDDO 646 : DO I1 = 1, NESP 647 : ALPHA.FRAMD(I1) = ALPHA.FRAMG(I1) 648 : ENDDO 649 : C 650 : C********** Fin cas mur 651 : C 652 : ELSE 653 : VALCEL = MPPHI.VPOCHA(NLCED, 1) 654 : LIMCEL = MPLPHI.VPOCHA(NLCED, 1) 655 : DCEL = (MPGPHI.VPOCHA(NLCED, 1) * DXD) + 656 : & (MPGPHI.VPOCHA(NLCED, 2) * DYD) 657 : PHID = VALCEL + (LIMCEL * DCEL) 658 : C 659 : C RN 660 : C 661 : VALCEL = MPRN1.VPOCHA(NLCED, 1) 662 : LIMCEL = MPLRN1.VPOCHA(NLCED, 1) 663 : DCEL = (MPGRN1.VPOCHA(NLCED, 1) * DXD) + 664 : & (MPGRN1.VPOCHA(NLCED, 2) * DYD) 665 : RN1D = VALCEL + (LIMCEL * DCEL) 666 : C 667 : C PN 668 : C 669 : VALCEL = MPPN1.VPOCHA(NLCED, 1) 670 : LIMCEL = MPLPN1.VPOCHA(NLCED, 1) 671 : DCEL = (MPGPN1.VPOCHA(NLCED, 1) * DXD) + 672 : & (MPGPN1.VPOCHA(NLCED, 2) * DYD) 673 : PN1D = VALCEL + (LIMCEL * DCEL) 674 : C 675 : C VN 676 : C 677 : VALCEL = MPVN1.VPOCHA(NLCED, 1) 678 : LIMCEL = MPLVN1.VPOCHA(NLCED, 1) 679 : DCEL = MPGVN1.VPOCHA(NLCED, 1)*DXD + 680 : & MPGVN1.VPOCHA(NLCED, 2)*DYD 681 : UX1D = VALCEL + LIMCEL * DCEL 682 : C 683 : VALCEL = MPVN1.VPOCHA(NLCED, 2) 684 : LIMCEL = MPLVN1.VPOCHA(NLCED, 2) 685 : DCEL = MPGVN1.VPOCHA(NLCED, 3)*DXD + 686 : & MPGVN1.VPOCHA(NLCED, 4)*DYD 687 : UY1D = VALCEL + LIMCEL * DCEL 688 : C 689 : C YN 690 : C 691 : DO I1 = 1, NESP 692 : INDCEL = 2 * I1 - 1 693 : VALCEL = MPYMA.VPOCHA(NLCED,I1) 694 : DCEL = MPGYMA.VPOCHA(NLCED, INDCEL)*DXD + 695 : & MPGYMA.VPOCHA(NLCED,INDCEL + 1 )*DYD 696 : LIMCEL = MPLYMA.VPOCHA(NLCED,I1) 697 : FRAMAS.FRAMD(I1) = VALCEL + (LIMCEL * DCEL) 698 : ENDDO 699 : C 700 : C ALPHAN 701 : C 702 : DO I1 = 1, NESP 703 : INDCEL = 2 * I1 - 1 704 : VALCEL = MPALC.VPOCHA(NLCED,I1) 705 : DCEL = MPGALC.VPOCHA(NLCED, INDCEL)*DXD + 706 : & MPGALC.VPOCHA(NLCED,INDCEL + 1 )*DYD 707 : LIMCEL = MPLALC.VPOCHA(NLCED,I1) 708 : ALPHA.FRAMD(I1) = VALCEL + (LIMCEL * DCEL) 709 : ENDDO 710 : 711 : C 712 : C 713 : C********** Si l'on fait pas de prediction, ce n'est pas necessaire de 714 : C controller la positivite' de la pression et de la densité; elle 715 : C est déjà garantie par la proprieté LED de limiteur. 716 : C If we want to check it, just uncomment LOGI1. 717 : C 718 : C LOGI1 = LOGI1 .OR. (RN1D .LT. 0.0D0) 719 : C $ .OR.(PN1D .LT. 0.0D0) 720 : C 721 : IF(LOGI1)THEN 722 : C 723 : C************* Premier ordre en espace local 724 : C 725 : PHIG = MPPHI.VPOCHA(NLCEG,1) 726 : RN1G = MPRN1.VPOCHA(NLCEG,1) 727 : PN1G = MPPN1.VPOCHA(NLCEG,1) 728 : UX1G = MPVN1.VPOCHA(NLCEG,1) 729 : UY1G = MPVN1.VPOCHA(NLCEG,2) 730 : DO I1 = 1, NESP 731 : FRAMAS.FRAMG(I1) = MPYMA.VPOCHA(NLCEG,I1) 732 : ENDDO 733 : DO I1 = 1, NESP 734 : ALPHA.FRAMG(I1) = MPALC.VPOCHA(NLCEG,I1) 735 : ENDDO 736 : C 737 : PHID = MPPHI.VPOCHA(NLCED,1) 738 : RN1D = MPRN1.VPOCHA(NLCED,1) 739 : PN1D = MPPN1.VPOCHA(NLCED,1) 740 : UX1D = MPVN1.VPOCHA(NLCED,1) 741 : UY1D = MPVN1.VPOCHA(NLCED,2) 742 : DO I1 = 1, NESP 743 : FRAMAS.FRAMD(I1) = MPYMA.VPOCHA(NLCED,I1) 744 : ENDDO 745 : DO I1 = 1, NESP 746 : ALPHA.FRAMD(I1) = MPALC.VPOCHA(NLCED,I1) 747 : ENDDO 748 : ENDIF 749 : C 750 : UN1G = UX1G * CNX + UY1G * CNY 751 : UT1G = UX1G * CTX + UY1G * CTY 752 : C 753 : UN1D = UX1D * CNX + UY1D * CNY 754 : UT1D = UX1D * CTX + UY1D * CTY 755 : C 756 : ENDIF 757 : C 758 : C 759 : C******** Les MELVALs 760 : C 761 : MELPHI.VELCHE(1,NLCF) = PHIG 762 : MELPHI.VELCHE(3,NLCF) = PHID 763 : C 764 : MELRN1.VELCHE(1,NLCF) = RN1G 765 : MELRN1.VELCHE(3,NLCF) = RN1D 766 : C 767 : MELPN1.VELCHE(1,NLCF) = PN1G 768 : MELPN1.VELCHE(3,NLCF) = PN1D 769 : C 770 : MEVUN1.VELCHE(1,NLCF) = UN1G 771 : MEVUN1.VELCHE(3,NLCF) = UN1D 772 : MEVUT1.VELCHE(1,NLCF) = UT1G 773 : MEVUT1.VELCHE(3,NLCF) = UT1D 774 : MELVNX.VELCHE(1,NLCF) = CNX 775 : MELVNY.VELCHE(1,NLCF) = CNY 776 : MELT1X.VELCHE(1,NLCF) = CTX 777 : MELT1Y.VELCHE(1,NLCF) = CTY 778 : C 779 : DO I1 = 1, NESP 780 : MELVA1 = MCHAMY.IELVAL(I1) 781 : MELVA1.VELCHE(1,NLCF) = FRAMAS.FRAMG(I1) 782 : MELVA1.VELCHE(3,NLCF) = FRAMAS.FRAMD(I1) 783 : ENDDO 784 : C 785 : DO I1 = 1, NESP 786 : MELVA1 = MCHAMA.IELVAL(I1) 787 : MELVA1.VELCHE(1,NLCF) = ALPHA.FRAMG(I1) 788 : MELVA1.VELCHE(3,NLCF) = ALPHA.FRAMD(I1) 789 : ENDDO 790 : C 791 : ENDDO 792 : C 793 : C**** Desactivation des SEGMENTs 794 : C 795 : SEGDES IPT1 796 : SEGDES IPT2 797 : C 798 : C MPOVALs 799 : C 800 : SEGDES MPNORM 801 : C 802 : SEGDES MPPHI 803 : SEGDES MPGPHI 804 : SEGDES MPLPHI 805 : C 806 : SEGDES MPRN1 807 : SEGDES MPGRN1 808 : SEGDES MPLRN1 809 : C 810 : SEGDES MPPN1 811 : SEGDES MPGPN1 812 : SEGDES MPLPN1 813 : C 814 : SEGDES MPVN1 815 : SEGDES MPGVN1 816 : SEGDES MPLVN1 817 : C 818 : C MELVALs 819 : C 820 : SEGDES MELVNX 821 : SEGDES MELVNY 822 : SEGDES MELT1X 823 : SEGDES MELT1Y 824 : SEGDES MEVUN1 825 : SEGDES MEVUT1 826 : C 827 : SEGDES MELPHI 828 : C 829 : SEGDES MELRN1 830 : C 831 : SEGDES MELPN1 832 : C 833 : IF (NESP .GE. 1)THEN 834 : SEGDES MPYMA 835 : SEGDES MPGYMA 836 : SEGDES MPLYMA 837 : SEGDES MPALC 838 : SEGDES MPGALC 839 : SEGDES MPLALC 840 : DO I1 = 1, NESP 841 : MELVA1 = MCHAMY.IELVAL(I1) 842 : SEGDES MELVA1 843 : MELVA1 = MCHAMA.IELVAL(I1) 844 : SEGDES MELVA1 845 : ENDDO 846 : SEGDES MCHAMA 847 : SEGDES MCHAMY 848 : SEGSUP FRAMAS 849 : SEGSUP ALPHA 850 : ENDIF 851 : C 852 : C**** Destruction du MELNTI correspondance local/global 853 : C 854 : SEGSUP MLENT1 855 : C 856 : 9999 CONTINUE 857 : C 858 : RETURN 859 : END 860 : 861 : 862 : 863 : 864 : 865 : 866 : 867 : 868 : 869 : 870 : 871 : 872 :
© Cast3M 2003 - All rights reserved.
Disclaimer