1 : C PERMAB SOURCE AM 12/01/26 21:15:01 7238 2 : SUBROUTINE PERMAB (MODORI,MCHORI,IPRIGI,IRET) 3 : *______________________________________________________________________ 4 : * 5 : * OPERATEUR PERMEABILITE (MILIEUX POREUX) APPELE PAR PERMEA 6 : * 7 : * Entrees : 8 : * --------- 9 : * 10 : * MODORI Pointeur sur un MMODEL 11 : * MCHORI Pointeur sur un MCHAML de materiau 12 : * 13 : * Sorties : 14 : * --------- 15 : * 16 : * IPRIGI Pointeur sur un objet RIGIDITE de permeabilite 17 : * IRET =1 ou 0 suivant succes ou non 18 : * 19 : * Passage aux nouveaux CHAMELEMs par jm CAMPENON le 07/91 20 : *______________________________________________________________________ 21 : * 22 : IMPLICIT INTEGER(I-N) 23 : IMPLICIT REAL*8(A-H,O-Z) 24 : * 25 : -INC CCHAMP 26 : -INC CCOPTIO 27 : -INC CCGEOME 28 : -INC SMRIGID 29 : -INC SMCHAML 30 : -INC SMELEME 31 : -INC SMCOORD 32 : -INC SMINTE 33 : -INC SMMODEL 34 : * 35 : SEGMENT WRK1 36 : REAL*8 REL(LRE,LRE),XE(3,NBBB) 37 : ENDSEGMENT 38 : * 39 : SEGMENT WRK2 40 : REAL*8 SHPWRK(6,NBNO),BGENE(NSTB,LRE) 41 : ENDSEGMENT 42 : * 43 : SEGMENT WRK3 44 : REAL*8 XGENE(NSTN,LRN) 45 : ENDSEGMENT 46 : * 47 : SEGMENT WRK4 48 : REAL*8 XLOC(IDIM,IDIM),XGLOB(IDIM,IDIM),TXR(IDIM,IDIM) 49 : REAL*8 VALMAT(NMATT) 50 : REAL*8 PMAT(NSTB,NSTB),PMAT1(IDIM,IDIM),PMAT2(IDIM,IDIM) 51 : ENDSEGMENT 52 : * 53 : SEGMENT WRK5 54 : REAL*8 BPSS(3,3),XEL(3,NBBB) 55 : REAL*8 XNTH(LRN,LRN),XNTB(LRN,LRN),XNTT(LRN) 56 : ENDSEGMENT 57 : * 58 : SEGMENT WRK6 59 : REAL*8 PKK(NSTB,NSTB) 60 : ENDSEGMENT 61 : * 62 : SEGMENT INFO 63 : INTEGER INFELL(JG) 64 : ENDSEGMENT 65 : * 66 : SEGMENT NOTYPE 67 : CHARACTER*16 TYPE(NBTYPE) 68 : ENDSEGMENT 69 : * 70 : SEGMENT MPTVAL 71 : INTEGER IPOS(NS),NSOF(NS) 72 : INTEGER IVAL(NCOSOU) 73 : CHARACTER*16 TYVAL(NCOSOU) 74 : ENDSEGMENT 75 : * 76 : 77 : INTEGER OOOVAL 78 : 79 : CHARACTER*8 CMATE 80 : CHARACTER*(NCONCH) CONM 81 : 82 : 83 : * INTTYP correspond au type de points d'integration utilise 84 : PARAMETER ( INTTYP=3 ) 85 : 86 : PARAMETER (NINF=3) 87 : INTEGER INFOS(NINF) 88 : LOGICAL lsupfo,lsupdp 89 : * 90 : NHRM=NIFOUR 91 : IRET = 0 92 : IPRIGI = 0 93 : * 94 : * Reduction du modele a la formulation poreuse 95 : * 96 : MMODE1 = MODORI 97 : SEGINI,MMODEL=MMODE1 98 : NSOUS = KMODEL(/1) 99 : N1 = 0 100 : DO isous = 1, NSOUS 101 : IMODEL = KMODEL(isous) 102 : SEGACT,IMODEL 103 : IF (FORMOD(1).EQ.'POREUX') THEN 104 : N1 = N1 + 1 105 : KMODEL(N1) = IMODEL 106 : ELSE 107 : SEGDES,IMODEL 108 : ENDIF 109 : ENDDO 110 : IF (N1.NE.NSOUS) SEGADJ,MMODEL 111 : IPMODL = MMODEL 112 : NSOUS = N1 113 : IF (NSOUS.LE.0) THEN 114 : MOTERR(1:8) = 'MMODEL ' 115 : INTERR(1) = MODORI 115 : CALL ERREUR(356) 116 : GOTO 9991 117 : ENDIF 118 : 119 : * 120 : * Reduction du champ au modele precedemment reduit 121 : * 122 : MCHELM = MCHORI 123 : SEGACT,MCHELM 124 : IF (TITCHE(1:8).NE.'CARACTER') THEN 125 : SEGDES,MCHELM 126 : MOTERR(1:16) = 'CARACTERISTIQUES' 126 : CALL ERREUR(291) 127 : GOTO 9991 128 : ENDIF 128 : CALL REDUAF(MCHORI,IPMODL,IPCHE1,0,IRET,KERRE) 129 : SEGDES,MCHELM 130 : IF (IRET.NE.1) THEN 130 : CALL ERREUR(KERRE) 131 : GOTO 9991 132 : ENDIF 133 : * 134 : * Verification du lieu support du MCHAML de materiau 135 : * 136 : ISUP=0 136 : CALL QUESUP(IPMODL,IPCHE1,INTTYP,0,ISUP,IRETMA) 137 : IF (ISUP.GT.1) GO TO 9991 138 : * 139 : * Activation du MMODEL 140 : * 141 : MMODEL=IPMODL 142 : SEGACT MMODEL 143 : NSOUS=KMODEL(/1) 144 : * 145 : * ON FABRIQUE LES MATRICES UNIQUEMENT POUR LES ZONES 146 : * DE MILIEU POREUX. 147 : * 148 : NRIGEL=NSOUS 149 : * 150 : * Initialisation du chapeau de l'objet RIGIDITE 151 : * 152 : SEGINI MRIGID 153 : ICHOLE=0 154 : IMGEO1=0 155 : IMGEO2=0 156 : IFORIG=IFOUR 157 : MTYMAT='PERMEABI' 158 : * 159 : * BOUCLE SUR LES SOUS ZONES DU MODELE 160 : * 161 : ISORI=0 162 : DO 500 ISOUS=1,NSOUS 163 : * 164 : * On recupere l'information generale 165 : * 166 : IMODEL=KMODEL(ISOUS) 167 : SEGACT IMODEL 168 : 169 : IF(FORMOD(1).NE.'POREUX') THEN 169 : CALL ERREUR(19) 170 : GO TO 9999 171 : ENDIF 172 : * 173 : * Traitement du modele 174 : * 175 : MELE=NEFMOD 176 : IPMAIL=IMAMOD 177 : CONM =CONMOD 178 : 179 : * 179 : CALL IDENT(IPMAIL,CONM,IPCHE1,0,INFOS,IRTD) 180 : IF (IRTD.EQ.0) GOTO 9999 181 : * 182 : * Nature du materiau 183 : * 184 : CMATE = CMATEE 185 : MATE = IMATEE 186 : INAT = INATUU 187 : * 188 : * Information sur l'element fini 189 : * 190 : IF (INFMOD(/1).LT.2+INTTYP) THEN 190 : CALL ELQUOI(MELE,0,INTTYP,IPINF,IMODEL) 191 : IF (IERR.NE.0) GO TO 9999 192 : INFO=IPINF 193 : MFR =INFELL(13) 194 : IELE =INFELL(14) 195 : IPORE=INFELL(8) 196 : MINTE=INFELL(11) 197 : segsup info 198 : else 199 : MFR =INFELE(13) 200 : IELE =INFELE(14) 201 : IPORE=INFELE(8) 202 : MINTE=infmod(5) 203 : endif 204 : IPMINT=MINTE 205 : * 206 : * Si necessaire PARTITIONNEMENT du segment XMATRI 207 : * 208 : IPT1=IPMAIL 209 : SEGACT,IPT1 210 : NBNN1 =IPT1.NUM(/1) 211 : NBELE1=IPT1.NUM(/2) 212 : * 213 : IF(MFR.EQ.33) THEN 214 : IDECAP=1 215 : ELSE IF(MFR.EQ.57) THEN 216 : IDECAP=2 217 : ELSE IF(MFR.EQ.59) THEN 218 : * 219 : IDECAP=3 220 : ENDIF 221 : LRE=NBNNE(IELE)* IDECAP 222 : LR1=NBNNE(IELE) 223 : * 224 : LTRK=OOOVAL(1,4) 225 : IF (LTRK.EQ.0) LTRK=OOOVAL(1,1) 226 : * Ajout a la taille en mots de la matrice des infos du segment 227 : LSEG=LRE*LRE*NBELE1 + 16 228 : NBLPRT=(LSEG-1)/LTRK+1 229 : NBLMAX=(NBELE1-1)/NBLPRT+1 230 : NBLPRT=(NBELE1-1)/NBLMAX+1 231 : * write(ioimp,*) ' PERMAB nblprt nblmax ',NBLPRT,NBLMAX,NBELE1 232 : MELEME=IPT1 233 : 234 : * BOUCLE DE PARTITIONNEMENT DU SEGMENT XMATRI 235 : 236 : DO 5000 IPRT = 1,NBLPRT 237 : ISORI= ISORI+1 238 : IF (ISORI.GT.IRIGEL(/2)) THEN 239 : NRIGEL=ISORI 240 : SEGADJ,MRIGID 241 : ENDIF 242 : IF (NBLPRT.GT.1) THEN 243 : JPRT=(IPRT-1)*NBLMAX 244 : SEGACT,IPT1 245 : NBSOUS=0 246 : NBREF=0 247 : NBNN=NBNN1 248 : NBELEM=MIN(NBLMAX,NBELE1-JPRT) 249 : * write(ioimp,*) ' creation segment ',nbnn,nbelem 250 : SEGINI,MELEME 251 : ITYPEL=IPT1.ITYPEL 252 : DO I=1,NBELEM 253 : IB=I+JPRT 254 : DO J=1,NBNN 255 : NUM(J,I)=IPT1.NUM(J,IB) 256 : ENDDO 257 : ICOLOR(I)=IPT1.ICOLOR(I) 258 : ENDDO 259 : ENDIF 260 : IPMAIL=MELEME 261 : * Fin du traitement particulier en cas de PARTITIONNEMENT du XMATRI 262 : * Quelques initialisations suite au partionnement 263 : IPDES = 0 264 : * 265 : NMATR = 0 266 : NMATF = 0 267 : IVAMAT = 0 268 : NCARA = 0 269 : NCARF = 0 270 : IVACAR = 0 271 : 272 : * Activation du MELEME support des rigidites 273 : MELEME=IPMAIL 274 : SEGACT,MELEME 275 : NBNN =NUM(/1) 276 : NBELEM=NUM(/2) 277 : 278 : NLIGRP = LRE 279 : NLIGRD = LRE 280 : SEGINI DESCR 281 : if(lnomid(1).ne.0) then 282 : nomid=lnomid(1) 283 : segact nomid 284 : modepl=nomid 285 : ndepl=lesobl(/2) 286 : ndum=lesfac(/2) 287 : lsupdp=.false. 288 : else 289 : lsupdp=.true. 289 : CALL IDDEPL(MFR,IFOUR,MODEPL,NDEPL,NDUM) 290 : endif 291 : if(lnomid(2).ne.0) then 292 : nomid=lnomid(2) 293 : segact nomid 294 : moforc=nomid 295 : nforc=lesobl(/2) 296 : lsupfo=.false. 297 : else 298 : lsupfo=.true. 298 : CALL IDFORC(MFR,IFOUR,MOFORC,NFORC,NDUM) 299 : endif 300 : * 301 : ID=1 302 : NOMID=MODEPL 303 : SEGACT NOMID 304 : NCP=LESOBL(/2) 305 : NOMID=MOFORC 306 : SEGACT NOMID 307 : * 308 : IF (MFR.EQ.33) THEN 309 : 310 : DO 4005 IB=1,NBSOM(IELE) 311 : NOMID=MODEPL 312 : LISINC(ID)=LESOBL(NCP) 313 : NOMID=MOFORC 314 : LISDUA(ID)=LESOBL(NCP) 315 : NOELEP(ID)=IBSOM(NSPOS(IELE)+IB-1) 316 : NOELED(ID)=IBSOM(NSPOS(IELE)+IB-1) 317 : ID=ID+1 318 : 4005 CONTINUE 319 : * 320 : IF (MELE.GE.108.AND.MELE.LE.110) THEN 321 : * 322 : LR1=(3*LRE-IPORE)/2 323 : DO 4008 INOEUD=LR1+1,LRE 324 : NOMID=MODEPL 325 : LISINC(ID)=LESOBL(NCP) 326 : NOMID=MOFORC 327 : LISDUA(ID)=LESOBL(NCP) 328 : NOELEP(ID)=INOEUD 329 : NOELED(ID)=INOEUD 330 : ID=ID+1 331 : 4008 CONTINUE 332 : END IF 333 : * 334 : DO 4006 IB=1,LR1 335 : DO 4007 INSOM=1,NBSOM(IELE) 336 : IF(IB.EQ.IBSOM(NSPOS(IELE)+INSOM-1)) GO TO 4006 337 : 4007 CONTINUE 338 : NOMID=MODEPL 339 : LISINC(ID)=LESOBL(NCP) 340 : NOMID=MOFORC 341 : LISDUA(ID)=LESOBL(NCP) 342 : NOELEP(ID)=IB 343 : NOELED(ID)=IB 344 : ID=ID+1 345 : 4006 CONTINUE 346 : * 347 : ELSE IF (MFR.EQ.57.OR.MFR.EQ.59) THEN 348 : 349 : DO 4205 IPR=1,IDECAP 350 : NCPDEC=NCP-IDECAP+IPR 351 : * 352 : DO 4105 IB=1,NBSOM(IELE) 353 : NOMID=MODEPL 354 : LISINC(ID)=LESOBL(NCPDEC) 355 : NOMID=MOFORC 356 : LISDUA(ID)=LESOBL(NCPDEC) 357 : NOELEP(ID)=IBSOM(NSPOS(IELE)+IB-1) 358 : NOELED(ID)=IBSOM(NSPOS(IELE)+IB-1) 359 : ID=ID+1 360 : 4105 CONTINUE 361 : * 362 : * JOINTS : A VOIR PLUS TARD 363 : * 364 : DO 4106 IB=1,LR1 365 : DO 4107 INSOM=1,NBSOM(IELE) 366 : IF(IB.EQ.IBSOM(NSPOS(IELE)+INSOM-1)) GO TO 4106 367 : 4107 CONTINUE 368 : NOMID=MODEPL 369 : LISINC(ID)=LESOBL(NCPDEC) 370 : NOMID=MOFORC 371 : LISDUA(ID)=LESOBL(NCPDEC) 372 : NOELEP(ID)=IB 373 : NOELED(ID)=IB 374 : ID=ID+1 375 : 4106 CONTINUE 376 : 377 : 4205 CONTINUE 378 : * 379 : ENDIF 380 : 381 : NOMID =MOFORC 382 : if(lsupfo)SEGSUP NOMID 383 : NOMID =MODEPL 384 : if(lsupdp)SEGSUP NOMID 385 : IPDES=DESCR 386 : SEGDES DESCR 387 : * 388 : * Initialisation de MINTE 389 : * 390 : SEGACT MINTE 391 : NBPGAU=POIGAU(/1) 392 : * 393 : 394 : CCCCCC LVAL=(LRE*(LRE+1))/2 395 : NELRIG=NBELEM 396 : SEGINI xMATRI 397 : * 398 : * Verification de la presence des composantes pour le materiau 399 : * 400 : NBROBL=0 401 : NBRFAC=0 402 : * cas isotrope 403 : IF (MATE.EQ.1) THEN 404 : * 405 : IF (MELE.GE.79.AND.MELE.LE.83) THEN 406 : NBROBL=2 407 : SEGINI NOMID 408 : LESOBL(1)='PERM' 409 : LESOBL(2)='VISC' 410 : ELSE IF (MELE.GE.108.AND.MELE.LE.110) THEN 411 : NBROBL=4 412 : SEGINI NOMID 413 : LESOBL(1)='PERT' 414 : LESOBL(2)='PERH' 415 : LESOBL(3)='PERB' 416 : LESOBL(4)='VISC' 417 : ELSE IF (MELE.GE.173.AND.MELE.LE.177) THEN 418 : NBROBL=4 419 : SEGINI NOMID 420 : LESOBL(1)='PK11' 421 : LESOBL(2)='PK12' 422 : LESOBL(3)='PK21' 423 : LESOBL(4)='PK22' 424 : ELSE IF (MELE.GE.178.AND.MELE.LE.182) THEN 425 : NBROBL=9 426 : SEGINI NOMID 427 : LESOBL(1)='PK11' 428 : LESOBL(2)='PK12' 429 : LESOBL(3)='PK13' 430 : LESOBL(4)='PK21' 431 : LESOBL(5)='PK22' 432 : LESOBL(6)='PK23' 433 : LESOBL(7)='PK31' 434 : LESOBL(8)='PK32' 435 : LESOBL(9)='PK33' 436 : ENDIF 437 : * cas orthotrope 438 : ELSE IF (MATE.EQ.2) THEN 439 : IF (IDIM.EQ.3) THEN 440 : NBROBL=10 441 : SEGINI NOMID 442 : LESOBL(1)='PER1' 443 : LESOBL(2)='PER2' 444 : LESOBL(3)='PER3' 445 : LESOBL(4)='VISC' 446 : LESOBL(5)='V1X ' 447 : LESOBL(6)='V1Y ' 448 : LESOBL(7)='V1Z ' 449 : LESOBL(8)='V2X ' 450 : LESOBL(9)='V2Y ' 451 : LESOBL(10)='V2Z ' 452 : ELSE IF(IDIM.EQ.2) THEN 453 : IF (IFOUR.LE.0) THEN 454 : NBROBL=5 455 : SEGINI NOMID 456 : LESOBL(1)='PER1' 457 : LESOBL(2)='PER2' 458 : LESOBL(3)='VISC' 459 : LESOBL(4)='V1X ' 460 : LESOBL(5)='V1Y ' 461 : ELSE IF (IFOUR.EQ.1) THEN 462 : NBROBL=6 463 : SEGINI NOMID 464 : LESOBL(1)='PER1' 465 : LESOBL(2)='PER2' 466 : LESOBL(3)='PER3' 467 : LESOBL(4)='VISC' 468 : LESOBL(5)='V1X ' 469 : LESOBL(6)='V1Y ' 470 : ENDIF 471 : ENDIF 472 : * cas anisotrope 473 : ELSE IF (MATE.EQ.3)THEN 474 : IF(IDIM.EQ.3)THEN 475 : NBROBL=13 476 : SEGINI NOMID 477 : LESOBL(1)='PER1' 478 : LESOBL(2)='PER2' 479 : LESOBL(3)='PER3' 480 : LESOBL(4)='PE12' 481 : LESOBL(5)='PE13' 482 : LESOBL(6)='PE23' 483 : LESOBL(7)='VISC' 484 : LESOBL(8)='V1X ' 485 : LESOBL(9)='V1Y ' 486 : LESOBL(10)='V1Z ' 487 : LESOBL(11)='V2X ' 488 : LESOBL(12)='V2Y ' 489 : LESOBL(13)='V2Z ' 490 : ELSE IF (IDIM.EQ.2) THEN 491 : IF (IFOUR.LE.0) THEN 492 : NBROBL=6 493 : SEGINI NOMID 494 : LESOBL(1)='PER1' 495 : LESOBL(2)='PER2' 496 : LESOBL(3)='PE12' 497 : LESOBL(4)='VISC' 498 : LESOBL(5)='V1X ' 499 : LESOBL(6)='V1Y ' 500 : ELSE IF (IFOUR.EQ.1) THEN 501 : NBROBL=7 502 : SEGINI NOMID 503 : LESOBL(1)='PER1' 504 : LESOBL(2)='PER2' 505 : LESOBL(3)='PE12' 506 : LESOBL(4)='PER3' 507 : LESOBL(5)='VISC' 508 : LESOBL(6)='V1X ' 509 : LESOBL(7)='V1Y ' 510 : ENDIF 511 : ENDIF 512 : * cas unidirectionnel 513 : ELSE IF (MATE.EQ.4) THEN 514 : IF (IDIM.EQ.3) THEN 515 : NBROBL=8 516 : SEGINI NOMID 517 : LESOBL(1)='PERM' 518 : LESOBL(2)='VISC' 519 : LESOBL(3)='V1X ' 520 : LESOBL(4)='V1Y ' 521 : LESOBL(5)='V1Z ' 522 : LESOBL(6)='V2X ' 523 : LESOBL(7)='V2Y ' 524 : LESOBL(8)='V2Z ' 525 : ELSE 526 : NBROBL=4 527 : SEGINI NOMID 528 : LESOBL(1)='PERM' 529 : LESOBL(2)='VISC' 530 : LESOBL(3)='V1X ' 531 : LESOBL(4)='V1Y ' 532 : ENDIF 533 : ENDIF 534 : * 535 : NMATR=NBROBL 536 : NMATF=NBRFAC 537 : NMATT = NMATR+NMATF 538 : MOMATR=NOMID 539 : * 540 : NBTYPE=1 541 : SEGINI NOTYPE 542 : MOTYPE=NOTYPE 543 : TYPE(1)='REAL*8' 544 : * 544 : CALL KOMCHA(IPCHE1,IPMAIL,CONM,MOMATR,MOTYPE,1,INFOS,3,IVAMAT) 545 : SEGSUP NOTYPE 546 : IF (IERR.NE.0) GOTO 9992 547 : IF (ISUP.EQ.1) THEN 547 : CALL VALCHE(IVAMAT,NMATR,IPMINT,IPPORE) 548 : ENDIF 549 : C 550 : C TRAITEMENT DES CHAMPS DE CARACTERISTIQUES * 551 : C 552 : NBROBL=0 553 : NBRFAC=0 554 : MOCARA=0 555 : IVECT=0 556 : * 557 : * EPAISSEUR DANS LE CAS DES CONTRAINTES PLANES 558 : * 559 : IF(IFOUR.EQ.-2.AND.((MELE.GE.79.AND.MELE.LE.83) 560 : & .OR.(MELE.GE.173.AND.MELE.LE.182)))THEN 561 : * 562 : NBROBL=0 563 : NBRFAC=1 564 : SEGINI NOMID 565 : MOCARA=NOMID 566 : LESFAC(1)='DIM3' 567 : * 568 : NBTYPE=1 569 : SEGINI NOTYPE 570 : TYPE(1)='REAL*8' 571 : ENDIF 572 : * 573 : NCARA=NBROBL 574 : NCARF=NBRFAC 575 : NCARR=NCARA+NCARF 576 : * 577 : IF (MOCARA.NE.0) THEN 578 : MOTYPE=NOTYPE 578 : CALL KOMCHA(IPCHE1,IPMAIL,CONM,MOCARA,MOTYPE,1,INFOS,3, 579 : $ IVACAR) 580 : SEGSUP NOTYPE 581 : IF (IERR.NE.0) GOTO 9992 582 : * 583 : IF (ISUP.EQ.1) THEN 583 : CALL VALCHE(IVACAR,NCARR,IPMINT,IPPORE) 584 : ENDIF 585 : ENDIF 586 : * 587 : *_____________________________________________________________________ 588 : * 589 : * NUMERO DES ETIQUETTES : 590 : * ETIQUETTES DE 1 A 98 POUR TRAITEMENT SPECIFIQUE A L ELEMENT 591 : * DANS LA ZONE SPECIFIQUE A CHAQUE ELEMENT COMMENCANT PAR : 592 : * 5 CONTINUE 593 : * ELEMENT 5 ETIQUETTES 1005 2005 3005 4005 ... 594 : * 44 CONTINUE 595 : * ELEMENT 44 ETIQUETTES 1044 2044 3044 4044 ... 596 : *_____________________________________________________________________ 597 : GOTO (99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99, 598 : 1 99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99, 599 : 2 99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99, 600 : 3 99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,79,79, 601 : 4 79,79,79,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99,99, 602 : 5 99,99,99,99,99,99,99,80,80,80),MELE 603 : * 604 : IF(MELE.GE.173.AND.MELE.LE.182) GO TO 173 605 : * 606 : 99 CONTINUE 607 : SEGSUP xMATRI 608 : MOTERR(1:4)=NOMTP(MELE) 609 : MOTERR(5:12)='PERMEABI' 609 : CALL ERREUR(86) 610 : GOTO 9990 611 : *_______________________________________________________________________ 612 : * 613 : * MILIEU POREUX 614 : *_______________________________________________________________________ 615 : * 616 : 79 CONTINUE 617 : * 618 : * Pour ces elements NBBB : Nb de noeuds 619 : * NBNO : Nb de fonctions de forme 620 : * 621 : DIM3=1.D0 622 : NBNO=IPORE 623 : NBBB=NBNN 624 : NSTN=1 625 : LRN=NBNO-NBBB 626 : NSTB=2 627 : IF(IFOUR.GT.0) NSTB=3 628 : * 629 : * CAS NON ISOTROPES 630 : * RECUPERATION DES FONCTIONS DE FORME ET LEURS DERIVEES 631 : * AU CENTRE DE L'ELEMENT POUR LE CALCUL DES AXES LOCAUX 632 : * 633 : IF(MATE.EQ.2.OR.MATE.EQ.3.OR.MATE.EQ.4) THEN 633 : CALL RESHPT(1,NBNO,IELE,MELE,0,IPT1,IRT1) 634 : MINTE2=IPT1 635 : SEGACT MINTE2 636 : SEGINI WRK4 637 : ENDIF 638 : * 639 : SEGINI WRK1,WRK2,WRK3 640 : I195=0 641 : I259=0 642 : I367=0 643 : DO 3079 IB=1,NBELEM 644 : * 645 : * On cherche les coordonnees des noeuds de l'element IB 646 : * 646 : CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) 647 : * 648 : * calcul des axes locaux dans les cas non isotropes 649 : * 650 : IF(MATE.EQ.2.OR.MATE.EQ.3.OR.MATE.EQ.4) THEN 651 : NBSH=MINTE2.SHPTOT(/2) 651 : CALL RLOCAL (XE,MINTE2.SHPTOT,NBSH,NBNN,TXR) 652 : if (nbsh.eq.-1) then 652 : call erreur(525) 653 : return 654 : endif 655 : ENDIF 656 : * 656 : CALL ZERO (REL,LRE,LRE) 657 : * 658 : * boucle sur les points de Gauss 659 : * 660 : ISDJC=0 661 : DO 4079 IGAU=1,NBPGAU 662 : 663 : * PRINT *, ' POINT DE GAUSS ',IGAU 664 : C 665 : C RECUPERATION DE L'EPAISSEUR 666 : C 667 : IF (IFOUR.EQ.-2)THEN 668 : MPTVAL=IVACAR 669 : IF (IVACAR.NE.0) THEN 670 : MELVAL=IVAL(1) 671 : IF (MELVAL.NE.0) THEN 672 : IGMN=MIN(IGAU,VELCHE(/1)) 673 : IBMN=MIN(IB,VELCHE(/2)) 674 : DIM3=VELCHE(IGMN,IBMN) 675 : ELSE 676 : DIM3=1.D0 677 : ENDIF 678 : ENDIF 679 : ENDIF 680 : * 680 : CALL BNPORE(IGAU,NBNO,NBBB,LRE,IFOUR,NSTB,NSTN,NHRM, 681 : & DIM3,XE,SHPTOT,SHPWRK,BGENE,XGENE,DJAC,2) 682 : IF(DJAC.LT.0.) ISDJC=ISDJC+1 683 : IF(DJAC.EQ.0.) I259 =IB 684 : DJAC=ABS(DJAC)*POIGAU(IGAU) 685 : 686 : * PRINT *,' MATRICE B LIGNE PAR LIGNE ' 687 : * DO 3367 IPZ = 1,NSTB 688 : ** PRINT *,' LIGNE ',IPZ 689 : * WRITE(6,3368) (BGENE(IPZ,JPZ), JPZ=1,LRE) 690 : *3368 FORMAT(8(1X,1PE10.3)/) 691 : *3367 CONTINUE 692 : 693 : EREF=1.D0 694 : MPTVAL=IVAMAT 695 : * 696 : * le cas isotrope 697 : * 698 : IF (MATE.EQ.1) THEN 699 : MELVAL=IVAL(1) 700 : IGMN=MIN(IGAU,VELCHE(/1)) 701 : IBMN=MIN(IB ,VELCHE(/2)) 702 : XK =VELCHE(IGMN,IBMN) 703 : * 704 : MELVAL=IVAL(2) 705 : IGMN=MIN(IGAU,VELCHE(/1)) 706 : IBMN=MIN(IB ,VELCHE(/2)) 707 : XMU =VELCHE(IGMN,IBMN) 708 : IF(XMU.EQ.0.D0) THEN 709 : I367=IB 710 : GO TO 4079 711 : ENDIF 712 : COMJAC=DJAC*EREF*EREF*XK/XMU 713 : DO 4279 I=1,LRN 714 : DO 4279 J=1,I 715 : DO 4279 K=1,NSTB 716 : REL(I,J)=REL(I,J)+COMJAC*BGENE(K,I)*BGENE(K,J) 717 : 4279 CONTINUE 718 : * 719 : * les cas non isotropes 720 : * 721 : ELSE IF(MATE.EQ.2.OR.MATE.EQ.3.OR.MATE.EQ.4) THEN 722 : * 723 : DO 4379 IM=1,NMATT 724 : IF (IVAL(IM).NE.0) THEN 725 : MELVAL=IVAL(IM) 726 : IBMN=MIN(IB ,VELCHE(/2)) 727 : IGMN=MIN(IGAU,VELCHE(/1)) 728 : VALMAT(IM)=VELCHE(IGMN,IBMN) 729 : ELSE 730 : VALMAT(IM)=0.D0 731 : ENDIF 732 : 4379 CONTINUE 733 : * 733 : CALL PERMAO(WRK4,IFOUR,MATE,EREF,KERRE) 734 : IF(KERRE.EQ.1) GO TO 99 735 : IF(KERRE.EQ.2) THEN 736 : I367=IB 737 : GO TO 4079 738 : ENDIF 739 : * 739 : CALL BDBST(BGENE,DJAC,PMAT,LRE,NSTB,REL) 740 : * 741 : * les cas non prévus 742 : * 743 : ELSE 744 : GO TO 99 745 : ENDIF 746 : * 747 : 4079 CONTINUE 748 : IF(ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) I195=IB 749 : * SEGINI XMATRI 750 : * IMATTT(IB)=XMATRI 751 : * 752 : * Remplissage de XMATRI 753 : * 753 : CALL REMPMT(REL,LRE,RE(1,1,ib)) 754 : * SEGDES XMATRI 755 : 3079 CONTINUE 756 : * 757 : SEGSUP WRK1,WRK2,WRK3 758 : IF(MATE.EQ.2.OR.MATE.EQ.3.OR.MATE.EQ.4) THEN 759 : SEGDES MINTE2 760 : SEGSUP WRK4 761 : ENDIF 762 : * 763 : IF(I195.NE.0) THEN 764 : INTERR(1)=I195 764 : CALL ERREUR(195) 765 : GOTO 9990 766 : ELSE IF(I259.NE.0) THEN 767 : INTERR(1)=I259 767 : CALL ERREUR(259) 768 : GOTO 9990 769 : ELSE IF(I367.NE.0) THEN 770 : INTERR(1)=I367 770 : CALL ERREUR(367) 771 : GOTO 9990 772 : ENDIF 773 : * 774 : SEGDES xMATRI 775 : * 776 : GOTO 610 777 : *_______________________________________________________________________ 778 : * 779 : * JOINTS EN FORMULATION MILIEUX POREUX 780 : *_______________________________________________________________________ 781 : * 782 : 80 CONTINUE 783 : * 784 : * Pour ces elements NBBB : Nb de noeuds 785 : * NBNO : Nb de fonctions de forme 786 : * 787 : NBNO=IPORE 788 : NBBB=NBNN 789 : NSTN=1 790 : LRN=(NBNO-NBBB)*3/2 791 : NSTB=2 792 : IF(IFOUR.EQ.2) NSTB=3 793 : * 794 : SEGINI WRK1,WRK2,WRK3,WRK5 795 : I195=0 796 : I259=0 797 : I367=0 798 : DO 3080 IB=1,NBELEM 799 : * 800 : * On cherche les coordonnees des noeuds de l'element IB 801 : * 801 : CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) 802 : * 803 : * calcul des axes locaux 804 : * 804 : CALL JOPLOC(XE,SHPTOT,NBBB,NBNO,IFOUR,XEL,BPSS) 805 : * 805 : CALL ZERO (REL,LRE,LRE) 806 : * 806 : CALL INTDEL(XNTH,XNTB,XNTT,LRN,MELE) 807 : * 808 : * boucle sur les points de Gauss 809 : * 810 : ISDJC=0 811 : DO 4080 IGAU=1,NBPGAU 812 : * 812 : CALL BNPORJ(IGAU,NBNO,NBBB,LRE,IFOUR,NSTB,NSTN,XE,XEL, 813 : . SHPTOT,SHPWRK,BPSS,BGENE,XGENE,DJAC,3) 814 : IF(DJAC.LT.0.) ISDJC=ISDJC+1 815 : IF(DJAC.EQ.0.) I259 =IB 816 : DJAC=ABS(DJAC)*POIGAU(IGAU) 817 : * 818 : EREF=1.D0 819 : MPTVAL=IVAMAT 820 : * 821 : * le cas isotrope (le seul) 822 : * 823 : MELVAL=IVAL(1) 824 : IGMN=MIN(IGAU,VELCHE(/1)) 825 : IBMN=MIN(IB ,VELCHE(/2)) 826 : XKT =VELCHE(IGMN,IBMN) 827 : * 828 : MELVAL=IVAL(2) 829 : IGMN=MIN(IGAU,VELCHE(/1)) 830 : IBMN=MIN(IB ,VELCHE(/2)) 831 : XKNH =VELCHE(IGMN,IBMN) 832 : * 833 : MELVAL=IVAL(3) 834 : IGMN=MIN(IGAU,VELCHE(/1)) 835 : IBMN=MIN(IB ,VELCHE(/2)) 836 : XKNB =VELCHE(IGMN,IBMN) 837 : * 838 : MELVAL=IVAL(4) 839 : IGMN=MIN(IGAU,VELCHE(/1)) 840 : IBMN=MIN(IB ,VELCHE(/2)) 841 : XMU =VELCHE(IGMN,IBMN) 842 : IF(XMU.EQ.0.D0) THEN 843 : I367=IB 844 : GO TO 4080 845 : ENDIF 846 : COMJAT=DJAC*EREF*EREF*XKT/XMU 847 : COMJNH=DJAC*EREF*EREF*XKNH/XMU 848 : COMJNB=DJAC*EREF*EREF*XKNB/XMU 849 : DO 4280 I=1,LRN 850 : DO 4280 J=1,I 851 : REL(I,J)=REL(I,J)+COMJAT*BGENE(1,I)*BGENE(1,J) 852 : . *XNTT(I)*XNTT(J) 853 : . +COMJNH*XGENE(1,I)*XGENE(1,J)*XNTH(I,J) 854 : . +COMJNB*XGENE(1,I)*XGENE(1,J)*XNTB(I,J) 855 : IF(IFOUR.EQ.2)THEN 856 : REL(I,J)=REL(I,J)+COMJAT*BGENE(2,I)*BGENE(2,J) 857 : . *XNTT(I)*XNTT(J) 858 : ENDIF 859 : 4280 CONTINUE 860 : * 861 : 4080 CONTINUE 862 : IF(ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) I195=IB 863 : * SEGINI XMATRI 864 : * IMATTT(IB)=XMATRI 865 : * 866 : * Remplissage de XMATRI 867 : * 867 : CALL REMPMT(REL,LRE,RE(1,1,ib)) 868 : * SEGDES XMATRI 869 : 3080 CONTINUE 870 : * 871 : SEGSUP WRK1,WRK2,WRK3,WRK5 872 : * 873 : IF(I195.NE.0) THEN 874 : INTERR(1)=I195 874 : CALL ERREUR(195) 875 : GOTO 9990 876 : ELSE IF(I259.NE.0) THEN 877 : INTERR(1)=I259 877 : CALL ERREUR(259) 878 : GOTO 9990 879 : ELSE IF(I367.NE.0) THEN 880 : INTERR(1)=I367 880 : CALL ERREUR(367) 881 : GOTO 9990 882 : ENDIF 883 : * 884 : SEGDES xMATRI 885 : * 886 : GOTO 610 887 : *_______________________________________________________________________ 888 : * 889 : * MILIEU POREUX - SUITE 890 : *_______________________________________________________________________ 891 : * 892 : 173 CONTINUE 893 : * 894 : * Pour ces elements NBBB : Nb de noeuds 895 : * NBNO : Nb de fonctions de forme 896 : * 897 : DIM3=1.D0 898 : NBNO=IPORE 899 : NBBB=NBNN 900 : NSTN=IDECAP 901 : LPP=NBNO-NBBB 902 : LRN=IDECAP*LPP 903 : NSTBE=2 904 : IF(IFOUR.GT.0) NSTBE=3 905 : NSTB=NSTBE*IDECAP 906 : 907 : * PRINT *,'NSTBE=',NSTBE 908 : * PRINT *,'NSTB=',NSTB 909 : * PRINT *,'IDECAP=',IDECAP 910 : * PRINT *,'LRE =',LRE 911 : 912 : * 913 : * CAS NON ISOTROPES 914 : * NON PREVU 915 : * 916 : IF(MATE.EQ.2.OR.MATE.EQ.3.OR.MATE.EQ.4) THEN 916 : CALL ERREUR(251) 917 : GO TO 9990 918 : ENDIF 919 : * 920 : NPKK = IDECAP 921 : SEGINI WRK1,WRK2,WRK3,WRK6 922 : I195=0 923 : I259=0 924 : I367=0 925 : DO 3173 IB=1,NBELEM 926 : 927 : * PRINT *,'ELEMENT ' , IB 928 : * 929 : * On cherche les coordonnees des noeuds de l'element IB 930 : * 930 : CALL DOXE(XCOOR,IDIM,NBNN,NUM,IB,XE) 931 : * 931 : CALL ZERO (REL,LRE,LRE) 932 : * 933 : * boucle sur les points de Gauss 934 : * 935 : ISDJC=0 936 : DO 4173 IGAU=1,NBPGAU 937 : 938 : * PRINT *, ' POINT DE GAUSS ',IGAU 939 : 940 : 941 : C 942 : C RECUPERATION DE L'EPAISSEUR 943 : C 944 : IF (IFOUR.EQ.-2)THEN 945 : MPTVAL=IVACAR 946 : IF (IVACAR.NE.0) THEN 947 : MELVAL=IVAL(1) 948 : IF (MELVAL.NE.0) THEN 949 : IGMN=MIN(IGAU,VELCHE(/1)) 950 : IBMN=MIN(IB,VELCHE(/2)) 951 : DIM3=VELCHE(IGMN,IBMN) 952 : ELSE 953 : DIM3=1.D0 954 : ENDIF 955 : ENDIF 956 : ENDIF 957 : * 957 : CALL BNQORE(IGAU,NBNO,NBBB,LRE,IFOUR,NSTB,NSTN,NHRM, 958 : & DIM3,XE,SHPTOT,SHPWRK,BGENE,XGENE,DJAC,IDECAP,2) 959 : IF(DJAC.LT.0.) ISDJC=ISDJC+1 960 : IF(DJAC.EQ.0.) I259 =IB 961 : DJAC=ABS(DJAC)*POIGAU(IGAU) 962 : 963 : * PRINT *,' MATRICE B LIGNE PAR LIGNE ' 964 : * DO 1367 IPZ = 1,NSTB 965 : * PRINT *,' LIGNE ',IPZ 966 : * WRITE(6,1368) (BGENE(IPZ,JPZ), JPZ=1,LRE) 967 : *1368 FORMAT(8(1X,1PE10.3)/) 968 : *1367 CONTINUE 969 : 970 : * 971 : EREF=1.D0 972 : MPTVAL=IVAMAT 973 : * 974 : * le cas isotrope 975 : * 976 : IF (MATE.EQ.1) THEN 977 : 978 : ICO=1 979 : IDECA2=IDECAP*IDECAP 980 : DO 1731 ICD = 1,IDECAP 981 : ICDA =(ICD -1) * NSTBE 982 : DO 1732 JCD = 1,IDECAP 983 : JCDA =(JCD -1) * NSTBE 984 : MELVAL=IVAL(ICO) 985 : IGMN=MIN(IGAU,VELCHE(/1)) 986 : IBMN=MIN(IB ,VELCHE(/2)) 987 : DO 1733 KCD = 1,NSTBE 988 : PKK(ICDA+KCD,JCDA+KCD) =VELCHE(IGMN,IBMN) 989 : 1733 CONTINUE 990 : ICO=ICO+1 991 : 1732 CONTINUE 992 : 1731 CONTINUE 993 : * 994 : 995 : * PRINT *,' MATRICE PKK' 996 : * IF (IDECAP.EQ.2) THEN 997 : * WRITE (6,1342) ((PKK(I,J),J=1,NSTB),I=1,NSTB) 998 : *1342 FORMAT(4(1X,1PE12.5)/) 999 : * 1000 : * ELSE IF (IDECAP.EQ.3) THEN 1001 : * WRITE (6,1343) ((PKK(I,J),J=1,NSTB),I=1,NSTB) 1002 : *1343 FORMAT(6(1X,1PE12.5)/) 1003 : * ENDIF 1004 : 1005 : COMJAC=DJAC*EREF*EREF 1005 : CALL BDBSTS(BGENE,COMJAC,PKK,LRE,NSTB,REL) 1006 : * 1007 : * les cas non prévus 1008 : * 1009 : ELSE 1010 : GO TO 99 1011 : ENDIF 1012 : * 1013 : 4173 CONTINUE 1014 : IF(ISDJC.NE.0.AND.ISDJC.NE.NBPGAU) I195=IB 1015 : * SEGINI XMATRI 1016 : * IMATTT(IB)=XMATRI 1017 : * 1018 : * Remplissage de XMATRI 1019 : * 1019 : CALL REMPMS(REL,LRE,RE(1,1,ib)) 1020 : * SEGDES XMATRI 1021 : 3173 CONTINUE 1022 : * 1023 : SEGSUP WRK1,WRK2,WRK3,WRK6 1024 : IF(MATE.EQ.2.OR.MATE.EQ.3.OR.MATE.EQ.4) THEN 1025 : SEGDES MINTE2 1026 : SEGSUP WRK4 1027 : ENDIF 1028 : * 1029 : IF(I195.NE.0) THEN 1030 : INTERR(1)=I195 1030 : CALL ERREUR(195) 1031 : GOTO 9990 1032 : ELSE IF(I259.NE.0) THEN 1033 : INTERR(1)=I259 1033 : CALL ERREUR(259) 1034 : GOTO 9990 1035 : ELSE IF(I367.NE.0) THEN 1036 : INTERR(1)=I367 1036 : CALL ERREUR(367) 1037 : GOTO 9990 1038 : ENDIF 1039 : * 1040 : SEGDES xMATRI 1041 : * 1042 : GOTO 610 1043 : * 1044 : * Desactivation des segment propres a la geometrie ISOUS 1045 : * 1046 : 610 CONTINUE 1047 : * 1048 : SEGDES MELEME 1049 : SEGDES MINTE 1050 : * 1051 : IF (ISUP.EQ.1) THEN 1051 : CALL DTMVAL(IVAMAT,3) 1052 : ELSE 1052 : CALL DTMVAL(IVAMAT,1) 1053 : ENDIF 1054 : * 1055 : NOMID=MOMATR 1056 : SEGSUP NOMID 1057 : IF (lsupdp) THEN 1058 : NOMID=MODEPL 1059 : SEGSUP,NOMID 1060 : ENDIF 1061 : IF (lsupfo) THEN 1062 : NOMID=MOFORC 1063 : SEGSUP,NOMID 1064 : ENDIF 1065 : * 1066 : COERIG(ISORI) = 1.D0 1067 : IRIGEL(1,ISORI)=IPMAIL 1068 : IRIGEL(2,ISORI)=0 1069 : IRIGEL(3,ISORI)=IPDES 1070 : IRIGEL(4,ISORI)=xMATRI 1071 : IRIGEL(5,ISORI)=NHRM 1072 : IRIGEL(6,ISORI)=0 1073 : IRIGEL(7,ISORI)=0 1074 : IRIGEL(8,ISORI)=0 1075 : 1076 : 1077 : * Fin de la boucle de PARTITIONNEMENT du segment XMATRI 1078 : 5000 CONTINUE 1079 : 1080 : SEGDES IMODEL 1081 : 1082 : 500 CONTINUE 1083 : IF(ISORI.NE.NRIGEL) GO TO 9999 1084 : SEGDES MRIGID 1085 : SEGSUP MMODEL 1086 : IRET = 1 1087 : IPRIGI = MRIGID 1088 : RETURN 1089 : * 1090 : * Erreur dans une sous zone desactivation et retour 1091 : * 1092 : 9990 CONTINUE 1093 : IF (ISUP.EQ.1) THEN 1093 : CALL DTMVAL(IVAMAT,3) 1094 : ELSE 1094 : CALL DTMVAL(IVAMAT,1) 1095 : ENDIF 1096 : * 1097 : 9992 CONTINUE 1098 : SEGSUP DESCR 1099 : SEGSUP xMATRI 1100 : SEGDES MELEME 1101 : SEGDES MINTE 1102 : * 1103 : NOMID=MOMATR 1104 : SEGSUP NOMID 1105 : 9999 CONTINUE 1106 : SEGDES IMODEL 1107 : SEGSUP MMODEL 1108 : SEGSUP MRIGID 1109 : IRET = 0 1110 : IPRIGI = 0 1111 : 1112 : 9991 CONTINUE 1113 : 1114 : MMODEL = IPMODL 1115 : DO isous = 1, NSOUS 1116 : IMODEL = KMODEL(isous) 1117 : SEGDES,IMODEL 1118 : ENDDO 1119 : SEGSUP,MMODEL 1120 : 1121 : RETURN 1122 : END 1123 : 1124 : 1125 : 1126 : 1127 : 1128 :
© Cast3M 2003 - All rights reserved.
Disclaimer