1 : C HYBP SOURCE CHAT 11/03/16 21:25:04 6902 2 : SUBROUTINE HYBP 3 : C----------------------------------------------------------------------- 4 : C Cet operateur permet de calculer la charge en fonction des traces de 5 : C charge dans le cadre d'une formulation variationnelle mixte hybride 6 : C pour les equations de DARCY. 7 : C 8 : C (n) (n-1) -1 t -1 9 : C H = (1 - BETA) * H + (DIV * M1 * DIV) * BETA * 10 : C -1 (n) (n-1) 11 : C ( DIV * M1 * (THETA * TH + (1-THETA) * TH - M * F) + SOUR ) 12 : C 13 : C--------------------------- 14 : C Phrase d'appel (GIBIANE) : 15 : C--------------------------- 16 : C 17 : C CHPO1 = HYBP MMODEL MRIGI1 CHPO2 (TABLE2) (CHPO3) (MRIGI2 CHPO4); 18 : C 19 : C------------------------ 20 : C Operandes et resultat : 21 : C------------------------ 22 : C 23 : C CHPO1 : CHPO centre de nom de composante H. 24 : C MRIGI1 : RIGIDITE masses hybrides élémentaires de sous-type DARCY. 25 : C MRIGI2 : Matrices masses hybrides elementaires de sous-type MASSE. 26 : C MMODEL : MODELE spécifiant la formulation. 27 : C CHPO2 : CHPO face de nom de composante TH. 28 : C TABLE2 : TABLE DARCY_TRANSITOIRE a inclure en transitoire; contient 29 : C Theta, le pas de temps, Ck|K|, TH(n-1) et H(n-1). 30 : C CHPO3 : CHPO centre de nom de composante SOUR. 31 : C CHPO4 : CHPO centre de composante FX,FY(,FZ), 32 : C densite de force moyenne par élément. 33 : C 34 : C----------------------------------------------------------------------- 35 : C 36 : C Langage : ESOPE + FORTRAN77 37 : C 38 : C Auteurs : 09/93 F.DABBENE - Cas permanent 39 : C 09/94 X.NOUVELLON - Extension au cas transitoire 40 : C 02/96 L.V.BENET - Prise en compte de forces de volume 41 : C 42 : C----------------------------------------------------------------------- 43 : IMPLICIT INTEGER(I-N) 44 : IMPLICIT REAL*8 (A-H,O-Z) 45 : * 46 : -INC CCOPTIO 47 : -INC SMCHAML 48 : -INC SMCHPOI 49 : -INC SMELEME 50 : -INC SMMODEL 51 : -INC SMRIGID 52 : -INC SMTABLE 53 : * 54 : SEGMENT IPMAHY 55 : INTEGER MAHYBR(NSOUS) 56 : ENDSEGMENT 57 : * 58 : LOGICAL LOGRE,LOGIN 59 : CHARACTER*4 NOMTOT(1) 60 : CHARACTER*4 NOMTO3(3) 61 : CHARACTER*6 CHAR6 62 : CHARACTER*8 TAPIND,TYPOBJ,CHARIN,CHARRE,LETYPE 63 : * 64 : * Initialisations 65 : * 66 : IVALIN = 0 67 : XVALIN = 0.D0 68 : LOGIN = .TRUE. 69 : IOBIN = 0 70 : TAPIND = 'MOT ' 71 : * 72 : * Lecture du MMODEL et test de la formulation 73 : * 73 : CALL LIROBJ('MMODEL',IPMODE,1,IRET) 74 : IF (IERR.NE.0) RETURN 75 : MMODEL = IPMODE 76 : * 77 : * Lecture de la TABLE domaine 78 : * 79 : IPTABL = 0 79 : C CALL LIRTAB('DOMAINE',IPTABL,1,IRET) 79 : CALL LEKMOD(MMODEL,IPTABL,IRET) 80 : IF (IERR.NE.0) RETURN 81 : CHARIN = 'MAILLAGE' 81 : CALL LEKTAB(IPTABL,CHARIN,IOBRE) 82 : IF (IERR.NE.0) RETURN 83 : IPGEOM = IOBRE 84 : CHARIN = 'ELTFA' 84 : CALL LEKTAB(IPTABL,CHARIN,IOBRE) 85 : IF (IERR.NE.0) RETURN 86 : IELTFA = IOBRE 87 : CHARIN = 'CENTRE' 87 : CALL LEKTAB(IPTABL,CHARIN,IOBRE) 88 : IF (IERR.NE.0) RETURN 89 : ICENTR = IOBRE 90 : CHARIN = 'FACE' 90 : CALL LEKTAB(IPTABL,CHARIN,IOBRE) 91 : IF (IERR.NE.0) RETURN 92 : IPFACE = IOBRE 93 : * 94 : * récup. MCHAML orientation normale 95 : * 95 : CALL LEKTAB(IPTABL,'XXNORMAE',IORIE) 96 : IF (IERR.NE.0) RETURN 97 : * 98 : * récup. CHPO orientation normale 99 : * 99 : CALL LEKTAB(IPTABL,'XXNORMAF',INORM) 100 : IF (IERR.NE.0) RETURN 101 : * 102 : * récup. CHPO surface des faces 103 : * 103 : CALL LEKTAB(IPTABL,'XXSURFAC',ISURF) 104 : IF (IERR.NE.0) RETURN 105 : * 106 : * Lecture de RIGIDITE obligatoire 107 : * 107 : CALL LIROBJ('RIGIDITE',IPRIGI,1,IRET) 108 : IF (IERR.NE.0) RETURN 109 : RI1 = IPRIGI 110 : * 111 : * Lecture du CHAMPOINT obligatoire 112 : * 112 : CALL LIROBJ('CHPOINT',ITPN,1,IRET) 113 : IF (IERR.NE.0) RETURN 114 : * 115 : * Lecture eventuelle du 1er CHAMPOINT optionnel 116 : * 117 : ICHP1 = 0 117 : CALL LIROBJ('CHPOINT',ICHP1,0,IRET1) 118 : IF (IERR.NE.0) RETURN 119 : * 120 : * Lecture eventuelle du 2me CHAMPOINT optionnel 121 : * 122 : ICHP2 = 0 122 : CALL LIROBJ('CHPOINT',ICHP2,0,IRET2) 123 : IF (IERR.NE.0) RETURN 124 : * 125 : * Lecture de la RIGIDITE complementaire 126 : * 127 : IPRIGI = 0 127 : CALL LIROBJ('RIGIDITE',IPRIGI,0,IRETR) 128 : IF (IERR.NE.0) RETURN 129 : RI2 = IPRIGI 130 : * 131 : * Lecture éventuelle de la TABLE Darcy_transitoire 132 : * 133 : IPTABL = 0 133 : CALL LIRTAB('DARCY_TRANSITOIRE',IPTABL,0,IRET) 134 : IF (IERR.NE.0) RETURN 135 : IF (IRET.EQ.1) THEN 136 : TYPOBJ = 'FLOTTANT' 136 : CALL ACCTAB(IPTABL,TAPIND,IVALIN,XVALIN,'THETA',LOGIN,IOBIN, 137 : S TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE) 138 : IF (IERR.NE.0) RETURN 139 : THETA = XVALRE 140 : THEMIN = -1.D-12 141 : THEMAX = 1.D0 - THEMIN 142 : IF ((THETA.LT.THEMIN).OR.(THETA.GT.THEMAX)) THEN 143 : REAERR(1) = REAL(THETA) 144 : REAERR(2) = REAL(0.D0) 145 : REAERR(3) = REAL(1.D0) 146 : MOTERR(1:8) = ' THETA ' 146 : CALL ERREUR(42) 147 : RETURN 148 : ENDIF 149 : IF (THETA.LT.0.D0) THETA=0.D0 150 : IF (THETA.GT.1.D0) THETA=1.D0 150 : CALL ACCTAB(IPTABL,TAPIND,IVALIN,XVALIN,'PAS',LOGIN,IOBIN, 151 : S TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE) 152 : IF (IERR.NE.0) RETURN 153 : DELTAT = XVALRE 154 : IF (DELTAT.LT.0.D0) THEN 155 : REAERR(1) = REAL(DELTAT) 156 : REAERR(2) = REAL(0.D0) 157 : MOTERR(1:8) = ' DELTAT ' 157 : CALL ERREUR(41) 158 : RETURN 159 : ENDIF 160 : TYPOBJ = 'MCHAML ' 160 : CALL ACCTAB(IPTABL,TAPIND,IVALIN,XVALIN,'SURF',LOGIN,IOBIN, 161 : S TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE) 162 : IF (IERR.NE.0) RETURN 163 : IPCK = IOBRE 164 : TYPOBJ = 'CHPOINT ' 164 : CALL ACCTAB(IPTABL,TAPIND,IVALIN,XVALIN,'TRACE',LOGIN,IOBIN, 165 : S TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE) 166 : IF (IERR.NE.0) RETURN 167 : ITP = IOBRE 167 : CALL ACCTAB(IPTABL,TAPIND,IVALIN,XVALIN,'CHARGE',LOGIN,IOBIN, 168 : S TYPOBJ,IVALRE,XVALRE,CHARRE,LOGRE,IOBRE) 169 : IF (IERR.NE.0) RETURN 170 : IH = IOBRE 171 : ELSE 172 : THETA = 0.D0 173 : DELTAT = 0.D0 174 : IPCK = 0 175 : ITP = 0 176 : IH = 0 177 : ENDIF 178 : * 179 : * Test CHAMPOINT centre H 180 : * 181 : IF (IH.NE.0) THEN 182 : NBCOMP = 1 183 : INDIC = 1 184 : NOMTOT(1) = 'H' 184 : CALL QUEPOI(IH,ICENTR,INDIC,NBCOMP,NOMTOT) 185 : IF (IERR.NE.0) RETURN 186 : ENDIF 187 : * (n-1) 188 : * Test du CHAMPOINT face TH 189 : * 190 : IF (ITP.NE.0) THEN 191 : NBCOMP = 1 192 : INDIC = 1 193 : NOMTOT(1) = 'TH' 193 : CALL QUEPOI(ITP,IPFACE,INDIC,NBCOMP,NOMTOT) 194 : IF (IERR.NE.0) RETURN 195 : ENDIF 196 : * 197 : * Test de la formulation contenu dans MMODEL 198 : * 199 : SEGACT MMODEL 200 : NSOUS = KMODEL(/1) 201 : SEGINI IPMAHY 202 : IDARCY = 0 203 : IPT1=IPGEOM 204 : IPT2=IPGEOM 205 : DO 10 ISOUS=1,NSOUS 206 : IF(NSOUS.GT.1)THEN 207 : SEGACT IPT2 208 : IPT1= IPT2.LISOUS(ISOUS) 209 : SEGDES IPT2 210 : ENDIF 211 : IMODEL = KMODEL(ISOUS) 212 : SEGACT IMODEL 213 : LETYPE = FORMOD(1) 214 : IF (LETYPE.EQ.'DARCY') THEN 215 : IDARCY = IDARCY + 1 216 : MAHYBR(ISOUS) = IPT1 217 : ENDIF 218 : SEGDES IMODEL 219 : 10 CONTINUE 220 : SEGDES MMODEL 221 : IF (IDARCY.EQ.0) THEN 222 : MOTERR = LETYPE 222 : CALL ERREUR(193) 223 : GOTO 100 224 : ENDIF 225 : * 226 : * Tri et test des CHAMPOINTs ARGUMENTS 227 : * 228 : ISOUR=0 229 : IFORC=0 230 : IF(IRET1.EQ.1)THEN 231 : NBCOMP = 0 232 : INDIC = 1 233 : NOMTO3(1) = ' ' 234 : NOMTO3(2) = ' ' 235 : NOMTO3(3) = ' ' 235 : CALL QUEPOI(ICHP1,ICENTR,INDIC,NBCOMP,NOMTO3) 236 : IF (IERR.NE.0) RETURN 237 : IF(NOMTO3(1).EQ.'SOUR')THEN 238 : ISOUR=ICHP1 239 : ELSE IF(NOMTO3(1).EQ.'FX')THEN 240 : IFORC=ICHP1 241 : ENDIF 242 : IF(IRET2.EQ.1)THEN 243 : NBCOMP = 0 244 : INDIC = 1 245 : NOMTO3(1) = ' ' 246 : NOMTO3(2) = ' ' 247 : NOMTO3(3) = ' ' 247 : CALL QUEPOI(ICHP2,ICENTR,INDIC,NBCOMP,NOMTO3) 248 : IF (IERR.NE.0) RETURN 249 : IF(NOMTO3(1).EQ.'SOUR')THEN 250 : ISOUR=ICHP2 251 : ELSE IF(NOMTO3(1).EQ.'FX')THEN 252 : IFORC=ICHP2 253 : ENDIF 254 : ENDIF 255 : ENDIF 256 : * 257 : * Récuperation des pointeurs ELTFA pour les zones ou DARCY est defini 258 : * 259 : MELEME = IELTFA 260 : SEGACT MELEME 261 : LZONES = LISOUS(/1) 262 : IF (LZONES.EQ.0) LZONES = 1 263 : IPT4 = IPGEOM 264 : SEGACT IPT4 265 : DO 30 ISOUS=1,NSOUS 266 : IMAMEL = MAHYBR(ISOUS) 267 : IF (IMAMEL.NE.0) THEN 268 : DO 20 ISZ=1,LZONES 269 : IF (LZONES.EQ.1) THEN 270 : IPT2 = IPT4 271 : IPT3 = MELEME 272 : ELSE 273 : IPT2 = IPT4.LISOUS(ISZ) 274 : IPT3 = LISOUS(ISZ) 275 : ENDIF 276 : IF (IPT2.EQ.IMAMEL) THEN 277 : MAHYBR(ISOUS) = IPT3 278 : GOTO 30 279 : ENDIF 280 : 20 CONTINUE 281 : IF (IMAMEL.EQ.MAHYBR(ISOUS)) THEN 282 : MOTERR(1:8) = ' MODELE ' 283 : MOTERR(9:16) = ' ELTFA ' 284 : INTERR(1) = ISOUS 284 : CALL ERREUR(698) 285 : SEGDES IPT4 286 : SEGDES MELEME 287 : GOTO 100 288 : ENDIF 289 : ENDIF 290 : 30 CONTINUE 291 : SEGDES IPT4 292 : SEGDES MELEME 293 : * 294 : * Test et tri du sous-type des matrices de rigiditée récupérées 295 : * 296 : IPFORC=0 297 : IPDARC=0 298 : SEGACT RI1 299 : LETYPE = RI1.MTYMAT 300 : IF(LETYPE.EQ.'DARCY')THEN 301 : IPDARC=RI1 302 : ELSE IF(LETYPE.EQ.'MASSE')THEN 303 : IPFORC=RI1 304 : ENDIF 305 : IF(RI2.NE.0)THEN 306 : SEGACT RI2 307 : LETYPE = RI2.MTYMAT 308 : IF(LETYPE.EQ.'MASSE')THEN 309 : IPFORC=RI2 310 : ELSE IF(LETYPE.EQ.'DARCY')THEN 311 : IPDARC=RI2 312 : ENDIF 313 : IF(IPFORC.EQ.0)THEN 314 : MOTERR(1:8) = 'RIGIDITE' 315 : MOTERR(9:16) = 'MASSE' 315 : CALL ERREUR(79) 316 : SEGDES RI1 317 : SEGDES RI2 318 : GOTO 100 319 : ELSE IF(IFORC.EQ.0)THEN 320 : MOTERR(1:8) = 'CHPOINT' 321 : MOTERR(9:16) = 'FORCE' 321 : CALL ERREUR(79) 322 : SEGDES RI1 323 : SEGDES RI2 324 : GOTO 100 325 : ENDIF 326 : ENDIF 327 : IF(IPDARC.EQ.0)THEN 328 : MOTERR(1:8) = 'RIGIDITE' 329 : MOTERR(9:16) = 'DARCY' 329 : CALL ERREUR(79) 330 : SEGDES RI1 331 : IF(IPFORC.NE.0)SEGDES RI2 332 : GOTO 100 333 : ENDIF 334 : * 335 : * Controle des pointeurs de MELEME de la rigidité 336 : * 337 : DO 40 ISOUS=1,NSOUS 338 : IMAMEL = MAHYBR(ISOUS) 339 : IF (IMAMEL.NE.0) THEN 340 : RI1=IPDARC 341 : IPTEST = RI1.IRIGEL(1,ISOUS) 341 : CALL CMPMAI(IMAMEL,IPTEST,LOGRE) 342 : IF (.NOT. LOGRE) THEN 343 : MOTERR(1:8) = 'DARCY' 344 : MOTERR(9:16) = ' ELTFA ' 345 : INTERR(1) = ISOUS 345 : CALL ERREUR(698) 346 : SEGDES RI1 347 : GOTO 100 348 : ENDIF 349 : IF(IPFORC.NE.0)THEN 350 : RI2=IPFORC 351 : IPTEST = RI2.IRIGEL(1,ISOUS) 351 : CALL CMPMAI(IMAMEL,IPTEST,LOGRE) 352 : IF (.NOT. LOGRE) THEN 353 : MOTERR(1:8) = 'MASSE' 354 : MOTERR(9:16) = ' ELTFA ' 355 : INTERR(1) = ISOUS 355 : CALL ERREUR(698) 356 : SEGDES RI1 357 : SEGDES RI2 358 : GOTO 100 359 : ENDIF 360 : ENDIF 361 : ENDIF 362 : 40 CONTINUE 363 : SEGDES RI1 364 : IF(IPFORC.NE.0)SEGDES RI2 365 : * 366 : * Vérification du support du MCHAML Ck|K| % au MMODEL 367 : * Controle du nombre de composantes reelles 368 : * 369 : IF (IPCK.NE.0) THEN 369 : IF(IPFORC.NE.0)CALL ERREUR(34) 370 : ISUP = 2 371 : ICOND = 1 371 : CALL QUESUP(IPMODE,IPCK,ISUP,ICOND,IRET,IRET2) 372 : IF (IRET.GT.1) GOTO 100 373 : MCHELM = IPCK 374 : SEGACT MCHELM 375 : DO 50 ISOUS=1,NSOUS 376 : IF (MAHYBR(ISOUS).NE.0) THEN 377 : MCHAML = ICHAML(ISOUS) 378 : SEGACT MCHAML 379 : N2 = TYPCHE(/2) 380 : IF (N2.GT.1) THEN 380 : CALL ERREUR(320) 381 : SEGDES MCHELM 382 : SEGDES MCHAML 383 : GOTO 100 384 : ENDIF 385 : CHAR6 = TYPCHE(1)(1:6) 386 : IF (CHAR6.NE.'REAL*8') THEN 387 : MOTERR(1:8) = NOMCHE(1) 387 : CALL ERREUR(679) 388 : SEGDES MCHELM 389 : SEGDES MCHAML 390 : GOTO 100 391 : ENDIF 392 : SEGDES MCHAML 393 : ENDIF 394 : 50 CONTINUE 395 : SEGDES MCHELM 396 : ENDIF 397 : * (n) 398 : * Test du CHAMPOINT face TH 399 : * On ne teste que la présence d'une composante de nom TH 400 : * 401 : MCHPO1 = ITPN 402 : SEGACT MCHPO1 403 : NSOUP1 = MCHPO1.IPCHP(/1) 404 : DO 70 ISOUS=1,NSOUP1 405 : MSOUPO = MCHPO1.IPCHP(ISOUS) 406 : SEGACT MSOUPO 407 : NC = NOCOMP(/2) 408 : DO 60 J=1,NC 409 : IF (NOCOMP(J).EQ.'TH ') GOTO 80 410 : 60 CONTINUE 411 : SEGDES MSOUPO 412 : 70 CONTINUE 413 : SEGDES MCHPO1 414 : MOTERR(1:4) = ' TH ' 414 : CALL ERREUR(181) 415 : GOTO 100 416 : 80 CONTINUE 417 : SEGDES MSOUPO 418 : SEGDES MCHPO1 419 : * 420 : * Construction du CHPO des charges 421 : * 422 : SEGDES IPMAHY 422 : CALL HYBP1(IPMAHY,ICENTR,IPFACE,IPDARC,ITPN,INORM,IORIE, 423 : S ISURF,THETA,DELTAT,IPCK,ITP,IH,ISOUR,IPFORC,IFORC,ICHP1) 423 : IF (ICHP1.NE.0) CALL ECROBJ('CHPOINT',ICHP1) 424 : * 425 : * Ménage 426 : * 427 : 100 CONTINUE 428 : SEGSUP IPMAHY 429 : RETURN 430 : END 431 : 432 : 433 : 434 : 435 : 436 : 437 : 438 : 439 : 440 : 441 : 442 : 443 : 444 : 445 :
© Cast3M 2003 - All rights reserved.
Disclaimer