1 : C SMTP1 SOURCE CHAT 09/10/09 21:24:23 6519 2 : SUBROUTINE SMTP1(IPMAHY,IPFACE,IPRIGI,THETA,THETAC,DELTAT,IPCK, 3 : S ITP,IH,ISOUR,IFLUX,IORIE,MCHPOI,IRI1) 4 : C----------------------------------------------------------------------- 5 : C 1) Calcul du second membre du systeme en trace de charge dans le cas 6 : C de la résolution des équations de Darcy par EFMH avec le modèle DARCY. 7 : C 2) Calcul du second membre du systeme en trace de charge et/ou d'une 8 : C contribution matricielle provenant de la convection dans le cas de la 9 : C résolution d'une équation de diffusion-convection par EFMH avec le 10 : C modèle DARCY. 11 : C----------------------------------------------------------------------- 12 : C 13 : C--------------------------- 14 : C Parametres Entree/Sortie : 15 : C--------------------------- 16 : C 17 : C E/ IPMAHY : Segment contenant le pointeur vers le meleme des 18 : C connectivites elements/faces pour les zones du MMODEL 19 : C ou on a defini DARCY. 20 : C E/ IPFACE : MELEME de type POI1 des faces 21 : C E/ IPRIGI : RIGIDITE de sous type 'DARCY' 22 : C E/ THETA : Parametre de discretisation temporelle (Theta-methode) 23 : C pour le terme de diffusion 24 : C E/ THETAC : Parametre de discretisation temporelle (Theta-methode) 25 : C pour le terme de convection 26 : C E/ DELTAT : Pas de temps 27 : C E/ IPCK : MCHAML donnant pour chaque element Ck|K| 28 : C E/ ITP : CHPO faces des traces de charge TH au temps precedent 29 : C E/ IH : CHPO centre des charges H au temps precedent 30 : C E/ ISOUR : CHPO centre des sources de composante SOUR 31 : C E/ IFLUX : CHPO face des flux de composante FLUX 32 : C E/ IORIE : MCHAML orientation des normales 33 : C /S MCHPOI : CHPOINT face de composante FLUX 34 : C /S IRI1 : RIGIDITE de sous type CONVEFMH 35 : C 36 : C---------------------- 37 : C Variables en COMMON : 38 : C---------------------- 39 : C 40 : C E/ IFOMOD : cf CCOPTIO 41 : C E/ NIFOUR : cf CCOPTIO 42 : C 43 : C--------------------- 44 : C Variables internes : 45 : C--------------------- 46 : C 47 : C LHS1 = 1 si on doit créer une matrice (Left Hand Side); 0 sinon. 48 : C RHS1 = 1 si on doit créer un second membre (Right Hand Side); 0 sinon. 49 : C 50 : C----------------------------------------------------------------------- 51 : C 52 : C Langage : ESOPE + FORTRAN77 53 : C 54 : C Auteurs : 09/93 F.DABBENE - Cas permanent 55 : C 09/94 X.NOUVELLON - Extension au cas transitoire 56 : C 12/95 F.DABBENE - Extension à la convection 57 : C 58 : C----------------------------------------------------------------------- 59 : IMPLICIT INTEGER(I-N) 60 : IMPLICIT REAL*8 (A-H,O-Z) 61 : * 62 : -INC CCOPTIO 63 : -INC CCHAMP 64 : -INC SMCOORD 65 : -INC SMCHPOI 66 : -INC SMCHAML 67 : -INC SMELEME 68 : -INC SMRIGID 69 : * 70 : SEGMENT IPMAHY 71 : INTEGER MAHYBR(NSOUS) 72 : ENDSEGMENT 73 : SEGMENT ICCPR 74 : INTEGER ICPR(NNGOT) 75 : ENDSEGMENT 76 : SEGMENT MTRAV 77 : REAL*8 RLIGNE(NBDDL),RCOLON(NBDDL),FLUORI(NBDDL) 78 : ENDSEGMENT 79 : * 80 : *- Initialisations 81 : * 82 : MCHPOI = 0 83 : IRI1 = 0 84 : IF (IPCK.EQ.0) THEN 85 : IF (ISOUR.EQ.0) THEN 86 : RHS1 = 0 87 : ELSE 88 : RHS1 = 1 89 : ENDIF 90 : IF (IFLUX.EQ.0) THEN 91 : LHS1 = 0 92 : ELSE 93 : LHS1 = 1 94 : ENDIF 95 : ELSE 96 : RHS1 = 1 97 : LHS1 = 0 98 : IF (THETAC.GT.1.D-4 .AND. IFLUX.NE.0) LHS1=1 99 : ENDIF 100 : * 101 : *- Creation du tableau ICPR pour le MELEME POI1 IPFACE 102 : * 103 : IK = 0 104 : NNGOT = XCOOR(/1)/(IDIM+1) 105 : SEGINI ICCPR 106 : MELEME = IPFACE 107 : SEGACT MELEME 108 : N2 = NUM(/2) 109 : DO 10 I2=1,N2 110 : K = NUM(1,I2) 111 : IF (ICPR(K).EQ.0) THEN 112 : IK = IK + 1 113 : ICPR(K) = IK 114 : ENDIF 115 : 10 CONTINUE 116 : SEGDES MELEME 117 : * 118 : *- Creation du CHAMPOINT face - On laisse actif le MPOVAL du CHPOINT 119 : * 120 : IF (RHS1.NE.0) THEN 121 : NSOUPO = 1 122 : NAT = 1 123 : SEGINI MCHPOI 124 : IRET = MCHPOI 125 : MTYPOI = 'FACE ' 126 : IFOPOI = IFOMOD 127 : JATTRI(1) = 2 128 : NC = 1 129 : SEGINI MSOUPO 130 : IPCHP(1) = MSOUPO 131 : SEGDES MCHPOI 132 : NOHARM(1) = NIFOUR 133 : IGEOC = IPFACE 134 : NOCOMP(1) = 'FLUX' 135 : N = N2 136 : SEGINI MPOVAL 137 : IPOVAL = MPOVAL 138 : SEGDES MSOUPO 139 : ENDIF 140 : * 141 : *- Recuperation du nombre de zone NBMAIL du MMODEL a partir IPMAHY 142 : * 143 : SEGACT IPMAHY 144 : NBMAIL = MAHYBR(/1) 145 : * 146 : *- Activation du segment RIGIDITE 147 : * 148 : MRIGID = IPRIGI 149 : SEGACT MRIGID 150 : * 151 : *- Activation du segment MPOVAL du CHAMPOINT centre de composante SOUR 152 : * 153 : IF (ISOUR.NE.0) THEN 154 : MCHPO1 = ISOUR 155 : SEGACT MCHPO1 156 : MSOUP1 = MCHPO1.IPCHP(1) 157 : SEGDES MCHPO1 158 : SEGACT MSOUP1 159 : MPOVA1 = MSOUP1.IPOVAL 160 : SEGDES MSOUP1 161 : SEGACT MPOVA1 162 : ENDIF 163 : * 164 : *- Données transitoire 165 : * 166 : IF (IPCK.NE.0) THEN 167 : * 168 : * Activation du segment MPOVAL du CHAMPOINT face de traces de charge 169 : * 170 : MCHPO2 = ITP 171 : SEGACT MCHPO2 172 : MSOUP2 = MCHPO2.IPCHP(1) 173 : SEGDES MCHPO2 174 : SEGACT MSOUP2 175 : MPOVA2 = MSOUP2.IPOVAL 176 : SEGDES MSOUP2 177 : SEGACT MPOVA2 178 : * 179 : * Activation du segment MPOVAL du CHAMPOINT centre de charges 180 : * 181 : MCHPO3 = IH 182 : SEGACT MCHPO3 183 : MSOUP3 = MCHPO3.IPCHP(1) 184 : SEGDES MCHPO3 185 : SEGACT MSOUP3 186 : MPOVA3 = MSOUP3.IPOVAL 187 : SEGDES MSOUP3 188 : SEGACT MPOVA3 189 : * 190 : * Activation du MCHAML contenant Ck|K| 191 : * 192 : MCHELM = IPCK 193 : SEGACT MCHELM 194 : ENDIF 195 : * 196 : *- Activation des données pour la convection 197 : * 198 : IF (IFLUX.NE.0) THEN 199 : * 200 : * Activation du MCHAML d'orientation des normales 201 : * 202 : MCHEL1 = IORIE 203 : SEGACT MCHEL1 204 : * 205 : * Activation du segment MPOVAL du CHAMPOINT flux aux faces 206 : * 207 : MCHPO4 = IFLUX 208 : SEGACT MCHPO4 209 : MSOUP4 = MCHPO4.IPCHP(1) 210 : SEGDES MCHPO4 211 : SEGACT MSOUP4 212 : MPOVA4 = MSOUP4.IPOVAL 213 : SEGDES MSOUP4 214 : SEGACT MPOVA4 215 : * 216 : * Initialisation du chapeau de rigidité 217 : * 218 : IF (LHS1.NE.0) THEN 219 : NRIGE = 7 220 : NRIGEL = NBMAIL 221 : SEGINI RI1 222 : IRI1 = RI1 223 : RI1.MTYMAT = 'CONVEFMH' 224 : RI1.IFORIG = IFORIG 225 : ENDIF 226 : ENDIF 227 : * 228 : *-------------------------------------------------- 229 : *= BOUCLE SUR LES MAILLAGES ELEMENTAIRES,ZONE IMAIL 230 : *-------------------------------------------------- 231 : * 232 : ITELEM = 0 233 : DO 110 IMAIL=1,NBMAIL 234 : C 235 : C- Activation de l'objet maillage ELTFA pour la zone IMAIL 236 : C 237 : MELEME = MAHYBR(IMAIL) 238 : IF (MELEME.EQ.0) GOTO 110 239 : SEGACT MELEME 240 : * 241 : *- Recuperation des caracteristiques de RIGIDITE de sous type DARCY 242 : *- pour la zone IMAIL en cours de traitement 243 : * 244 : xMATRI = IRIGEL(4,IMAIL) 245 : SEGACT xMATRI 246 : NBELEM = re(/3) 247 : NELRIG = NBELEM 248 : * XMATRI = IMATTT(1) 249 : * SEGACT XMATRI 250 : NBDDL = RE(/1) 251 : NLIGRP = NBDDL 252 : NLIGRD = NBDDL 253 : * SEGDES XMATRI 254 : SEGINI MTRAV 255 : * 256 : *- Activation du segment contenant les valeurs 257 : *- du MCHAML Ck|K| pour la zone IMAIL 258 : * 259 : IF (IPCK.NE.0) THEN 260 : MCHAML = ICHAML(IMAIL) 261 : SEGACT MCHAML 262 : MELVAL = IELVAL(1) 263 : SEGACT MELVAL 264 : ENDIF 265 : * 266 : *- Convection 267 : * 268 : IF (IFLUX.NE.0) THEN 269 : * 270 : * Activation du segment contenant les valeurs du MCHAML d'orientation 271 : * des normales par face pour la zone IMAIL 272 : * 273 : MCHAM1 = MCHEL1.ICHAML(IMAIL) 274 : SEGACT MCHAM1 275 : MELVA1 = MCHAM1.IELVAL(1) 276 : SEGACT MELVA1 277 : * 278 : * Initialisation du chapeau de rigidité pour la zone IMAIL 279 : * 280 : IF (LHS1.NE.0) THEN 281 : RI1.COERIG(IMAIL) = 1.D0 282 : RI1.IRIGEL(1,IMAIL) = IRIGEL(1,IMAIL) 283 : RI1.IRIGEL(3,IMAIL) = IRIGEL(3,IMAIL) 284 : SEGINI xMATR1 285 : RI1.IRIGEL(4,IMAIL) = xMATR1 286 : RI1.IRIGEL(7,IMAIL) = 2 287 : ENDIF 288 : ENDIF 289 : * 290 : *------------------------------------------------------- 291 : * BOUCLE SUR LES ELEMENTS DU MAILLAGE ELEMENTAIRE IMAIL 292 : *------------------------------------------------------- 293 : * 294 : DO 100 IEL=1,NBELEM 295 : ITELEM = ITELEM + 1 296 : * 297 : *- Orientation du flux et calcul de la source = -div(u*Th) 298 : * 299 : DIVUTH = 0.D0 300 : IF (IFLUX.NE.0) THEN 301 : DO 20 I=1,NBDDL 302 : IPOPTS = ICPR(NUM(I,IEL)) 303 : FLUORI(I)= MPOVA4.VPOCHA(IPOPTS,1)*MELVA1.VELCHE(I,IEL) 304 : DIVUTH = DIVUTH - (FLUORI(I)*MPOVA2.VPOCHA(IPOPTS,1)) 305 : 20 CONTINUE 306 : ENDIF 307 : * 308 : *- Calcul de la somme des coefs pour une ligne ; une colonne 309 : *- -1 t -1 310 : *- LIGNE = RE * DIV ; COLON = DIV * RE 311 : *- -1 t 312 : *- Calcul de CONSD = DIV * RE * DIV 313 : * 314 : CONSD = 0.D0 315 : * XMATRI = IMATTT(IEL) 316 : * SEGACT XMATRI 317 : DO 40 I=1,NBDDL 318 : RCOLON(I) = 0.D0 319 : RLIGNE(I) = 0.D0 320 : DO 30 J=1,NBDDL 321 : RCOLON(I) = RCOLON(I) + RE(J,I,iel) 322 : RLIGNE(I) = RLIGNE(I) + RE(I,J,iel) 323 : 30 CONTINUE 324 : CONSD = CONSD + RLIGNE(I) 325 : 40 CONTINUE 326 : * SEGDES XMATRI 327 : * 328 : *- Recuperation de la valeur SOUR au centre de l'element 329 : *- ET ajout du terme de convection éventuelle (1-thetac)*(-div(u*Th)) 330 : * 331 : IF (IPCK.EQ.0) THEN 332 : VSOU = MPOVA1.VPOCHA(ITELEM,1) 333 : ELSE 334 : IF (ISOUR.NE.0) THEN 335 : VSOU = MPOVA1.VPOCHA(ITELEM,1)+((1.D0-THETAC)*DIVUTH) 336 : ELSE 337 : VSOU = (1.D0-THETAC) * DIVUTH 338 : ENDIF 339 : ENDIF 340 : * 341 : * 342 : IF (IPCK.EQ.0) THEN 343 : IF (RHS1.NE.0) THEN 344 : DO 45 IDDL=1,NBDDL 345 : VALFLU = RLIGNE(IDDL) * VSOU / CONSD 346 : IPOPTS = ICPR(NUM(IDDL,IEL)) 347 : VPOCHA(IPOPTS,1) = VPOCHA(IPOPTS,1) - VALFLU 348 : 45 CONTINUE 349 : ENDIF 350 : IF (LHS1.NE.0) THEN 351 : * SEGINI XMATR1 352 : * IMATR1.IMATTT(IEL) = XMATR1 353 : DO 55 J=1,NBDDL 354 : REMU1 = -1.D0 / CONSD 355 : DO 50 I=1,NBDDL 356 : XMATR1.RE(I,J,iel) = REMU1*RLIGNE(I)*FLUORI(J) 357 : 50 CONTINUE 358 : 55 CONTINUE 359 : * SEGDES XMATR1 360 : ENDIF 361 : ELSE 362 : BETA = CONSD*DELTAT/(VELCHE(1,IEL)+THETA*CONSD*DELTAT) 363 : VH = MPOVA3.VPOCHA(ITELEM,1) 364 : DO 70 IDDL=1,NBDDL 365 : SMFLU = RLIGNE(IDDL) * VSOU * BETA/CONSD 366 : SMP = RLIGNE(IDDL) * VH * (1.D0-BETA) 367 : STP = 0.D0 368 : DO 60 J=1,NBDDL 369 : VTH = MPOVA2.VPOCHA(ICPR(NUM(J,IEL)),1) 370 : STP = STP + RCOLON(J) * VTH 371 : 60 CONTINUE 372 : SMTP = RLIGNE(IDDL) * STP * (1.D0-THETA)*BETA/CONSD 373 : SMDDL = SMFLU + SMP + SMTP 374 : IPOPTS = ICPR(NUM(IDDL,IEL)) 375 : VPOCHA(IPOPTS,1) = VPOCHA(IPOPTS,1) - SMDDL 376 : 70 CONTINUE 377 : IF (LHS1.NE.0) THEN 378 : * SEGINI XMATR1 379 : * IMATR1.IMATTT(IEL) = XMATR1 380 : REMU1 = (-1.D0) * THETAC * BETA / CONSD 381 : DO 90 J=1,NBDDL 382 : DO 80 I=1,NBDDL 383 : XMATR1.RE(I,J,iel)=REMU1*RLIGNE(I) * FLUORI(J) 384 : 80 CONTINUE 385 : 90 CONTINUE 386 : * SEGDES XMATR1 387 : ENDIF 388 : ENDIF 389 : 100 CONTINUE 390 : SEGDES xMATRI 391 : SEGDES MELEME 392 : IF (IPCK.NE.0) THEN 393 : SEGDES MCHAML 394 : SEGDES MELVAL 395 : ENDIF 396 : IF (IFLUX.NE.0) THEN 397 : SEGDES MCHAM1 398 : SEGDES MELVA1 399 : IF (LHS1.NE.0) SEGDES xMATR1 400 : ENDIF 401 : SEGSUP MTRAV 402 : 110 CONTINUE 403 : * 404 : *- Desactivation des segments 405 : * 406 : IF (RHS1.NE.0) THEN 407 : SEGDES MPOVAL 408 : ENDIF 409 : IF (IPCK.NE.0) THEN 410 : SEGDES MCHELM 411 : SEGDES MPOVA2 412 : SEGDES MPOVA3 413 : ENDIF 414 : IF (IFLUX.NE.0) THEN 415 : SEGDES MCHEL1 416 : SEGDES MPOVA4 417 : IF (LHS1.NE.0) SEGDES RI1 418 : ENDIF 419 : IF (ISOUR.NE.0) THEN 420 : SEGDES MPOVA1 421 : ENDIF 422 : SEGDES IPMAHY 423 : SEGDES MRIGID 424 : SEGSUP ICCPR 425 : * 426 : RETURN 427 : END 428 : 429 : 430 :
© Cast3M 2003 - All rights reserved.
Disclaimer