1 : * UPDAEFMH PROCEDUR CHAT 12/08/07 21:15:14 7481 2 : ********************************************************************** 3 : 'DEBP' UPDAEFMH MoDARCY*'MMODEL' Porosite*'CHPOINT' MateDiff*'CHPOINT' 4 : Diffdisp*'CHPOINT' ChPSour*'CHPOINT' Cini*'CHPOINT' 5 : Tcini*'CHPOINT' Deltat*'FLOTTANT' Qface/'CHPOINT' 6 : nomespec*'LISTMOTS' nbespece*'ENTIER' nbsource*'ENTIER' 7 : LMLump*'LOGIQUE' DECENTR*'LOGIQUE' massEFMH*'RIGIDITE' 8 : mattr/'MATRIK' mattm/'RIGIDITE' Tbdartra*'TABLE' 9 : TABMODI*'TABLE' CHCLIM*'TABLE'; 10 : * ATTENTION La vitesse est optionnelle, L'ordre est important 11 : * et les types d'arguments qui se suivent aussi pour tester leur 12 : * présence 13 : * 14 : * Attention il faudra transformer les vitesses en débits face 15 : * et sortir le champ 16 : * 17 : * |-----------------------------------------------------------------| 18 : * | Phrase d'appel (en GIBIANE) | 19 : * |-----------------------------------------------------------------| 20 : * | | 21 : * |SMTr MatrTr TbDarTra MassEFMH Difftot = UPDAEFMH MoDARCY Porosite| 22 : * | MateDiff difftot ChPSour cini tcini deltat | 23 : * | (Qface) nomespec nbespece nbsource LMLump | 24 : * | DECENTR massEFMH | mattr tbdartra TABMODI; | 25 : * | | mattm | 26 : * |-----------------------------------------------------------------| 27 : * | Généralités : MATTEFMH construit la matrice de discrétisation | 28 : * | du problème de transport convection-diffusion pour| 29 : * | le premier pas de tps d'un algorithme transitoire.| 30 : * | Le second membre et les Conditions limites de flux| 31 : * | sont pris en compte. | 32 : * | RESTE TCINI, DECENTR et TERME LIN | 33 : * |-----------------------------------------------------------------| 34 : * | | 35 : * |-----------------------------------------------------------------| 36 : * | ENTREES | 37 : * |-----------------------------------------------------------------| 38 : * | MoDARCY : modele Darcy. | 39 : * | | 40 : * | Porosite : champ par elements de composante 'CK' | 41 : * | | 42 : * | MateDiff : Tenseur de diffusion (type iso, ..) champ par | 43 : * | points de composante 'K' en isotrope, 'K11', 'K21', | 44 : * | 'K22' en anisotrope 2d et 'K11', 'K21', 'K22', 'K31'| 45 : * | 'K32', 'K33' en anisotrope 3d. Type 'CARACTERISTIQUE'| 46 : * | | 47 : * | Diffdisp : Tenseur de dispersion (type iso, ..) champ par | 48 : * | points de composante 'K' en isotrope, 'K11', 'K21', | 49 : * | 'K22' en anisotrope 2d et 'K11', 'K21', 'K22', 'K31'| 50 : * | 'K32', 'K33' en anisotrope 3d. Type 'CARACTERISTIQUE'| 51 : * | | 52 : * | ChPSour : Champ par points des sources volumiques par unité de | 53 : * | temps (support maillage centre). Composante ?????? | 54 : * | | 55 : * | Cini : Concentration initiale, CHPOINT centre. | 56 : * | Composante 'H'. | 57 : * | | 58 : * | Tcini : Trace de concentration aux faces (eventuellement à | 59 : * | plusieurs composantes (espèces) | 60 : * | | 61 : * | Deltat : Pas de temps | 62 : * | | 63 : * | Qface : vitesse aux faces, CHPO face de composantes Vx, Vy | 64 : * | en 2d et Vx, Vy, Vz en 3d. Il s'agit plus exatement | 65 : * | de (V.n)n, c'est à dire de la composante normale de | 66 : * | la vitesse aux faces. ???????? (je pressens que | 67 : * | castem va sortir des flux, cad intégrés sur surfaces)| 68 : * | | 69 : * | nomespec : liste des noms de composante des espèces dans Cini | 70 : * | | 71 : * | nbespece : nombre de composante de Cini, soit nombre d'especes | 72 : * | | 73 : * | nbsource : nombre de composantes du terme source qd X especes | 74 : * | | 75 : * | LMLump : Logique. Si vrai on effectue une condensation de | 76 : * | masse de la matrice EFMH | 77 : * | | 78 : * | DECENTR : Logique. Vrai veut dire schémas décentrés et faux | 79 : * | veut dire schéma convectif centré. | 80 : * | | 81 : * | MatTm : matrice globale sur les traces. MATRIK en entrée | 82 : * | sort MATRIK si non modifiée, RIGIDITE sinon | 83 : * | Soit on rentre cet argument soit le suivant Mattr | 84 : * | | 85 : * | MatTr : idem mais rigidité en entrée on ressort cette matrice| 86 : * | inchangée si les options MATMODI indiquent aucune | 87 : * | modif. Optionnel. On rentre Mattm si absent. | 88 : * | | 89 : * | TbDarTra : table Darcy transitoire utilisée par MHYB, SMTP ... | 90 : * | | 91 : * | TABMODI : table contenant des logiques indiquant la nécessité | 92 : * | ou non de reclalculer certains termes. | 93 : * | 'POROSITE' : VRAI si le coefficient devant D/DT | 94 : * | (porosité) est modifié depuis le dernier| 95 : * | appel | 96 : * | 'DELTAT' : VRAI si le pas de tps a changé | 97 : * | 'CONVECTI' : VRAI si la vitesse a changé | 98 : * | 'COEF_LIN' : VRAI si le coef en facteur de C a changé| 99 : * | 'DIFFUSI' : VRAI si les diffusivités ont changé | 100 : * | | 101 : * | CHCLIM : table d'indice 'NEUMANN' et 'DIRICHLET' contenant les| 102 : * | Chpoint à n composantes contenant les conditions aux | 103 : * | limites de Neumann et Dirichlet par espece. | 104 : * | | 105 : * | | 106 : * |-----------------------------------------------------------------| 107 : * | ENTREES-SORTIES | 108 : * |-----------------------------------------------------------------| 109 : * | | 110 : * | MassEFMH : matrice elementaire EFMH | 111 : * | | 112 : * | Remarque | 113 : * | -------- | 114 : * | On a toujours intéret à rentrer Mattm si on l'a et qu'il n'y a | 115 : * | pas de modification, afin de conserver les factorisations LU | 116 : * | dans solvEFMH en cas de résolution par méthode directe. | 117 : * | | 118 : * |-----------------------------------------------------------------| 119 : * | SORTIES | 120 : * |-----------------------------------------------------------------| 121 : * | | 122 : * | SMTr : second membre sur les traces | 123 : * | | 124 : * | Matrtr : Matrice du problème, rigidité si nouvelle, égale à la| 125 : * | matrice entrée (mattm ou mattr) sinon | 126 : * | | 127 : * |-----------------------------------------------------------------| 128 : * | VARIABLES INTERNES | 129 : * |-----------------------------------------------------------------| 130 : * | | 131 : * | LCALMATP : Logique, VRAI si on recalcule la matrice du systeme | 132 : * | avec diffusion produite par MATP | 133 : * | | 134 : * | LCALCONV : Logique, VRAI si on recalcule la matrice du systeme | 135 : * | avec convection produite par SMTP | 136 : * | | 137 : * | PCONV : Logique indiquant VRAI si présence de convection | 138 : * | | 139 : * | CoefDt : coeff devant dC/dt integre sur les elements | 140 : * | | 141 : * | Numdiff : diffusivité numérique due au décentrement | 142 : * | | 143 : * | CCini : concentration aux centres (une composante) | 144 : * | | 145 : * | Tccini : Trace de concentration aux faces (une composante) | 146 : * | | 147 : * | SSource : Source aux centre (une composante) | 148 : * | | 149 : * | SSMTr : second membre sur les traces pour une espèce | 150 : * | | 151 : * | MatConv : matrice globale sur les traces pour la convection | 152 : * | | 153 : * | Matrdiff : matrice globale sur les traces pour la diffusion | 154 : * | | 155 : * | toltheta : 1.D-4 seuil en dessous duquel on considère que la | 156 : * | valeur de theta du theta-schéma est nulle (schéma | 157 : * | explicite) ou au contraire euler-implicite si | 158 : * | theta > 0.9999 | 159 : * | | 160 : * | tetaconv : valeur du coef d'implicitation pour la convection | 161 : * | | 162 : * | MODDIFF : Logique, si VRAI on recalcule les termes liés à la | 163 : * | diffusion, en supposant qu'elle a changé, donnés par | 164 : * | Matediff. En implicite cela affecte la matrice | 165 : * | rigidité du problème. | 166 : * | | 167 : * | MODPOROS : Si VRAI on recalcule la contribution de la porosité | 168 : * | pour la matrice du problème et le coeficient devant | 169 : * | dC/dt (EFMH) intégré sur les mailles. la matrice | 170 : * | Identité des VF est recalculée. Logique. | 171 : * | | 172 : * | MODELTT : Logique. Si vrai on recalcule la matrice rigidité du | 173 : * | problème. Le pas de temps a changé. | 174 : * | | 175 : * | MODCONV : Logique. Si vrai on recalcule la matrice rigidité du | 176 : * | problème. La vitesse a changé. | 177 : * | | 178 : * | MODCOLIN : Logique. Si vrai on recalcule la matrice rigidité du | 179 : * | problème. Le coef linéaire devant C a changé. | 180 : * | | 181 : * | FLUNEU : LOGIQUE valant VRAI si conditions de Neumann | 182 : * | | 183 : * | CLFLUX : Chpoint à n composantes contenant les flux imposés | 184 : * | pour chaque espece chimique. nul si pas de flux | 185 : * | OPTIONNEL | 186 : * | | 187 : ********************************************************************** 188 : 189 : 190 : *--------------------------------------------------------------------- 191 : *---------- On récupere les conditions limites ------------------ 192 : *--------------------------------------------------------------------- 193 : FLUNEU = FAUX ; 194 : FLUTOT = FAUX ; 195 : FLUMIX = FAUX ; 196 : FLUCLI = FAUX ; 197 : 198 : * Neumann 199 : 'SI' ('EXISTE' CHCLIM 'NEUMANN') ; 200 : CLFLUX1 = CHCLIM . 'NEUMANN' ; 201 : FLUNEU = VRAI ; 202 : 'SINON' ; 203 : * on crée un champ vide 204 : CLFLUX1 dum = 'KOPS' MATRIK ; 205 : 'OUBLIER' dum ; 206 : 'FINSI' ; 207 : 208 : * Flux total 209 : 'SI' ('EXISTE' CHCLIM 'FLUTOTAL') ; 210 : CLFLUX2 = CHCLIM . 'FLUTOTAL' ; 211 : FLUTOT = VRAI ; 212 : 'SINON' ; 213 : * on crée un champ vide 214 : CLFLUX2 dum = 'KOPS' MATRIK ; 215 : 'OUBLIER' dum ; 216 : 'FINSI' ; 217 : 218 : * Flux mixte 219 : 'SI' ('EXISTE' CHCLIM 'FLUMIXTE') ; 220 : * comme on impose A Dgrad C + B C = flumix, on le traite sous 221 : * la forme D grad C + (B/A) C = flumix/A plus naturelle en EFMH car 222 : * D grad C est le flux diffusif 223 : COFA = -1.D0 * CHCLIM . 'FLUMIXTE' . 'COEFA' ; 224 : CLFLUX3 = CHCLIM . 'FLUMIXTE' . 'VAL' '/' COFA ; 225 : CLFLUX3 = CHAN 'ATTRIBUT' CLFLUX3 NATURE DISCRET ; 226 : FLUMIX = VRAI ; 227 : 'SINON' ; 228 : * on crée un champ vide 229 : CLFLUX3 dum = 'KOPS' MATRIK ; 230 : 'OUBLIER' dum ; 231 : 'FINSI' ; 232 : 233 : 234 : * On fabrique le terme de flux complet 235 : 'SI' (FLUNEU 'OU' FLUTOT 'OU' FLUMIX) ; 236 : CLFLUX = CLFLUX1 'ET' CLFLUX2 'ET' CLFLUX3 ; 237 : FLUCLI = VRAI ; 238 : 'FINSI' ; 239 : 240 : 241 : 242 : *--------------------------------------------------------------------- 243 : *---------- Initialisations de tables, coefficients ------------------ 244 : *--------------------------------------------------------------------- 245 : 246 : toltheta = 1.D-4 ; 247 : 248 : MODPOROS = TABMODI . 'POROSITE' ; 249 : MODCONV = TABMODI . 'CONVECTI' ; 250 : MODDELTT = TABMODI . 'DELTAT' ; 251 : MODCOLIN = TABMODI . 'COEF_LIN' ; 252 : MODDIFF = TABMODI . 'DIFFUSIV' ; 253 : 254 : * On ne recalcule pas a priori la matrice du probleme 255 : LCALMATP = FAUX ; 256 : * On ne recalculera pas a priori la matrice de convection du probleme 257 : LCALCONV = FAUX ; 258 : 259 : 'SI' ('EXISTE' QFACE) ; 260 : PCONV = VRAI ; 261 : 'SINON' ; 262 : PCONV = FAUX ; 263 : 'FINSI' ; 264 : 265 : 'SI' (PCONV) ; 266 : tetaconv = TbDarTra . 'THETA_CONVECTION' ; 267 : 'FINSI' ; 268 : 269 : 270 : * Si la porosité est modifiée (ou plus exactement le coef 271 : * devant dC/dt) 272 : 273 : 274 : 'SI' (MODPOROS) ; 275 : * Calcul du terme devant le dC/dt integré sur le volume 276 : * CHAMELEM 'SCAL'. 277 : CoefDt = ('NOMC' 'SCAL' Porosite) ; 278 : COefDt = Porosite * ('DOMA' modarcy VOLUME) ; 279 : CoefDt = 'KCHA' Modarcy Coefdt 'CHAM' ; 280 : TbDarTra . 'SURF' = CoefDt ; 281 : * on recalculera la matrice totale du probleme 282 : LCALMATP = VRAI ; 283 : 'FINSI' ; 284 : 285 : * 286 : * Si la convection est modifiée (vitesse modifiée, voir apparition 287 : * ou suppression d'un champ de vitesse) en cours de simulation. 288 : * 289 : 290 : 'SI' (PCONV 'ET' MODCONV) ; 291 : * TbDarTra . 'THETA_CONVECTION'= TetaConv ; 292 : * On recalculera la matrice de convection du probleme 293 : LCALCONV = VRAI ; 294 : 'FINSI' ; 295 : 296 : * 297 : * Si le pas de temps est modifié 298 : * 299 : 300 : 'SI' (MODDELTT) ; 301 : TbDarTra . 'PAS' = DeltaT ; 302 : * on recalculera la matrice totale du probleme 303 : LCALMATP = VRAI ; 304 : LCALCONV = VRAI ; 305 : 'FINSI' ; 306 : 307 : * 308 : * Si le coefficient de décroissance (ou linaire) en facteur de C est 309 : * modifié ????????????????, 310 : * 311 : 312 : 'SI' (MODCOLIN) ; 313 : * A priori sera géré dans le coeff de porosité 314 : 'FINSI' ; 315 : 316 : 317 : *--------------------------------------------------------------------- 318 : *--------------------- CALCUL DU DECENTREMENT ------------------------ 319 : *--------------------------------------------------------------------- 320 : 321 : compmat = 'EXTRAIRE' Modarcy 'MATERIAU' ; 322 : ANISO = ('DIME' compmat) > 1 ; 323 : 324 : * Seulement si présence de convection et option décentrée 325 : * Modifié seulement si la convection ou la diffusion changent (Péclet 326 : * de maille = u dx / diff) 327 : 328 : 'SI' (PCONV 'ET' MODCONV) ; 329 : * la convection change d'où recalcule de la matrice du problème 330 : * plus la matrice de convection. La matrice de diffusion numérique et 331 : * de dispersion change aussi. 332 : LCALCONV = VRAI ; 333 : LCALMATP = VRAI ; 334 : * les dispersivités changent d'où recalcule de matrice diffusion 335 : * plus la matrice du problème plus la matrice de convection 336 : 337 : * GBM IL RESTE A CODER LE RECALCUL DE LA DISPERSION (IDEM EN VF) 338 : * ET DU DECENTREMENT (fait en VF). 339 : 340 : 'FINSI' ; 341 : 342 : 'SI' (MODDIFF) ; 343 : * recalcul de la matrice complete du probleme + de la diffusion 344 : LCALMATP = VRAI ; 345 : zozo = DIFFANIS Matediff 'EFMH' ANISO ; 346 : Matediff = zozo ; 347 : 'FINSI' ; 348 : 349 : * Difftot calculé que si MODDIF ou dispersion change (a faire alors) 350 : 351 : SI (MODDIFF) ; 352 : Difftot = MateDiff '+' diffdisp ; 353 : difftot = 'KCHA' Modarcy difftot 'CHAM' ; 354 : difftot = CHANGER difftot TYPE CARACTERISTIQUES ; 355 : 356 : SI (ANISO) ; 357 : 'SI' ('EGA' ('VALEUR' 'DIME') 2) ; 358 : difftot = 'MATE' Modarcy 'DIRECTION' 359 : (1. 0. ) 'PARALLELE' 'K11' ('EXCO' 'K11' difftot) 360 : 'K21' ('EXCO' 'K21' difftot) 361 : 'K22' ('EXCO' 'K22' difftot) ; 362 : 'FINSI' ; 363 : 'SI' ('EGA' ('VALEUR' 'DIME') 3) ; 364 : difftot = 'MATE' Modarcy 'DIRECTION' 365 : (1. 0. 0.) (0. 1. 0.) 'PARALLELE' 366 : 'K11' ('EXCO' 'K11' difftot) 367 : 'K21' ('EXCO' 'K21' difftot) 368 : 'K22' ('EXCO' 'K22' difftot) 369 : 'K31' ('EXCO' 'K31' difftot) 370 : 'K32' ('EXCO' 'K32' difftot) 371 : 'K33' ('EXCO' 'K33' difftot) ; 372 : 'FINSI' ; 373 : SINON ; 374 : difftot = 'MATE' Modarcy 'K' difftot ; 375 : FINSI ; 376 : FINSI ; 377 : 378 : 379 : * 'SI' ('EGA' ('VALEUR' 'DIME') 2) ; 380 : * difftot = 'MATE' Modarcy 'DIRECTION' 381 : * (1. 0. ) 'PARALLELE' 'K11' ('EXCO' 'K11' difftot) 382 : * 'K21' ('EXCO' 'K21' difftot) 383 : * 'K22' ('EXCO' 'K22' difftot) ; 384 : 385 : * 'FINSI' ; 386 : * 'SI' ('EGA' ('VALEUR' 'DIME') 3) ; 387 : * difftot = 'MATE' Modarcy 'DIRECTION' 388 : * (1. 0. 0.) (0. 1. 0.) 'PARALLELE' 389 : * 'K11' ('EXCO' 'K11' difftot) 390 : * 'K21' ('EXCO' 'K21' difftot) 391 : * 'K22' ('EXCO' 'K22' difftot) 392 : * 'K31' ('EXCO' 'K31' difftot) 393 : * 'K32' ('EXCO' 'K32' difftot) 394 : * 'K33' ('EXCO' 'K33' difftot) ; 395 : * 'FINSI' ; 396 : *'FINSI' ; 397 : 398 : *--------------------------------------------------------------------- 399 : *-------------- Matrice masse inverse des EFMH ----------------------- 400 : *--------------------------------------------------------------------- 401 : 402 : 403 : * Calcul des matrices de masse elementaires inverses 404 : * affectées par un changement de la diffusivité totale. 405 : * Cela ce porduit si la diffusivité change ou, en présence 406 : * de décentrement, si la diffusivité numérique change. 407 : 'SI' (MODDIFF 'OU' (DECENTR 'ET' PCONV 'ET' MODCONV)) ; 408 : 'SI' (LMLump) ; 409 : * masse lumping 410 : MassEFMH = 'MHYB' MoDARCY Difftot 'LUMP' ; 411 : 'SINON' ; 412 : MassEFMH = 'MHYB' MoDARCY Difftot ; 413 : 'FINSI' ; 414 : LCALMATP = VRAI ; 415 : LCALCONV = VRAI ; 416 : 'FINSI' ; 417 : 418 : 419 : *---------------------------------------------------------------------- 420 : *---------- CALCUL DE LA MATRICE DE DIFFUSION DU PROBL ---------------- 421 : *---------------------------------------------------------------------- 422 : 423 : 424 : * Calcul de la matrice du probleme diffusion transitoire 425 : * Elle est modifiée si la matrice mass EFMH change, 426 : * ou si le pas de temps change, ou si le coefficient 427 : * en facteur de la dérivée temporelle change (porosité*retard), 428 : * ou si le coefficient du terme linéaire a(x,t) * C(x,t) change. 429 : 430 : * Pour des raisons pratiques on n'a pas stocké les matrices 431 : * de convection et de diffusion séparément. Si l'une est modifiée 432 : * on doit donc recalculer les deux matrices. Ce choix est discutable, 433 : * on gagne un facteur 3 en stockage (mat totale au lieu de mat totale 434 : * + diffusion + convection). Par ailleurs, le cout de l'assemblage + 435 : * resolution + créer une des deux matrice rend souvent négligeable 436 : * le recalcul inutile d'une des deux matrices précitées. On réduit 437 : * également le nombre d'entrées-sorties. 438 : 439 : 'SI' (LCALMATP 'OU' LCALCONV) ; 440 : Matrdiff = 'MATP' MoDARCY MassEFMH TbDarTra ; 441 : 'FINSI' ; 442 : 443 : 444 : 445 : *---------------------------------------------------------------------- 446 : *--- CALCUL DE LA MATRICE DE CONVECTION ET DES SECONDS MEMBRES -------- 447 : *---------------------------------------------------------------------- 448 : 449 : 450 : * Nombres d'espèces à gérer est fixé et ne peut changer au cours d'une 451 : * simulation. Il sera d'ailleurs mis dans le modele 452 : 453 : 'SI' (nbespece > 0) ; 454 : 455 : * pour un schéma en temps non euler implicite, il faut 456 : * la trace à l'instant précédent pour le second membre 457 : * pour la convection ou la diffusion. Elle est issue 458 : * des calculs précédents en EFMH. C'est impératif. 459 : * On ne prévoit pas de la recalculer car c'est très cher. 460 : 461 : 'REPETER' bloc1 nbespece ; 462 : * On extrait la composante de Cini, Tcini et de la source 463 : CCini = 'NOMC' 'H' ('EXCO' (extr &bloc1 nomespec) Cini) ; 464 : TCCini = 'NOMC' 'TH' ('EXCO' (extr &bloc1 nomespec) Tcini) ; 465 : 'SI' (nbsource > 1) ; 466 : SSource = 'NOMC' 'SOUR' ('EXCO' (extr &bloc1 nomespec) 467 : Chpsour) ; 468 : 'SINON' ; 469 : SSource = 'NOMC' 'SOUR' Chpsour ; 470 : 'FINSI' ; 471 : * Conditions initiales 472 : TbDarTra . 'CHARGE' = CCini ; 473 : TbDarTra . 'TRACE'= TCCini ; 474 : * Prise en compte du terme source et eventuellement 475 : * de la convection avec le schema centre 476 : 'SI' (PCONV) ; 477 : * convection 478 : 'SI' (TetaConv 'NEG' 0.0D0 toltheta) ; 479 : * schéma partiellement implicite, matrice MatConv 480 : * retirer pour voir Matconv à gauche si non nécessaire 481 : * ???????????? (avec test sur MODCCONV) 482 : MatConv SSMTr = 'SMTP' MoDARCY MassEFMH TbDarTra SSource 483 : (NOMC 'FLUX' QFace) ; 484 : 'SINON' ; 485 : * schéma explicite, calcul du second membre uniquement 486 : SSMTr = 'SMTP' MoDARCY MassEFMH TbDarTra SSource 487 : (NOMC 'FLUX' QFace) ; 488 : 'FINSI' ; 489 : 'SINON' ; 490 : * pas de convection, calcul du second membre restant 491 : SSMTr = 'SMTP' MoDARCY MassEFMH TbDarTra SSource ; 492 : 'FINSI' ; 493 : * On reconstitue un champ de second membre 494 : 'SI' ((&bloc1) 'EGA' 1) ; 495 : SSMTr = 'NOMC' ('EXTRAIRE' &bloc1 nomespec) SSMTr ; 496 : SMTR = SSMTr ; 497 : 'SINON' ; 498 : SSMTr = 'NOMC' ('EXTRAIRE' &bloc1 nomespec) ('COPIER' SSMTr) ; 499 : SMTr = SSMTr 'ET' SMTr; 500 : 'FINSI' ; 501 : 'FIN' bloc1 ; 502 : 503 : * on stoque la matrice en assemblant la matrice de convection 504 : * du dernier calcul (elles sont toutes identiques). 505 : 506 : 'SI' (PCONV 'ET' (LCALMATP 'OU' LCALCONV)) ; 507 : 'SI' (TetaConv 'NEG' 0.0D0 toltheta) ; 508 : * calcul partiellement implicite pour la convection 509 : MatrTr = Matrdiff 'ET' MatConv ; 510 : 'DETRUIT' MatConv ; 511 : 'DETRUIT' Matrdiff ; 512 : 'SINON' ; 513 : MatrTr = Matrdiff ; 514 : 'FINSI' ; 515 : 'SINON' ; 516 : 'SI' (PCONV 'ET' (TetaConv 'NEG' 0.0D0 toltheta)) ; 517 : 'DETRUIT' MatConv ; 518 : 'FINSI' ; 519 : 'SI' (('NON' PCONV) 'ET' LCALMATP) ; 520 : MatrTr = Matrdiff ; 521 : 'FINSI' ; 522 : 'FINSI' ; 523 : 524 : 'FINSI' ; 525 : 526 : 527 : *--------------------------------------------------------------------- 528 : *------ Conditions aux limites de Flux (mixtes et Neumann) ----------- 529 : *--------------------------------------------------------------------- 530 : 531 : 532 : 'SI' (FLUCLI) ; 533 : SMTr = SMTR 'ET' CLFLUX ; 534 : 'FINSI' ; 535 : 536 : 537 : *--------------------------------------------------------------------- 538 : *------ Matrice sortante --------------------------------- ----------- 539 : *--------------------------------------------------------------------- 540 : 541 : * si pas de modif et matrice rigidité en entrée 542 : 'SI' (('EXISTE' mattr) 'ET' ('NON' (LCALMATP 'OU' LCALCONV))) ; 543 : Matrtr = mattr; 544 : 'FINSI' ; 545 : 546 : * si pas de modif et matrice matrik en entrée 547 : 'SI' (('EXISTE' mattm) 'ET' ('NON' (LCALMATP 'OU' LCALCONV))) ; 548 : Matrtr = mattm; 549 : 'FINSI' ; 550 : 551 : 552 : 'FINP' SMTr MatrTr TbDarTra MassEFMH ; 553 : 554 : 555 : 556 :
© Cast3M 2003 - All rights reserved.
Disclaimer