1 : C HYBP1 SOURCE CHAT 09/10/09 21:18:48 6519 2 : SUBROUTINE HYBP1(IPMAHY,ICENTR,IPFACE,IPDARC,ITPN,INORM,IORIE, 3 : S ISURF,THETA,DELTAT,IPCK,ITP,IH,ISOUR,IPFORC,IFORC,IRET) 4 : C----------------------------------------------------------------------- 5 : C Calcul de la charge au centre de chaque maille 6 : C----------------------------------------------------------------------- 7 : C 8 : C--------------------------- 9 : C Parametres Entree/Sortie : 10 : C--------------------------- 11 : C 12 : C E/ IPMAHY : Segment contenant le pointeur vers le meleme des 13 : C connectivites elements/faces pour les zones du MMODEL 14 : C ou on a defini DARCY. 15 : C E/ ICENTR : Pointeur vers l'objet maillage CENTRE 16 : C E/ IPFACE : Pointeur vers l'objet maillage FACE -1 17 : C E/ IPDARC : Matrice de sous type DARCY (contient RE ). 18 : C E/ ITPN : CHPO face des Traces de charge TH au temps n 19 : C E/ THETA : Parametre de discretisation temporelle (Theta-methode) 20 : C E/ DELTAT : Pas de temps 21 : C E/ IPCK : MCHAML donnant pour chaque element Ck|K| 22 : C E/ ITP : CHPO face des traces de charge TH au temps n-1 23 : C E/ IH : CHPO centre des charges H au temps n-1 24 : C E/ ISOUR : CHPO centre SOUR contenant un eventuel terme source 25 : C E/ IPFORC : Matrice de sous type MASSE (contient RE ). 26 : C E/ IFORC : CHPO centre force contenant la force volumique 27 : C /S IRET : CHPO centre H resultat contenant la charge au temps n 28 : C 29 : C---------------------- 30 : C Variables en COMMON : 31 : C---------------------- 32 : C 33 : C E/ IDIM : cf CCOPTIO 34 : C E/ IFOMOD : cf CCOPTIO 35 : C E/ NIFOUR : cf CCOPTIO 36 : C 37 : C---------------------- 38 : C Tableaux de travail : 39 : C---------------------- 40 : C 41 : C ICPR(I)=J : Le noeud I a le numero local J 42 : C Correspondance numerotation globale/locale 43 : C NNGOT : Nombre de noeuds total du domaine 44 : C 45 : C----------------------------------------------------------------------- 46 : C 47 : C Langage : ESOPE + FORTRAN77 48 : C 49 : C Auteurs : 09/93 F.DABBENE - Cas permanent 50 : C 09/94 X.NOUVELLON - Extension au cas transitoire 51 : C 02/96 L.V.BENET - Prise en compte de forces de volume 52 : C 53 : C----------------------------------------------------------------------- 54 : IMPLICIT INTEGER(I-N) 55 : IMPLICIT REAL*8 (A-H,O-Z) 56 : * 57 : -INC CCOPTIO 58 : -INC SMCHAML 59 : -INC SMCHPOI 60 : POINTEUR MCHPO5.MCHPOI,MCHPO6.MCHPOI,MSOUP6.MSOUPO 61 : -INC SMCOORD 62 : -INC SMELEME 63 : -INC SMRIGID 64 : * 65 : CHARACTER*4 MOREFD, MOREFP 66 : CHARACTER*1 BLAN1 67 : SEGMENT ICCPR 68 : INTEGER ICPR(NNGOT) 69 : ENDSEGMENT 70 : SEGMENT IORGA 71 : REAL*8 VAA(ITES) 72 : ENDSEGMENT 73 : SEGMENT IPMAHY 74 : INTEGER MAHYBR(NSOUS) 75 : ENDSEGMENT 76 : SEGMENT MTRAV 77 : REAL*8 VECT(NBDDL),FLFOR(NBDDL),RFOR(NBDDL) 78 : ENDSEGMENT 79 : * 80 : BLAN1=' ' 81 : * 82 : * Initialisations 83 : * 84 : IRET = 0 85 : IPT1 = IPFACE 86 : MRIGID = IPDARC 87 : RI1 = IPFORC 88 : MCHPOI = ITPN 89 : MCHELM = IPCK 90 : MCHEL1 = IORIE 91 : MCHPO1 = ISOUR 92 : MCHPO2 = ITP 93 : MCHPO3 = IH 94 : MCHPO4 = IFORC 95 : MCHPO5 = INORM 96 : MCHPO6 = ISURF 97 : VSOUR = 0.D0 98 : * 99 : * Initialisations se servant de IPMAHY et MRIGID 100 : * 101 : SEGACT IPMAHY 102 : NBMAIL = MAHYBR(/1) 103 : SEGACT MRIGID 104 : DO 5 IMAIL=1,NBMAIL 105 : IF (MAHYBR(IMAIL).NE.0) THEN 106 : DESCR = IRIGEL(3,IMAIL) 107 : SEGACT DESCR 108 : MOREFP = LISINC(1) 109 : SEGDES DESCR 110 : GOTO 10 111 : ENDIF 112 : 5 CONTINUE 113 : 10 CONTINUE 114 : NNGOT = XCOOR(/1)/(IDIM+1) 115 : SEGINI ICCPR 116 : * 117 : *= Creation du tableau ICPR pour le maillage IPT1 118 : * 119 : SEGACT IPT1 120 : N2 = IPT1.NUM(/2) 121 : IK = 0 122 : DO 15 I2=1,N2 123 : K = IPT1.NUM(1,I2) 124 : IF (ICPR(K).EQ.0) THEN 125 : IK = IK + 1 126 : ICPR(K) = IK 127 : ENDIF 128 : 15 CONTINUE 129 : SEGDES IPT1 130 : ITES = IK 131 : * 132 : *= Extraction des TH et réduction sur le support IPT1 du CHPO ITPN 133 : *= Les valeurs du CHPOINT en TH sont dans VAA avec l'ordre IPT1 134 : * 135 : SEGINI IORGA 136 : SEGACT MCHPOI 137 : NSOUPO = IPCHP(/1) 138 : DO 30 ISOU=1,NSOUPO 139 : MSOUPO = IPCHP(ISOU) 140 : SEGACT MSOUPO 141 : NC = NOCOMP(/2) 142 : NUMHAR = NOHARM(1) 143 : ITTP = 0 144 : DO 20 I=1,NC 145 : IF (NOCOMP(I).EQ.MOREFP) THEN 146 : ITTP = I 147 : ENDIF 148 : 20 CONTINUE 149 : IF (ITTP.NE.0) THEN 150 : MELEME = IGEOC 151 : SEGACT MELEME 152 : N2 = NUM(/2) 153 : MPOVAL = IPOVAL 154 : SEGACT MPOVAL 155 : DO 25 I2=1,N2 156 : K = NUM(1,I2) 157 : IF (ICPR(K).EQ.0) THEN 158 : INTERR(1) = K 159 : MOTERR(1:8) = 'FACE ' 159 : CALL ERREUR(64) 160 : SEGDES MPOVAL, MELEME, MSOUPO 161 : SEGDES MCHPOI, MRIGID, IPMAHY 162 : SEGSUP ICCPR, IORGA 163 : RETURN 164 : ELSE 165 : VAA(ICPR(K)) = VPOCHA(I2,ITTP) 166 : ENDIF 167 : 25 CONTINUE 168 : SEGDES MPOVAL, MELEME 169 : ENDIF 170 : SEGDES MSOUPO 171 : 30 CONTINUE 172 : SEGDES MCHPOI 173 : * 174 : *= Création des chapeaux du CHPO centre P résultat 175 : * 176 : NSOUPO = 1 177 : NAT = 1 178 : SEGINI MCHPOI 179 : IRET = MCHPOI 180 : DO 31 ITIT=1,72 181 : MOCHDE(ITIT:ITIT)=BLAN1 182 : 31 CONTINUE 183 : MTYPOI = 'CENTRE ' 184 : IFOPOI = IFOMOD 185 : JATTRI(1) = 1 186 : NC = 1 187 : SEGINI MSOUPO 188 : IPCHP(1) = MSOUPO 189 : SEGDES MCHPOI 190 : NOHARM(1) = NUMHAR 191 : IGEOC = ICENTR 192 : NOCOMP(1) = 'H ' 193 : IPT2 = ICENTR 194 : SEGACT IPT2 195 : N = IPT2.NUM(/2) 196 : SEGDES IPT2 197 : SEGINI MPOVAL 198 : IPOVAL = MPOVAL 199 : SEGDES MSOUPO 200 : * 201 : * Activation du MPOVAL du CHPO centre SOUR 202 : * 203 : IF (ISOUR.NE.0) THEN 204 : SEGACT MCHPO1 205 : MSOUP1 = MCHPO1.IPCHP(1) 206 : SEGDES MCHPO1 207 : SEGACT MSOUP1 208 : MPOVA1 = MSOUP1.IPOVAL 209 : SEGDES MSOUP1 210 : SEGACT MPOVA1 211 : ENDIF 212 : * 213 : * Activation du MPOVAL du CHPO face trace de charge au temps (n-1) 214 : * 215 : IF (ITP.NE.0) THEN 216 : SEGACT MCHPO2 217 : MSOUP2 = MCHPO2.IPCHP(1) 218 : SEGDES MCHPO2 219 : SEGACT MSOUP2 220 : MPOVA2 = MSOUP2.IPOVAL 221 : SEGDES MSOUP2 222 : SEGACT MPOVA2 223 : ENDIF 224 : * 225 : * Activation du MPOVAL du CHPO centre charge au temps (n-1) 226 : * 227 : IF (IH.NE.0) THEN 228 : SEGACT MCHPO3 229 : MSOUP3 = MCHPO3.IPCHP(1) 230 : SEGDES MCHPO3 231 : SEGACT MSOUP3 232 : MPOVA3 = MSOUP3.IPOVAL 233 : SEGDES MSOUP3 234 : SEGACT MPOVA3 235 : ENDIF 236 : * 237 : * activation des objets liés à la présence d'une force volumique 238 : * 239 : IF (IFORC.NE.0) THEN 240 : * 241 : * Activation du MPOVAL du CHPO force appuyé au centre des éléments volumiques 242 : * 243 : SEGACT MCHPO4 244 : MSOUP4 = MCHPO4.IPCHP(1) 245 : SEGDES MCHPO4 246 : SEGACT MSOUP4 247 : MPOVA4 = MSOUP4.IPOVAL 248 : SEGDES MSOUP4 249 : SEGACT MPOVA4 250 : * 251 : * Activation du MPOVAL du CHPO des vecteurs normales appuyé au centre des faces 252 : * 253 : SEGACT MCHPO5 254 : MSOUP5 = MCHPO5.IPCHP(1) 255 : SEGDES MCHPO5 256 : SEGACT MSOUP5 257 : MPOVA5 = MSOUP5.IPOVAL 258 : SEGDES MSOUP5 259 : SEGACT MPOVA5 260 : * 261 : * Activation du MPOVAL du CHPO des surfaces appuyé au centre des faces 262 : * 263 : SEGACT MCHPO6 264 : MSOUP6 = MCHPO6.IPCHP(1) 265 : SEGDES MCHPO6 266 : SEGACT MSOUP6 267 : MPOVA6 = MSOUP6.IPOVAL 268 : SEGDES MSOUP6 269 : SEGACT MPOVA6 270 : * 271 : * Activation du MCHEL des orientations des faces 272 : * 273 : SEGACT MCHEL1 274 : * 275 : * Activation du MRIGI de la matrice masse hybride 276 : * 277 : SEGACT RI1 278 : ENDIF 279 : * 280 : * Activation de Ck|K| 281 : * 282 : IF (IPCK.NE.0) SEGACT MCHELM 283 : * 284 : *--------------------------------------- 285 : *= Boucle sur les maillages elementaires 286 : *--------------------------------------- 287 : * 288 : ITELEM = 0 289 : DO 110 IMAIL=1,NBMAIL 290 : IF (MAHYBR(IMAIL).EQ.0) GOTO 110 291 : * 292 : * Récupération des infos pour la zone IMAIL ou Darcy est définie 293 : * 294 : MELEME = IRIGEL(1,IMAIL) 295 : xMATRI = IRIGEL(4,IMAIL) 296 : SEGACT MELEME 297 : NBDDL = NUM(/1) 298 : NBELEM = NUM(/2) 299 : SEGACT xMATRI 300 : SEGINI MTRAV 301 : * 302 : * Activation des objets necessaires à la prise en compte des forces de volumes 303 : * 304 : IF (IFORC.NE.0) THEN 305 : MCHAM1 = MCHEL1.ICHAML(IMAIL) 306 : SEGACT MCHAM1 307 : MELVA1 = MCHAM1.IELVAL(1) 308 : SEGDES MCHAM1 309 : SEGACT MELVA1 310 : xMATR1 = RI1.IRIGEL(4,IMAIL) 311 : SEGACT xMATR1 312 : ELSE 313 : DO 35 I=1,NBDDL 314 : RFOR(I)=0.D0 315 : 35 CONTINUE 316 : ENDIF 317 : * 318 : C 319 : C Activation du MELVAL du MCHAML Ck|K| 320 : C 321 : IF (IPCK.NE.0) THEN 322 : MCHAML = ICHAML(IMAIL) 323 : SEGACT MCHAML 324 : MELVAL = IELVAL(1) 325 : SEGDES MCHAML 326 : SEGACT MELVAL 327 : ENDIF 328 : * 329 : *------------------------------------------------------- 330 : * BOUCLE SUR LES ELEMENTS DU MAILLAGE ELEMENTAIRE IMAIL 331 : *------------------------------------------------------- 332 : * 333 : DO 100 IEL=1,NBELEM 334 : ITELEM = ITELEM + 1 335 : * 336 : *- Calcul de VECT = DIV * RE , CONSD = VECT * DIV 337 : * 338 : * XMATRI = IMATTT(IEL) 339 : * SEGACT XMATRI 340 : CONSD = 0.D0 341 : DO 45 J=1,NBDDL 342 : VECT(J) = 0.D0 343 : DO 40 I=1,NBDDL 344 : VECT(J) = VECT(J) + RE(I,J,iel) 345 : 40 CONTINUE 346 : CONSD = CONSD + VECT(J) 347 : 45 CONTINUE 348 : * SEGDES XMATRI 349 : * 350 : *- Recuperation de SOUR et de H(n-1) 351 : * 352 : IF (ISOUR.NE.0) VSOUR=MPOVA1.VPOCHA(ITELEM,1) 353 : IF (IFORC.NE.0) THEN 354 : * 355 : *- calcul des flux de forces aux faces de l'element 356 : * 357 : DO 55 IDDL=1,NBDDL 358 : FLFOR(IDDL)= 0.D0 359 : IPOPTS = ICPR(NUM(IDDL,IEL)) 360 : DO 50 I=1,IDIM 361 : FLFOR(IDDL) = FLFOR(IDDL) + MPOVA5.VPOCHA(IPOPTS,I) 362 : $ *MELVA1.VELCHE(IDDL,IEL) * MPOVA4 363 : $ .VPOCHA(ITELEM,I) *MPOVA6.VPOCHA(IPOPTS,1) 364 : 50 CONTINUE 365 : 55 CONTINUE 366 : * 367 : *- Construction du tableau aux faces M.FORCE 368 : * 369 : * XMATR1 = IMATR1.IMATTT(IEL) 370 : * SEGACT XMATR1 371 : DO 65 I=1,NBDDL 372 : RFOR(I)=0.D0 373 : DO 60 J=1,NBDDL 374 : RFOR(I) = RFOR(I) + XMATR1.RE(I,J,iel)*FLFOR(J) 375 : 60 CONTINUE 376 : 65 CONTINUE 377 : * SEGDES XMATR1 378 : ENDIF 379 : * 380 : *- Calcul de p pour l'element IEL 381 : * 382 : VALP = 0.D0 383 : IF (IPCK.EQ.0) THEN 384 : DO 80 I=1,NBDDL 385 : NUMFA = ICPR(NUM(I,IEL)) 386 : VALP = VALP + VECT(I) * (VAA(NUMFA) - RFOR(I)) 387 : 80 CONTINUE 388 : VPOCHA(ITELEM,1) = (VALP + VSOUR) / CONSD 389 : ELSE 390 : BETA = CONSD*DELTAT/(VELCHE(1,IEL)+THETA*CONSD*DELTAT) 391 : DO 90 I=1,NBDDL 392 : NUMFA = ICPR(NUM(I,IEL)) 393 : VALP = VALP + VECT(I) * ( THETA*VAA(NUMFA) 394 : S + (1.D0-THETA)*MPOVA2.VPOCHA(NUMFA,1)) 395 : 90 CONTINUE 396 : VALP = (VALP + VSOUR) * BETA / CONSD 397 : VH = MPOVA3.VPOCHA(ITELEM,1) 398 : VPOCHA(ITELEM,1) = VALP + (1.D0-BETA)*VH 399 : ENDIF 400 : 100 CONTINUE 401 : IF (IFORC.NE.0) SEGDES xMATR1, MELVA1 402 : SEGDES xMATRI, MELEME 403 : IF (IPCK.NE.0) SEGDES MELVAL 404 : SEGSUP MTRAV 405 : 110 CONTINUE 406 : SEGDES MPOVAL 407 : SEGDES MRIGID, IPMAHY 408 : IF (IPCK.NE.0) SEGDES MCHELM, MPOVA2, MPOVA3 409 : IF (ISOUR.NE.0) SEGDES MPOVA1 410 : IF (IFORC.NE.0) 411 : S SEGDES MPOVA4, MPOVA5, MPOVA6, 412 : S MCHEL1, RI1 413 : SEGSUP ICCPR , IORGA 414 : * 415 : RETURN 416 : END 417 : 418 : 419 : 420 : 421 : 422 : 423 : 424 : 425 : 426 :
© Cast3M 2003 - All rights reserved.
Disclaimer