1 : C HDEBIT SOURCE CHAT 11/03/16 21:24:42 6902 2 : SUBROUTINE HDEBIT 3 : C----------------------------------------------------------------------- 4 : C Calcul du débit à travers chaque face dans le cas d'une 5 : C formulation DARCY en éléments finis mixte hybride. 6 : C----------------------------------------------------------------------- 7 : C 8 : C Une face appartient en général à deux éléments. On ne définit qu'une 9 : C normale à la face. Le ddl aux faces est VN * S * signe(normale) 10 : C avec signe=1 si normale sortante et -1 sinon, S=Surface de la face et 11 : C VN=produit scalaire entre la vitesse à la face et la normale sortante. 12 : C 13 : C On obtient le débit en effectuant le produit du champoint de trace de 14 : C charge TH par une matrice - exprimée en fonction de la matrice de 15 : C rigidité de sous type DARCY et de la matrice DIV (matrice uniligne 16 : C composé de 1) - d'une part, par le MCHAML des orientations de normale 17 : C d'autre part. 18 : C 19 : C Comme nous utilisons des élémennts finis non conforme, il ne faut pas 20 : C assembler les quantitées calculées. Comme le flux est continu au signe 21 : C près, on ne stocke que le flux dans le sens de la normale calculée. 22 : C 23 : C -1 t (n) -1 (n) 24 : C Q = M1 * DIV * H - M1 * TH 25 : C 26 : C 27 : C avec M1 : Matrice massse hybride 28 : C DIV : Matrice uniligne representant la divergence 29 : C 30 : C 31 : C--------------------------- 32 : C Phrase d'appel (GIBIANE) : 33 : C--------------------------- 34 : C 35 : C CHP1 = HDEB MMODEL MRIGI1 CHP2 CHP3 (MRIGI2 CHPO4) (CHP5); 36 : C 37 : C------------------------ 38 : C Opérandes et résultat : 39 : C------------------------ 40 : C 41 : C CHP1 : CHAMPOINT contenant le débit porté par la normale à la face 42 : C MMODEL : Objet MMODEL spécifiant la formulation 43 : C MRIGI1 : RIGIDITE masses hybrides élémentaires de sous-type DARCY. 44 : C MRIGI2 : Matrices masses hybrides elementaires de sous-type MASSE. 45 : C CHP2 : CHPOINT centre contenant la charge H (ou concentration) 46 : C CHP3 : CHPOINT contenant la trace de charge TH 47 : C (ou les traces de concentration) 48 : C CHPO4 : CHPO centre de composante FX,FY(,FZ), 49 : C densite de force moyenne par élément 50 : C CHPO5 : CHPOIN face contenant le flux de la vitesse convective 51 : C 52 : C 53 : C----------------------------------------------------------------------- 54 : C 55 : C Langage : ESOPE + FORTRAN77 56 : C 57 : C Auteurs : 09/93 F.DABBENE - Cas permanent 58 : C 09/94 X.NOUVELLON - Extension au cas transitoire 59 : C 02/96 L.V.BENET - Prise en compte de forces de volume 60 : C 06/98 F.AURIOL - formulation prenant en compte les valeurs 61 : C des charges ou des concentrations au centre 62 : C des elements. Prise en compte de la vitesse 63 : C convective 64 : C 65 : C----------------------------------------------------------------------- 66 : IMPLICIT INTEGER(I-N) 67 : IMPLICIT REAL*8 (A-H,O-Z) 68 : * 69 : -INC CCOPTIO 70 : -INC SMCHAML 71 : -INC SMCHPOI 72 : -INC SMELEME 73 : -INC SMMODEL 74 : -INC SMRIGID 75 : -INC SMTABLE 76 : * 77 : SEGMENT IPMAHY 78 : INTEGER MAHYBR(NSOUS) 79 : ENDSEGMENT 80 : * 81 : LOGICAL LOGRE,LOGIN 82 : CHARACTER*4 NOMTOT(1) 83 : CHARACTER*4 NOMTO3(3) 84 : CHARACTER*4 NOMTH 85 : CHARACTER*6 CHAR6 86 : CHARACTER*8 TAPIND,TYPOBJ,CHARIN,CHARRE,LETYPE 87 : * 88 : * Initialisations 89 : * 90 : IVALIN = 0 91 : XVALIN = 0.D0 92 : LOGIN = .TRUE. 93 : IOBIN = 0 94 : TAPIND = 'MOT ' 95 : * 96 : * Lecture du MMODEL et test de la formulation 97 : * 97 : CALL LIROBJ('MMODEL',IPMODE,1,IRET) 98 : IF (IERR.NE.0) RETURN 99 : MMODEL = IPMODE 100 : SEGACT MMODEL 101 : N1=KMODEL(/1) 102 : DO 7 I=1,N1 103 : IMODEL=KMODEL(I) 104 : SEGACT IMODEL 105 : IF(FORMOD(1).NE.'DARCY')THEN 106 : MOTERR(1:16) = 'DARCY ' 106 : CALL ERREUR(719) 107 : RETURN 108 : ENDIF 109 : SEGDES IMODEL 110 : 7 CONTINUE 111 : SEGDES MMODEL 112 : * 113 : * Lecture de la TABLE domaine 114 : * 115 : IPTABL = 0 115 : CALL LEKMOD(MMODEL,IPTABL,IRET) 116 : IF (IERR.NE.0) RETURN 117 : CHARIN = 'MAILLAGE' 118 : TYPOBJ = 'MAILLAGE' 118 : CALL LEKTAB(IPTABL,CHARIN,IOBRE) 119 : IF (IERR.NE.0) RETURN 120 : IPGEOM = IOBRE 120 : CALL LEKTAB(IPTABL,'ELTFA',IOBRE) 121 : IF (IERR.NE.0) RETURN 122 : IELTFA = IOBRE 122 : CALL LEKTAB(IPTABL,'CENTRE',IOBRE) 123 : IF (IERR.NE.0) RETURN 124 : ICENTR = IOBRE 124 : CALL LEKTAB(IPTABL,'FACE',IOBRE) 125 : IF (IERR.NE.0) RETURN 126 : IPFACE = IOBRE 127 : * 128 : * récup. MCHAML orientation normale 129 : * 129 : CALL LEKTAB(IPTABL,'XXNORMAE',IORIE) 130 : IF (IERR.NE.0) RETURN 131 : * 132 : * récup. CHPO orientation normale 133 : * 133 : CALL LEKTAB(IPTABL,'XXNORMAF',INORM) 134 : IF (IERR.NE.0) RETURN 135 : * 136 : * récup. CHPO surface des faces 137 : * 137 : CALL LEKTAB(IPTABL,'XXSURFAC',ISURF) 138 : IF (IERR.NE.0) RETURN 139 : * 140 : * Lecture de RIGIDITE obligatoire 141 : * 141 : CALL LIROBJ('RIGIDITE',IPRIGI,1,IRET) 142 : IF (IERR.NE.0) RETURN 143 : RI1 = IPRIGI 144 : * 145 : * Lecture du CHAMPOINT des charges (ou concentrations au centre) 146 : * 146 : CALL LIROBJ('CHPOINT',ITPP,1,IRET) 147 : * 148 : ICHP1 = ITPP 149 : * 150 : * Lecture du CHAMPOINT de concentration aux facees 151 : * 151 : CALL LIROBJ('CHPOINT',ITPN,1,IRET) 152 : IF (IERR.NE.0) RETURN 153 : * 154 : * Lecture eventuelle du 1er CHAMPOINT optionnel 155 : * 156 : ICHP2 = 0 156 : CALL LIROBJ('CHPOINT',ICHP2,0,IRET2) 157 : IF (IERR.NE.0) RETURN 158 : * 159 : * Lecture eventuelle du 2em CHAMPOINT optionnel 160 : * 161 : ICHP3 = 0 161 : CALL LIROBJ('CHPOINT',ICHP3,0,IRET3) 162 : IF (IERR.NE.0) RETURN 163 : * 164 : * Lecture de la RIGIDITE complementaire 165 : * 166 : IPRIGI = 0 166 : CALL LIROBJ('RIGIDITE',IPRIGI,0,IRETR) 167 : IF (IERR.NE.0) RETURN 168 : RI2 = IPRIGI 169 : * 170 : * Test de la formulation contenu dans MMODEL 171 : * 172 : SEGACT MMODEL 173 : NSOUS = KMODEL(/1) 174 : SEGINI IPMAHY 175 : IDARCY = 0 176 : IPT1= IPGEOM 177 : IPT2= IPGEOM 178 : DO 10 ISOUS=1,NSOUS 179 : IF(NSOUS.GT.1) THEN 180 : SEGACT IPT2 181 : IPT1= IPT2.LISOUS(ISOUS) 182 : SEGDES IPT2 183 : ENDIF 184 : IMODEL = KMODEL(ISOUS) 185 : SEGACT IMODEL 186 : LETYPE = FORMOD(1) 187 : IF (LETYPE.EQ.'DARCY') THEN 188 : IDARCY = IDARCY + 1 189 : MAHYBR(ISOUS) = IPT1 190 : ENDIF 191 : SEGDES IMODEL 192 : 10 CONTINUE 193 : SEGDES MMODEL 194 : IF (IDARCY.EQ.0) THEN 195 : MOTERR = LETYPE 195 : CALL ERREUR(193) 196 : GOTO 100 197 : ENDIF 198 : * 199 : * Récuperation des pointeurs ELTFA pour les zones ou DARCY est defini 200 : * 201 : MELEME = IELTFA 202 : SEGACT MELEME 203 : LZONES = LISOUS(/1) 204 : IF (LZONES.EQ.0) LZONES = 1 205 : IPT4 = IPGEOM 206 : SEGACT IPT4 207 : DO 30 ISOUS=1,NSOUS 208 : IMAMEL = MAHYBR(ISOUS) 209 : IF (IMAMEL.NE.0) THEN 210 : DO 20 ISZ=1,LZONES 211 : IF (LZONES.EQ.1) THEN 212 : IPT2 = IPT4 213 : IPT3 = MELEME 214 : ELSE 215 : IPT2 = IPT4.LISOUS(ISZ) 216 : IPT3 = LISOUS(ISZ) 217 : ENDIF 218 : IF (IPT2.EQ.IMAMEL) THEN 219 : MAHYBR(ISOUS) = IPT3 220 : GOTO 30 221 : ENDIF 222 : 20 CONTINUE 223 : IF (IMAMEL.EQ.MAHYBR(ISOUS)) THEN 224 : INTERR(1) = ISOUS 224 : CALL ERREUR(698) 225 : SEGDES IPT4 226 : SEGDES MELEME 227 : GOTO 100 228 : ENDIF 229 : ENDIF 230 : 30 CONTINUE 231 : SEGDES IPT4 232 : SEGDES MELEME 233 : * 234 : * Test et tri du sous-type des matrices de rigiditée récupérées 235 : * 236 : IPFORC=0 237 : IPDARC=0 238 : SEGACT RI1 239 : LETYPE = RI1.MTYMAT 240 : IF(LETYPE.EQ.'DARCY')THEN 241 : IPDARC=RI1 242 : ELSE IF(LETYPE.EQ.'MASSE')THEN 243 : IPFORC=RI1 244 : ENDIF 245 : IF(RI2.NE.0)THEN 246 : SEGACT RI2 247 : LETYPE = RI2.MTYMAT 248 : IF(LETYPE.EQ.'MASSE')THEN 249 : IPFORC=RI2 250 : ELSE IF(LETYPE.EQ.'DARCY')THEN 251 : IPDARC=RI2 252 : ENDIF 253 : IF(IPFORC.EQ.0)THEN 254 : MOTERR(1:8) = 'RIGIDITE' 255 : MOTERR(9:16) = 'MASSE' 255 : CALL ERREUR(79) 256 : SEGDES RI1 257 : SEGDES RI2 258 : GOTO 100 259 : ENDIF 260 : ENDIF 261 : IF(IPDARC.EQ.0)THEN 262 : MOTERR(1:8) = 'RIGIDITE' 263 : MOTERR(9:16) = 'DARCY' 263 : CALL ERREUR(79) 264 : SEGDES RI1 265 : IF(IPFORC.NE.0)SEGDES RI2 266 : GOTO 100 267 : ENDIF 268 : * 269 : * Controle des pointeurs de MELEME de la rigidité 270 : * 271 : DO 40 ISOUS=1,NSOUS 272 : IMAMEL = MAHYBR(ISOUS) 273 : IF (IMAMEL.NE.0) THEN 274 : RI1=IPDARC 275 : IPTEST = RI1.IRIGEL(1,ISOUS) 275 : CALL CMPMAI(IMAMEL,IPTEST,LOGRE) 276 : IF(.NOT. LOGRE)THEN 277 : MOTERR(1:8) = 'DARCY' 278 : MOTERR(9:16) = ' ELTFA ' 279 : INTERR(1) = ISOUS 279 : CALL ERREUR(698) 280 : SEGDES RI1 281 : GOTO 100 282 : ENDIF 283 : IF(IPFORC.NE.0)THEN 284 : RI2=IPFORC 285 : IPTEST = RI2.IRIGEL(1,ISOUS) 285 : CALL CMPMAI(IMAMEL,IPTEST,LOGRE) 286 : IF(.NOT. LOGRE)THEN 287 : MOTERR(1:8) = 'MASSE' 288 : MOTERR(9:16) = ' ELTFA ' 289 : INTERR(1) = ISOUS 289 : CALL ERREUR(698) 290 : SEGDES RI1 291 : SEGDES RI2 292 : GOTO 100 293 : ENDIF 294 : ENDIF 295 : ENDIF 296 : 40 CONTINUE 297 : SEGDES RI1 298 : IF(IPFORC.NE.0)SEGDES RI2 299 : * 300 : * Test du CHAMPOINT de valeurs aux faces 301 : * 302 : IFORC=0 303 : MCHPOI=ITPN 304 : SEGACT MCHPOI 305 : NBSOUS=IPCHP(/1) 306 : SEGDES MCHPOI 307 : ICHP0= ITPN 308 : IITH=0 309 : IVCONV=0 310 : IF(NBSOUS.NE.1)THEN 311 : C c'est le cas en presence de multiplicateurs de Lagrange 312 : C on va extraire la composante TH 313 : NOMTH='TH ' 313 : CALL EXCOPP(ITPN,NOMTH,ITPTH,NOMTH,0) 314 : IF(IERR.NE.0)GO TO 100 315 : ICHP0= ITPTH 316 : IITH=1 317 : ENDIF 318 : NOMTOT(1)=' ' 319 : INDIC=0 320 : NBCOMP=-1 320 : CALL QUEPOI(ICHP0,IPFACE,INDIC,NBCOMP,NOMTOT) 321 : IF(IERR.NE.0)GO TO 100 322 : IF(INDIC.LT.0)THEN 323 : IRETOU=0 324 : IPOINT=0 325 : CHARIN=' ' 325 : CALL SKNAME(ITPN,CHARIN,IRETOU,IPOINT) 326 : MOTERR(1:8)= CHARIN(1:8) 326 : IF(IRETOU.EQ.1)CALL ERREUR(-301) 326 : CALL ERREUR(609) 327 : GO TO 100 328 : ENDIF 329 : IF(IITH.EQ.0)THEN 330 : IF(NBSOUS.EQ.1)THEN 331 : MCHPOI=ICHP0 332 : SEGACT MCHPOI 333 : MSOUPO=IPCHP(1) 334 : SEGDES MCHPOI 335 : SEGACT MSOUPO 336 : NC=NOCOMP(/2) 337 : IF(NC.EQ.1)THEN 338 : IF(NOCOMP(1).EQ.'TH ')IITH=1 339 : ENDIF 340 : SEGDES MSOUPO 341 : ENDIF 342 : ENDIF 343 : * 344 : * Tri et test des CHAMPOINTs ARGUMENTS optionnels 345 : * 346 : IF(IPFORC.NE.0) THEN 347 : IF(IITH.NE.1)THEN 348 : IRETOU=0 349 : IPOINT=0 350 : CHARIN=' ' 350 : CALL SKNAME(ITPN,CHARIN,IRETOU,IPOINT) 351 : MOTERR(1:4)= 'TH ' 352 : MOTERR(5:12)= CHARIN(1:8) 352 : CALL ERREUR(77) 353 : ENDIF 354 : IF(IRET2.EQ.1)THEN 355 : IFORC=ICHP2 356 : NOMTOT(1)=' ' 357 : INDIC=1 358 : NBCOMP=-1 358 : CALL QUEPOI(ICHP1,ICENTR,INDIC,NBCOMP,NOMTOT) 359 : IF(IERR.NE.0)GO TO 100 360 : IF(INDIC.LT.0)THEN 360 : CALL ERREUR(348) 361 : GO TO 100 362 : ENDIF 363 : ELSE 363 : CALL ERREUR(641) 364 : GO TO 100 365 : ENDIF 366 : NBCOMP = 0 367 : INDIC = 1 368 : NOMTO3(1) = ' ' 369 : NOMTO3(2) = ' ' 370 : NOMTO3(3) = ' ' 370 : CALL QUEPOI(IFORC,ICENTR,INDIC,NBCOMP,NOMTO3) 371 : IF(NOMTO3(1).NE.'FX')THEN 372 : MOTERR(1:4)= NOMTO3(1) 372 : CALL ERREUR(197) 373 : ENDIF 374 : IF(IRET3.NE.0)THEN 375 : NBCOMP = -1 376 : INDIC = 0 377 : NOMTOT(1) = ' ' 377 : CALL QUEPOI(ICHP3,IPFACE,INDIC,NBCOMP,NOMTOT) 378 : IF(IERR.NE.0)GO TO 100 379 : IF(INDIC.LT.0)THEN 379 : CALL ERREUR(348) 380 : GO TO 100 381 : ENDIF 382 : IVCONV=ICHP3 383 : ENDIF 384 : C Cas ou il n y a pas de rigidite RI2 385 : ELSE 386 : NOMTOT(1)=' ' 387 : INDIC=1 388 : NBCOMP=-1 388 : CALL QUEPOI(ICHP1,ICENTR,INDIC,NBCOMP,NOMTOT) 389 : IF(IERR.NE.0)GO TO 100 390 : IF(INDIC.LT.0)THEN 390 : CALL ERREUR(348) 391 : GO TO 100 392 : ENDIF 393 : MCHPOI=ICHP0 394 : SEGACT MCHPOI 395 : MSOUPO=IPCHP(1) 396 : SEGDES MCHPOI 397 : SEGACT MSOUPO 398 : NC=NOCOMP(/2) 399 : MCHPO1= ICHP1 400 : SEGACT MCHPO1 401 : MSOUP1=MCHPO1.IPCHP(1) 402 : SEGDES MCHPO1 403 : SEGACT MSOUP1 404 : NC1=MSOUP1.NOCOMP(/2) 405 : IF(NC1.NE.NC)THEN 406 : MOTERR(1:40)= 407 : & 'pas le meme nombre de composantes ' 407 : CALL ERREUR(-301) 407 : CALL ERREUR(21) 408 : GO TO 100 409 : ENDIF 410 : IF(IITH.EQ.1) THEN 411 : IF((MSOUP1.NOCOMP(1)).NE.'H ')THEN 412 : MOTERR(1:40)= 413 : & 'le nom de la composante doit etre H ' 413 : CALL ERREUR(-301) 413 : CALL ERREUR(21) 414 : GO TO 100 415 : ENDIF 416 : IF(IRET2.NE.0)THEN 417 : NBCOMP = -1 418 : INDIC = 0 419 : NOMTOT(1) = ' ' 419 : CALL QUEPOI(ICHP2,IPFACE,INDIC,NBCOMP,NOMTOT) 420 : IF(IERR.NE.0)GO TO 100 421 : IF(INDIC.LT.0)THEN 421 : CALL ERREUR(348) 422 : GO TO 100 423 : ENDIF 424 : IVCONV=ICHP2 425 : ENDIF 426 : ELSE 427 : NNC1=0 428 : DO 50 I=1,NC 429 : DO 35 J=1,NC1 430 : IF(NOCOMP(I).EQ.MSOUP1.NOCOMP(J))THEN 431 : NNC1=NNC1+1 432 : GO TO 32 433 : ENDIF 434 : 35 CONTINUE 435 : MOTERR(1:40)= 436 : & 'les noms de composantes sont differents ' 436 : CALL ERREUR(-301) 436 : CALL ERREUR(21) 437 : GO TO 100 438 : 32 CONTINUE 439 : 50 CONTINUE 440 : IF(NNC1.NE.NC1)THEN 441 : MOTERR(1:40)= 442 : & 'un des CHPO a plus de comp. que l autre ' 442 : CALL ERREUR(21) 443 : GO TO 100 444 : ENDIF 445 : IVCONV=ICHP2 446 : ENDIF 447 : SEGDES MSOUPO,MSOUP1 448 : ENDIF 449 : * 450 : * Construction du CHAMPOINT des debits aux faces 451 : * 452 : SEGDES IPMAHY 452 : CALL HDEBI2(IPMAHY,IPFACE,IPDARC,ICHP0,IORIE, 453 : S ICHP1,IVCONV,IITH,INORM,IORIE,ISURF,IPFORC,IFORC,IRET) 454 : 455 : 455 : IF (IRET.NE.0) CALL ECROBJ('CHPOINT',IRET) 456 : * 457 : * Ménage 458 : * 459 : 100 CONTINUE 460 : SEGSUP IPMAHY 461 : RETURN 462 : END 463 : 464 : 465 : 466 : 467 : 468 : 469 : 470 : 471 : 472 : 473 : 474 : 475 : 476 : 477 : 478 : 479 : 480 : 481 : 482 :
© Cast3M 2003 - All rights reserved.
Disclaimer