1 : C SIGMAX SOURCE FANDEUR 13/01/29 21:16:28 7683 2 : C RIGIX SOURCE PV 08/09/23 21:15:05 6164 3 : subroutine SIGMAX (MATE,IMAT,NBGMAT,NELMAT,NMATT,CMATE, 4 : & IVAMAT,IMODEL,IREPS2,IVADEP, 5 : & IVASTR,UZDPG,RYDPG,RXDPG,IIPDPG,IRETER) 6 : 7 : c 8 : C PROCEDURE UTILISEE DANS LE CAS D'ELEMENTS XFEM 9 : c POUR LE CALCUL DE la contrainte (élastique) 10 : C 11 : C 12 : C********************************************************* 13 : C PARTIE DECLARATIVE 14 : C********************************************************* 15 : 16 : C 17 : IMPLICIT REAL*8(A-H,O-Z) 18 : C 19 : PARAMETER (NDDLMAX=30,NBNIMAX=10) 20 : CTY PARAMETER (NDDLMAX=20,NBNIMAX=10) 21 : 22 : PARAMETER (NBENRMAX=5) 23 : DIMENSION MLRE(NBENRMAX+1) 24 : C 25 : -INC CCOPTIO 26 : -INC SMCOORD 27 : -INC SMMODEL 28 : -INC SMCHAML 29 : -INC SMINTE 30 : -INC SMELEME 31 : -INC SMLREEL 32 : C 33 : POINTEUR MCHEX1.MCHELM,MINT1.MINTE,MINT2.MINTE 34 : C 35 : C Segment (type LISTENTI) contenant les informations sur un element 36 : SEGMENT INFO 37 : INTEGER INFELL(JG) 38 : ENDSEGMENT 39 : c 40 : SEGMENT WRK1 41 : REAL*8 XE(3,NBBB) 42 : REAL*8 DDHOOK(LHOOK,LHOOK) 43 : REAL*8 XDDL(LRE),XSTRS(LHOOK) 44 : ENDSEGMENT 45 : c 46 : SEGMENT WRK2 47 : REAL*8 SHPWRK(6,NBNO),BGENE(LHOOK,LRE) 48 : c REAL*8 LV7WRK(NBENRMA2,2,6,NBNO) 49 : REAL*8 LV7WRK(NBENRMA2,2,6,NBBB) 50 : ENDSEGMENT 51 : C 52 : SEGMENT,MVELCH 53 : REAL*8 VALMAT(NV1) 54 : ENDSEGMENT 55 : c 56 : SEGMENT MPTVAL 57 : INTEGER IPOS(NS),NSOF(NS) 58 : INTEGER IVAL(NCOSOU) 59 : CHARACTER*16 TYVAL(NCOSOU) 60 : ENDSEGMENT 61 : C 62 : SEGMENT MRACC 63 : INTEGER TLREEL(NBENRMA2,NBI) 64 : ENDSEGMENT 65 : C 66 : SEGMENT MTRACE 67 : REAL*8 TRACE(3,NBPTEL) 68 : ENDSEGMENT 69 : c 70 : CHARACTER*8 CMATE 71 : DIMENSION UDPGE(3) 72 : LOGICAL BDPGE 73 : C 74 : c write (*,*) '#############################' 75 : c write (*,*) '##### DEBUT DE SIGMAX #####' 76 : c write (*,*) '#############################' 77 : C 78 : C********************************************************* 79 : c Introduction du point autour duquel se fait le mouvement 80 : c de la section en defo plane generalisee 81 : C********************************************************* 82 : C En 1D : pas de rotation 83 : BDPGE=.FALSE. 84 : NDPGE=0 85 : XDPGE=0.D0 86 : YDPGE=0.D0 87 : IF (IFOUR.EQ.-3) THEN 88 : BDPGE=.TRUE. 89 : NDPGE=3 90 : UDPGE(1)=UZDPG 91 : UDPGE(2)=RYDPG 92 : UDPGE(3)=RXDPG 93 : SEGACT,MCOORD 94 : IREF=(IIPDPG-1)*(IDIM+1) 95 : XDPGE=XCOOR(IREF+1) 96 : YDPGE=XCOOR(IREF+2) 97 : ELSE IF (IDIM.EQ.1) THEN 98 : IF (IFOUR.EQ. 7 .OR. IFOUR.EQ. 8 .OR. IFOUR.EQ. 9 .OR. 99 : . IFOUR.EQ.10 .OR. IFOUR.EQ.14) THEN 100 : BDPGE=.TRUE. 101 : NDPGE=1 102 : UDPGE(1)=UZDPG 103 : ELSE IF (IFOUR.EQ.11) THEN 104 : BDPGE=.TRUE. 105 : NDPGE=2 106 : UDPGE(1)=UZDPG 107 : UDPGE(2)=RXDPG 108 : ENDIF 109 : ENDIF 110 : C 111 : NV1=NMATT 112 : SEGINI,MVELCH 113 : C 114 : C 115 : C********************************************************* 116 : c 117 : C RECUPERATION + ACTIVATIONS + VALEURS UTILES 118 : c 119 : C********************************************************* 120 : C SEGACT MMODEL,IMODEL deja activé dans epsi1 121 : c 122 : C++++ Activation au cas ou ++++++++++++++++++++++++++++++ 123 : SEGACT MCOORD 124 : 125 : C++++ Recup + Activation de la geometrie ++++++++++++++++ 126 : MELEME= IMAMOD 127 : c SEGACT MELEME deja activé dans epsi1 128 : 129 : 130 : C++++ RECUP DES INFOS EF ++++++++++++++++++++++++++++++++ 131 : c + OBTENUES PAR ELQUOI DANS RIGI1 PENDANT PHASE 1 132 : C segment INFO deja actif dans RIGI1 133 : c + rigi1 n appelle pas elquoi, c'est modeli qui l'a fait 134 : c mais du coup, on na pas de segment minte 135 : c (car depend de si pt de g pour rigi, pour sigma....) 136 : c c'est + simple de rappeler elquoi ici 137 : MELE = NEFMOD 137 : call elquoi(MELE,0,3,IPINF,IMODEL) 138 : INFO = IPINF 139 : c MELE = INFELL(1) 140 : c NBPGA2= INFELL(2) 141 : c NBPGAU= INFELL(3) 142 : c NBPGAU= INFELL(4) 143 : c ICARA = INFELL(5) 144 : NGAU1 = INFELL(6) 145 : c LW = INFELL(7) 146 : LRE = INFELL(9) 147 : LHOOK = INFELL(10) 148 : MINT1 = INFELL(11) 149 : segact,MINT1 150 : MINT2 = INFELL(12) 151 : if(MINT2.ne.0) segact,MINT2 152 : MFR = INFELL(13) 153 : IELE = INFELL(14) 154 : NDDL = INFELL(15) 155 : NSTRS = INFELL(16) 156 : c write(6,*)'-> EPSIX infell',(infell(iou),iou=1,16) 157 : 158 : c + AUTRES INFOS 159 : C nbre de noeuds par element 160 : NBNN1 = NUM(/1) 161 : C nbre d elements 162 : NBEL1 = NUM(/2) 163 : 164 : c REM: pour se passer du dimensionnement du nbre d'enrichissement dans 165 : c elquoi et le realiser localement , on pourrait ecrire: 166 : c LRE = NDDLMAX*NBNN1 167 : c NDDL= NDDLMAX 168 : 169 : C sous decoupage et points de Gauss de l element geometrique de base 170 : CTY if(MELE.eq.263) then 171 : if(MELE.EQ.263.OR.MELE.EQ.264) then 172 : NGAU2 = MINT2.POIGAU(/1) 173 : endif 174 : c write(*,*) 'dim de MINT2=6,',(MINT2.SHPTOT(/2)),(MINT2.SHPTOT(/3)) 175 : c write(*,*) 'MINT2',(MINT2.QSIGAU(iou),iou=1,NGAU) 176 : c segdes,MINT2 177 : 178 : 179 : c nbre maxi de fonction de forme par noeud (fonction std comprise) 180 : c NBNI = NDDL/IDIM inutile! 181 : 182 : 183 : C++++ Recup des infos d enrichissement +++++++++++++++++++ 184 : c recup du MCHEX1 d'enrichissement 185 : NOBMO1 = IVAMOD(/1) 186 : if(NOBMO1.ne.0) then 187 : do iobmo1=1,NOBMO1 188 : if((TYMODE(iobmo1)).eq.'MCHAML') then 189 : MCHEX1 = IVAMOD(iobmo1) 190 : segact,MCHEX1 191 : if((MCHEX1.TITCHE).eq.'ENRICHIS') then 192 : MCHAM1 = MCHEX1.ICHAML(1) 193 : segact,MCHAM1 194 : segdes,MCHEX1 195 : goto 1000 196 : endif 197 : segdes,MCHEX1 198 : endif 199 : enddo 200 : write(*,*) 'Le modele est vide (absence d enrichissement)' 201 : * return 202 : else 203 : write(*,*) 'Il n y a pas de MCHEML enrichissement dans le Modele' 204 : * return 205 : endif 206 : 207 : 1000 continue 208 : c niveau d enrichissement(s) du modele (ddls std u exclus) 209 : c NBENR1= 0 si std, 1 si H seul, 2 si H+F1, 3 si H+F1+F2, etc... 210 : if(NOBMO1.ne.0) then 211 : NBENR1= MCHAM1.IELVAL(/1) 212 : else 213 : NBENR1 = 0 214 : endif 215 : c write(*,*) 'niveau d enrichissement(s) du modele',NBENR1 216 : c 217 : C********************************************************* 218 : C INITIALISATIONS... 219 : C********************************************************* 220 : IRETER = 0 221 : c 222 : c preparation des tables avec: 223 : 224 : if(NBENR1.ne.0) then 225 : do ienr=1,NBENR1 226 : c -le nombre d'inconnues de chaque sous-zone 227 : c determinee depuis le nombre de fonction de forme 228 : c ienr= 1: U+H(1+1=2), 2: U+H+F1(2+4=6), 3: U+H+F1+F2(6+4=10) 229 : nbniJ = 2 + ((ienr-1)*4) 230 : MLRE(1+ienr) = IDIM*NBNN1*nbniJ 231 : enddo 232 : endif 233 : C Tables + longues car 1er indice correspond au fontion de forme std 234 : MLRE(1) = IDIM*NBNN1*1 235 : 236 : if(NBENR1.lt.(NBENRMAX+1)) then 237 : do ienr=(NBENR1+1),(NBENRMAX) 238 : MLRE(1+ienr) = 0 239 : enddo 240 : endif 241 : c 242 : c ...DU SEGMENT WRK1 243 : NBENRMA2 = NBENRMAX 244 : NBBB = NBNN1 245 : SEGINI,WRK1 246 : 247 : c ...DU SEGMENT WRK2 248 : c NBNO = NBNI 249 : NBNO = LRE/IDIM 250 : SEGINI,WRK2 251 : C 252 : c ...DU SEGMENT MRACC 253 : NBENRMA2 = NBENRMAX 254 : NBI = NBNN1 255 : segini,MRACC 256 : C 257 : C du nombre d erreur sur le noms de composantes 258 : NBERR1=0 259 : 260 : C********************************************************* 261 : C 262 : C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> BOUCLE SUR LES ELEMENTS 263 : C 264 : DO 2000 J=1,NBEL1 265 : c write(*,*) '===================================' 266 : c write(*,*) '========= EF',J,' /NBEL1 =========' 267 : C 268 : C 269 : C********************************************************* 270 : C POUR CHAQUE ELEMENT, ... 271 : C********************************************************* 272 : C 273 : C++++ ON RECUPERE LES COORDONNEES DES NOEUDS DE L ELEMENT 273 : CALL DOXE(XCOOR,IDIM,NBNN1,NUM,J,XE) 274 : C 275 : c 276 : C++++ NBENRJ = niveau d enrichissement de l element ++++ 277 : C =0 si EF std =1 si U+H =2 si U+H+F1 =3 si U+H+F1+F2 278 : NBENRJ=0 279 : if(NBENR1.ne.0) then 280 : do IENR=1,NBENR1 281 : MELVA1 = MCHAM1.IELVAL(IENR) 282 : segact,MELVA1 283 : do I=1,NBNN1 284 : mlree1 = MELVA1.IELCHE(I,J) 285 : C on en profite pour remplir MRACC table de raccourcis pour cet element 286 : TLREEL(IENR,I) = mlree1 287 : if(mlree1.ne.0) then 288 : NBENRJ = max(NBENRJ,IENR) 289 : C et on active la listreel 290 : segact,mlree1 291 : endif 292 : enddo 293 : segdes,MELVA1 294 : enddo 295 : endif 296 : if(NBENRMA2.gt.NBENR1) then 297 : do IENR2=(NBENR1+1),NBENRMA2 298 : do I=1,NBNN1 299 : TLREEL(IENR2,I) = 0 300 : enddo 301 : enddo 302 : endif 303 : C 304 : c 305 : c++++ quelques indices pour dimensionner au plus juste 306 : c nbre total de ddl de l'élément considéré 307 : NLIGRD = MLRE(1+NBENRJ) 308 : NLIGRP = MLRE(1+NBENRJ) 309 : c nbre fonction de forme=((Ni_std+Ni_enrichi)*nbnoeud)=Ninconnues/idim 310 : NBNI = (MLRE(1+NBENRJ)) / IDIM 311 : 312 : c write(*,*) 'EF',J,' NBENRJ=',NBENRJ, 313 : c &'donc',NLIGRD,' ddls et ',NBNI,' fonctions de forme' 314 : c 315 : C++++ ON CALCULE XDDL ++++ 316 : MPTVAL = IVADEP 317 : NCOSOU = IVAL(/1) 318 : 319 : INITOT = 0 320 : C-->> BOUCLE SUR LES niveaux d'ENRICHISSEMENTS (0:U 1:A 2:BCDE1 3:BCDE2) 321 : DO IENR=0,NBENRJ 322 : *nbre de fonction(s) de ce niveau d'enrichissement (=1 si std ou H, =4 pour F1,2,...) 323 : if(IENR .le. 1) then 324 : NBNIENR = 1 325 : else 326 : NBNIENR = 4 327 : endif 328 : C---->> BOUCLE SUR LES fonctions de forme de ce type d'enrichissement 329 : DO INI=1,NBNIENR 330 : INITOT = INITOT + 1 331 : C------>> BOUCLE SUR LA DIMENSION 332 : DO 2220 KDIM=1,IDIM 333 : ICOMP = (INITOT-1)*IDIM + KDIM 334 : 335 : c --cas ou on n'a pas trouvé assez de composantes-- 336 : if(ICOMP.GT.NCOSOU) GOTO 2221 337 : 338 : MELVAL = IVAL(ICOMP) 339 : c --cas ou on a pas trouvé cette composante dans cette zone du 340 : c chpoint solution devenu mchaml -- 341 : if(MELVAL.eq.0)then 342 : c Avait on besoin de cette composante? 343 : c oui, si c'est une composante obligatoire 344 : if(IENR.eq.0) goto 2991 345 : c oui, si l'un des noeuds est enrichi 346 : do I=1,NBNN1 347 : if(TLREEL(IENR,I).ne.0) goto 2991 348 : enddo 349 : c non, si c'est facultatif et qu'on n'est pas enrichi -> on saute 350 : goto 2220 351 : c ->AVERTISSEMENT puis on saute 352 : 2991 NBERR1=NBERR1+1 353 : if(IIMPI.lt.1) goto 2220 354 : c write(IOIMP,*) 'PB OPERATEUR SIGM :' 355 : write(IOIMP,991) ICOMP,IENR,INI,KDIM 356 : 991 format(2X,'ABSENCE DANS LE CHPOINT DEPLACEMENT DE LA ',I3, 357 : $ 'ieme COMPOSANTE (enrichissement',I3, 358 : $ ', fonction',I3,', direction ',I3, 359 : $ ') NECESSAIRE POUR L UN DES NOEUDS SUIVANTS :') 360 : write(IOIMP,*)' noeuds :',(NUM(iou,J),iou=1,NBNN1) 361 : goto 2220 362 : endif 363 : 364 : C---------->> BOUCLE SUR LES NOEUDS 365 : DO I=1,NBNN1 366 : IQ = ((INITOT-1)*NBNN1*IDIM) + ((I-1)*IDIM) + KDIM 367 : XDDL(IQ) = VELCHE(I,J) 368 : ENDDO 369 : 370 : 2220 CONTINUE 371 : ENDDO 372 : ENDDO 373 : 374 : c --cas normal (toutes les composantes souhaitees etaient presentes)-- 375 : GOTO 2223 376 : 377 : c --cas ou on n'a pas trouvé assez de composantes-- 378 : 2221 CONTINUE 379 : if (IIMPI.ge.1) then 380 : WRITE(IOIMP,2222) J,NCOSOU,ICOMP 381 : 2222 FORMAT(2X,'ATTENTION : ELEMENT ',I6, 382 : & ' LE CHAMP DE DEPLACEMENT CONTIENT ',I3,' COMPOSANTES', 383 : & ' PAR NOEUD AU LIEU DE ',I3,' ATTENDUES') 384 : endif 385 : NDDL=NCOSOU*NBNN1 386 : NBENRJ=IENR 387 : 388 : 2223 CONTINUE 389 : c --cas ou on a une ou des erreurs-- 390 : IF (NBERR1.gt.0.and.J.eq.NBEL1) THEN 391 : write(IOIMP,*) 'OPERATEUR GRAD : ABSENCE DANS LE CHPOINT ', 392 : & 'DEPLACEMENT DE CERTAINES INCONNUES ATTENDUES PAR LE MODELE' 393 : ENDIF 394 : 395 : c 396 : c rem: il serait probablement interessant au niveau du tems cpu 397 : c d'utiliser moins de pts de Gauss lorsque l element est élastique 398 : c On pourrait par ex. utiliser MINT2 = infell(12) qui contient 399 : c le segment d'integration de l'EF std (QUA4 par ex.) 400 : * if((NBENRJ.eq.0).and.(MINT2.ne.0)) then 401 : * MINTE = MINT2 402 : * NBPGAU= NGAU2 403 : * else 404 : MINTE = MINT1 405 : NBPGAU= NGAU1 406 : * endif 407 : c 408 : c pour les def quadratiques 409 : ISDJC=0 410 : NBPTEL=NBPGAU 411 : IF (IREPS2.EQ.1) SEGINI MTRACE 412 : c 413 : C 414 : C********************************************************* 415 : C 416 : C>>>>>>>>>> BOUCLE SUR LES POINTS DE GAUSS 417 : C 418 : DO 2500 KGAU=1,NBPGAU 419 : C 420 : C********************************************************* 421 : C Initialisation à 0 422 : C********************************************************* 423 : 424 : c ZERO ne serait-il pas facultatif? 424 : CALL ZERO(SHPWRK,6,NBNO) 425 : XGAU = 0.D0 426 : YGAU = 0.D0 427 : ZGAU = 0.D0 428 : C 429 : i6zz = 3 430 : IF (IDIM.EQ.3) i6zz = 4 431 : c do ienr7=1,NBENRMAX 432 : do ienr7=1,NBENRJ 433 : do inod7=1,NBNN1 434 : c do i6=1,6 435 : do i6=1,i6zz 436 : LV7WRK(ienr7,1,i6,inod7) = 0 437 : LV7WRK(ienr7,2,i6,inod7) = 0 438 : enddo 439 : enddo 440 : enddo 441 : 442 : 443 : c********************************************************* 444 : c Calcul des fonction de forme std dans repere local 445 : c********************************************************* 446 : 447 : ccccc BOUCLE SUR LES NOEUDS ccccccccccccccccccccccccccccc 448 : c (et donc sur les Ni std) 449 : DO 2510 I=1,NBNN1 450 : 451 : C++++ Calcul des Ni std 452 : c (rappel: 1:Ni 2:Ni,qsi 3:Ni,eta avec i=1,4) 453 : SHPWRK(1,I) = SHPTOT(1,I,KGAU) 454 : SHPWRK(2,I) = SHPTOT(2,I,KGAU) 455 : SHPWRK(3,I) = SHPTOT(3,I,KGAU) 456 : IF (IDIM.EQ.3) SHPWRK(4,I) = SHPTOT(4,I,KGAU) 457 : 458 : C++++ Calcul des coordonnees dans le repere global 459 : XGAU = XGAU + (SHPWRK(1,I)*XE(1,I)) 460 : YGAU = YGAU + (SHPWRK(1,I)*XE(2,I)) 461 : IF (IDIM.EQ.3) ZGAU = ZGAU + (SHPWRK(1,I)*XE(3,I)) 462 : 463 : 2510 CONTINUE 464 : ccccc fin de BOUCLE SUR LES NOEUDS ccccccccccccccccccccccc 465 : 466 : 467 : 468 : c********************************************************* 469 : c Passage des fonctions de forme std dans repere global 470 : c********************************************************* 471 : 472 : C++++ CALCUL DES Ni,x Ni,y (i=1,4) + CALCUL DE det(J) 472 : CALL JACOBI(XE,SHPWRK,IDIM,NBNN1,DJAC) 473 : c if(J.eq.1.and.KGAU.eq.1) 474 : c &write(*,*) 'Ni(i=1,4)=',(SHPWRK(1,iou),iou=1,NBNN1) 475 : 476 : c********************************************************* 477 : c Si on est pas du tout enrichi on peut sauter une grosse partie 478 : c********************************************************* 479 : if(NBENRJ.eq.0) goto 2999 480 : 481 : c********************************************************* 482 : c Calcul des level set + leurs derivees dans repere global 483 : c********************************************************* 484 : 485 : ccccc BOUCLE SUR LES enrichissements ccccccccccccccccccc 486 : do 2520 ienr=1,NBENRJ 487 : 488 : c MELVA1=MCHAM1.IELVAL(IENR) 489 : c segact,MELVA1 490 : 491 : ccccc BOUCLE SUR LES NOEUDS ccccccccccccccccccccccccccc 492 : DO 2521 I=1,NBNN1 493 : 494 : C++++ Le I eme noeud est-il ienr-enrichi? 495 : mlree1= TLREEL(IENR,I) 496 : if(mlree1.eq.0) goto 2521 497 : 498 : 499 : c Calcul du repere local de fissure(=PSI,PHI) 500 : c (rappel: 1,1:psi 1,2:phi 2,1 psi,x ...etc...) 501 : do 2522 inode=1,NBNN1 502 : c pour le H-enrichissement, on n a pas gardé PSI (inutile) 503 : if(ienr.ne.1) then 504 : c valeur de PSI au inode^ieme noeud 505 : xpsi1 = mlree1.PROG(inode) 506 : c qu on multiplie par la valeur de Ni^std au pt de G considéré 507 : LV7WRK(ienr,1,1,I)= LV7WRK(ienr,1,1,I) 508 : & + (SHPWRK(1,inode)*xpsi1) 509 : LV7WRK(ienr,1,2,I)= LV7WRK(ienr,1,2,I) 510 : & + (SHPWRK(2,inode)*xpsi1) 511 : LV7WRK(ienr,1,3,I)= LV7WRK(ienr,1,3,I) 512 : & + (SHPWRK(3,inode)*xpsi1) 513 : IF (IDIM.EQ.3) LV7WRK(ienr,1,4,I)= LV7WRK(ienr,1,4,I) 514 : & + (SHPWRK(4,inode)*xpsi1) 515 : c valeur de PHI au inode^ieme noeud 516 : xphi1 = mlree1.PROG(NBNN1+inode) 517 : else 518 : xphi1 = mlree1.PROG(inode) 519 : endif 520 : LV7WRK(ienr,2,1,I)= LV7WRK(ienr,2,1,I) 521 : & + (SHPWRK(1,inode)*xphi1) 522 : LV7WRK(ienr,2,2,I)= LV7WRK(ienr,2,2,I) 523 : & + (SHPWRK(2,inode)*xphi1) 524 : LV7WRK(ienr,2,3,I)= LV7WRK(ienr,2,3,I) 525 : & + (SHPWRK(3,inode)*xphi1) 526 : IF (IDIM.EQ.3) LV7WRK(ienr,2,4,I)= LV7WRK(ienr,2,4,I) 527 : & + (SHPWRK(4,inode)*xphi1) 528 : 2522 continue 529 : 530 : 2521 continue 531 : ccccc fin de BOUCLE SUR LES NOEUDS ccccccccccccccccccccccc 532 : 533 : 534 : 2520 CONTINUE 535 : ccccc fin de BOUCLE SUR LES enrichissements cccccccccccccccc 536 : 537 : c on a construit 538 : C LV7WRK(ienr, PSI/PHI, valeur/deriveeparqsi/pareta, iNOEUD) 539 : 540 : c********************************************************* 541 : c Ajout des fonctions de forme d enrichissement 542 : c + leurs derivees dans repere global 543 : c********************************************************* 543 : CTY CALL SHAPX(XGAU,YGAU,MELE,L7WRK,NBENRMAX,NBENRJ, 543 : CALL SHAPX(XGAU,YGAU,ZGAU,MELE,LV7WRK,NBENRMAX,NBENRJ, 544 : & TLREEL,J,SHPWRK,IRET) 545 : 546 : c retour a la partie commune aux EF enrichi et non enrichi 547 : 2999 continue 548 : 549 : C********************************************************* 550 : C CALCUL DE B' 551 : C********************************************************* 552 : c ZERO ne serait-il pas facultatif? 552 : c call ZERO(BGENE,LHOOK,NLIGRP) 553 : KB=1 554 : C boucle sur tous les Ni 555 : DO 3001 II=1,NBNI 556 : 557 : BGENE(1,KB) = SHPWRK(2,II) 558 : BGENE(2,KB+1) = SHPWRK(3,II) 559 : BGENE(4,KB) = SHPWRK(3,II) 560 : BGENE(4,KB+1) = SHPWRK(2,II) 561 : 562 : IF(IDIM.EQ.3) THEN 563 : BGENE(3,KB+2)=SHPWRK(4,II) 564 : BGENE(5,KB)=SHPWRK(4,II) 565 : BGENE(5,KB+2)=SHPWRK(2,II) 566 : BGENE(6,KB+1)=SHPWRK(4,II) 567 : BGENE(6,KB+2)=SHPWRK(3,II) 568 : ENDIF 569 : 570 : KB = KB + IDIM 571 : 572 : 3001 CONTINUE 573 : C 574 : c if(J.eq.5.and.KGAU.eq.1) then 575 : c if(KGAU.eq.1) then 576 : c write(*,*) 'BGENE(1,..)=',(BGENE(1,iou),iou=1,2*NBNI) 577 : c write(*,*) 'BGENE(2,..)=',(BGENE(2,iou),iou=1,2*NBNI) 578 : c write(*,*) 'BGENE(4,..)=',(BGENE(3,iou),iou=1,2*NBNI) 579 : c endif 580 : 581 : C********************************************************* 582 : C CALCUL DE D 583 : C********************************************************* 584 : 585 : c on cherche les matrices de Hooke 586 : MPTVAL=IVAMAT 587 : IF (IMAT.EQ.2) THEN 588 : MELVAL=IVAL(1) 589 : IBMN=MIN(J ,IELCHE(/2)) 590 : IGMN=MIN(KGAU,IELCHE(/1)) 591 : MLREEL=IELCHE(IGMN,IBMN) 592 : SEGACT MLREEL 593 : IF (KGAU.LE.NBGMAT.AND.(J.LE.NELMAT.OR.NBGMAT.GT.1)) 593 : 1 CALL DOHOOO(PROG,LHOOK,DDHOOK) 594 : SEGDES MLREEL 595 : 596 : ELSE IF (IMAT.EQ.1) THEN 597 : DO 9004 IM=1,NMATT 598 : IF (IVAL(IM).NE.0) THEN 599 : MELVAL=IVAL(IM) 600 : IBMN=MIN(J ,VELCHE(/2)) 601 : IGMN=MIN(KGAU,VELCHE(/1)) 602 : VALMAT(IM)=VELCHE(IGMN,IBMN) 603 : ELSE 604 : VALMAT(IM)=0.D0 605 : ENDIF 606 : 9004 CONTINUE 607 : c 608 : IF (MATE.EQ.2.OR.MATE.EQ.3.OR.MATE.EQ.4) THEN 609 : * IF (IGAU.LE.NBGMAT) 609 : * 1 CALL DOHMAO(VALMAT,CMATE,IFOUR,IDIM,TXR,XLOC,XGLOB,D1HO, 610 : * 2 ROTH,DDHOOK,LHOOK,1,IRTD) 611 : write(*,*) 'cas non prevu a ce jour pour les EF xfem' 612 : return 613 : ELSE 614 : IF (KGAU.LE.NBGMAT.AND.(J.LE.NELMAT.OR.NBGMAT.GT.1)) 614 : 1 CALL DOHMAS(VALMAT,CMATE,IFOUR,LHOOK,1,DDHOOK,IRTD) 615 : ENDIF 616 : ENDIF 617 : c 618 : c 619 : C********************************************************* 620 : C CALCUL DE sigma = D*B*u 621 : C********************************************************* 622 : c 622 : CALL DBST(BGENE,DDHOOK,XDDL,(NBNI*IDIM),LHOOK,XSTRS) 623 : c 624 : c calcul des eps 2 (termes quadratiques de la deformation) 625 : c 626 : IF (IREPS2.EQ.1) 626 : 1 CALL DBST2(SHPWRK,DDHOOK,XDDL,XE,NBNI,IFOUR,LHOOK,XSTRS, 627 : 2 TRACE,KGAU,XDPGE,YDPGE,UDPGE,NIFOUR) 628 : c 629 : c remplissage du segment contenant les contraintes 630 : c 631 : MPTVAL=IVASTR 632 : DO 7004 ICOMP=1,LHOOK 633 : MELVAL=IVAL(ICOMP) 634 : IBMN=MIN(J,VELCHE(/2)) 635 : VELCHE(KGAU,IBMN)=XSTRS(ICOMP) 636 : 7004 CONTINUE 637 : c 638 : * if(J.eq.5 .and. KGAU.eq.1) then 639 : c if(KGAU.eq.1) then 640 : c write(*,*) J,KGAU,'sigma(..)=',(XSTRS(iou),iou=1,LHOOK) 641 : c endif 642 : c 643 : 2500 CONTINUE 644 : C FIN DE BOUCLE SUR LES POINTS DE GAUSS <<<<<<<<<<<<<< 645 : C 646 : c quelques suppressions 647 : c (ici element non-incompressible=> pas besoin de MTRACE (cf epsi2) 648 : IF (IREPS2.EQ.1) THEN 649 : SEGSUP MTRACE 650 : ENDIF 651 : c 652 : c quelques desactivations 653 : do IENR=1,NBENR1 654 : do I=1,NBNN1 655 : mlree1 = TLREEL(IENR,I) 656 : if(mlree1.ne.0) segdes,mlree1 657 : enddo 658 : enddo 659 : c 660 : 2000 CONTINUE 661 : C FIN DE BOUCLE SUR LES ELEMENTS <<<<<<<<<<<<<<<<<<<<< 662 : c 663 : C********************************************************* 664 : C SUPPRESSION ET DESACTIVATION DE SEGMENTS 665 : C********************************************************* 666 : C 667 : SEGSUP,WRK1,WRK2,MVELCH 668 : 669 : SEGSUP,MRACC 670 : 671 : segdes,MINT1 672 : if(MINT2.ne.0) segdes,MINT2 673 : C 674 : RETURN 675 : END 676 : 677 :
© Cast3M 2003 - All rights reserved.
Disclaimer