1 : C MASSX SOURCE BP208322 13/10/02 21:15:20 7830 2 : C RIGIX SOURCE BP208322 08/12/10 21:15:26 6216 3 : subroutine MASSX (ivamat,ivacar,NMATT,CMATE, 4 : & IMODEL,LOCIRI) 5 : c 6 : C PROCEDURE UTILISEE DANS LE CAS D'ELEMENTS XFEM 7 : c POUR LE CALCUL DE RIGIDITE ELEMENTAIRES 8 : C 9 : C 10 : C********************************************************* 11 : C PARTIE DECLARATIVE 12 : C********************************************************* 13 : 14 : C 15 : IMPLICIT REAL*8(A-H,O-Z) 16 : C 17 : CHARACTER*8 CMATE 18 : CYT PARAMETER (NDDLMAX=20,NBNIMAX=10) 19 : PARAMETER (NDDLMAX=30,NBNIMAX=10) 20 : CHARACTER*4 MOTINC(NDDLMAX),MOTDUA(NDDLMAX) 21 : CYT DATA MOTINC/'UX ','UY ','AX ','AY ', 22 : CYT >'B1X ','B1Y ','C1X ','C1Y ','D1X ','D1Y ','E1X ','E1Y ', 23 : CYT >'B2X ','B2Y ','C2X ','C2Y ','D2X ','D2Y ','E2X ','E2Y '/ 24 : CYT DATA MOTDUA/ 'FX ','FY ','FAX ','FAY ', 25 : CYT >'FB1X','FB1Y','FC1X','FC1Y','FD1X','FD1Y','FE1X','FE1Y', 26 : CYT >'FB2X','FB2Y','FC2X','FC2Y','FD2X','FD2Y','FE2X','FE2Y'/ 27 : 28 : DATA MOTINC/'UX ','UY ','UZ ','AX ','AY ','AZ ', 29 : >'B1X ','B1Y ','B1Z ','C1X ','C1Y ','C1Z ','D1X ','D1Y ', 30 : >'D1Z ','E1X ','E1Y ','E1Z ','B2X ','B2Y ','B2Z ','C2X ', 31 : >'C2Y ','C2Z ','D2X ','D2Y ','D2Z ','E2X ','E2Y ','E2Z '/ 32 : DATA MOTDUA/'FX ','FY ', 'FZ ','FAX ','FAY ','FAZ ', 33 : >'FB1X','FB1Y','FB1Z','FC1X','FC1Y','FC1Z','FD1X','FD1Y', 34 : >'FD1Z','FE1X','FE1Y','FE1Z','FB2X','FB2Y','FB2Z','FC2X', 35 : >'FC2Y','FC2Z','FD2X','FD2Y','FD2Z','FE2X','FE2Y','FE2Z'/ 36 : C 37 : PARAMETER (NBENRMAX=5) 38 : C NBENRMAX deja defini dans rigixr.eso 39 : INTEGER LOCIRI(10,(1+NBENRMAX)) 40 : DIMENSION MLRE(NBENRMAX+1) 41 : C 42 : -INC CCOPTIO 43 : -INC SMCOORD 44 : -INC SMMODEL 45 : -INC SMCHAML 46 : -INC SMINTE 47 : -INC SMELEME 48 : -INC SMRIGID 49 : -INC SMLREEL 50 : C 51 : POINTEUR MCHEX1.MCHELM,MINT1.MINTE,MINT2.MINTE 52 : C 53 : C Segment (type LISTENTI) contenant les informations sur un element 54 : SEGMENT INFO 55 : INTEGER INFELL(JG) 56 : ENDSEGMENT 57 : 58 : SEGMENT WRK1 59 : REAL*8 XE(3,NBBB) 60 : REAL*8 REL(LRE,LRE),RINT(LRE,LRE) 61 : ENDSEGMENT 62 : C 63 : SEGMENT WRK2 64 : REAL*8 SHPWRK(6,NBNO),BGENE(LHOOK,LRE) 65 : c REAL*8 LV7WRK(NBENRMA2,2,6,NBNO) 66 : REAL*8 LV7WRK(NBENRMA2,2,6,NBBB) 67 : ENDSEGMENT 68 : C 69 : SEGMENT MVELCH 70 : REAL*8 VALMAT(NV1) 71 : ENDSEGMENT 72 : C 73 : SEGMENT MPTVAL 74 : INTEGER IPOS(NS),NSOF(NS) 75 : INTEGER IVAL(NCOSOU) 76 : CHARACTER*16 TYVAL(NCOSOU) 77 : ENDSEGMENT 78 : C 79 : SEGMENT MRACC 80 : INTEGER TLREEL(NBENRMA2,NBI) 81 : ENDSEGMENT 82 : C 83 : c write (ioimp,*) '############################' 84 : c write (ioimp,*) '# DEBUT DE MASSX #' 85 : c write (ioimp,*) '############################' 86 : C 87 : C********************************************************* 88 : c 89 : C RECUPERATION + ACTIVATIONS + VALEURS UTILES 90 : c 91 : C********************************************************* 92 : C sont deja actifs dans rigi1.eso : 93 : C SEGACT MMODEL,IMODEL,MELVAL 94 : c IVAMAT 95 : MPTVAL=IVAMAT 96 : c 97 : C++++ Activation au cas ou ++++++++++++++++++++++++++++++ 98 : SEGACT MCOORD 99 : 100 : C++++ Recup + Activation de la geometrie ++++++++++++++++ 101 : MELEME= IMAMOD 102 : c SEGACT MELEME deja actif 103 : 104 : C++++ RECUP DES INFOS EF ++++++++++++++++++++++++++++++++ 105 : c + OBTENUES PAR ELQUOI DANS MASSE1 106 : C segment INFO OBTENUES PAR ELQUOI DANS MASSE1 mais supprimé dans MASSE1 107 : MELE = NEFMOD 107 : call elquoi(MELE,0,4,IPINF,IMODEL) 108 : INFO = IPINF 109 : c MELE = INFELL(1) 110 : c NBPGA2= INFELL(2) 111 : c NBPGAU= INFELL(3) 112 : c NBPGAU= INFELL(4) 113 : c ICARA = INFELL(5) 114 : NGAU1 = INFELL(6) 115 : c LW = INFELL(7) 116 : LRE = INFELL(9) 117 : LHOOK = INFELL(10) 118 : MINT1 = INFELL(11) 119 : segact,MINT1 120 : MINT2 = INFELL(12) 121 : if(MINT2.ne.0) segact,MINT2 122 : MFR = INFELL(13) 123 : IELE = INFELL(14) 124 : NDDL = INFELL(15) 125 : NSTRS = INFELL(16) 126 : c write(ioimp,*)'-> MASSX infell',(infell(iou),iou=1,16) 127 : 128 : c + AUTRES INFOS 129 : C nbre de noeuds par element 130 : NBNN1 = NUM(/1) 131 : C nbre d elements 132 : NBEL1 = NUM(/2) 133 : 134 : c REM: pour se passer du dimensionnement du nbre d'enrichissement dans 135 : c elquoi et le realiser localement , on pourrait ecrire: 136 : c LRE = NDDLMAX*NBNN1 137 : c NDDL= NDDLMAX 138 : 139 : C sous decoupage et points de Gauss de l element geometrique de base 140 : CYT if(MELE.eq.263) then 141 : if(MELE.eq.263.or.MELE.eq.264) then 142 : c NGAU2 = 4 143 : NGAU2 = MINT2.POIGAU(/1) 144 : c NBSSEF = (NGAU1/NGAU2)**0.5 145 : endif 146 : c write(ioimp,*) 'dim de MINT2=6,',(MINT2.SHPTOT(/2)),(MINT2.SHPTOT(/3)) 147 : c write(ioimp,*) 'MINT2',(MINT2.QSIGAU(iou),iou=1,NGAU) 148 : c segdes,MINT2 149 : 150 : 151 : c nbre maxi de fonction de forme par noeud (fonction std comprise) 152 : c NBNI = NDDL/IDIM inutile! 153 : 154 : 155 : C++++ Recup des infos d enrichissement +++++++++++++++++++ 156 : c recup du MCHEX1 d'enrichissement 157 : NOBMO1 = IVAMOD(/1) 158 : if(NOBMO1.ne.0) then 159 : do iobmo1=1,NOBMO1 160 : if((TYMODE(iobmo1)).eq.'MCHAML') then 161 : MCHEX1 = IVAMOD(iobmo1) 162 : segact,MCHEX1 163 : if((MCHEX1.TITCHE).eq.'ENRICHIS') then 164 : MCHAM1 = MCHEX1.ICHAML(1) 165 : segact,MCHAM1 166 : segdes,MCHEX1 167 : goto 1000 168 : endif 169 : segdes,MCHEX1 170 : endif 171 : enddo 172 : write(ioimp,*) 'Le modele est vide (absence d enrichissement)' 173 : * return 174 : else 175 : write(ioimp,*) 'Aucun MCHAML enrichissement dans le Modele' 176 : * return 177 : endif 178 : 179 : 1000 continue 180 : c niveau d enrichissement(s) du modele (ddls std u exclus) 181 : c NBENR1= 0 si std, 1 si H seul, 2 si H+F1, 3 si H+F1+F2, etc... 182 : if(NOBMO1.ne.0) then 183 : NBENR1= MCHAM1.IELVAL(/1) 184 : else 185 : NBENR1 = 0 186 : endif 187 : c write(ioimp,*) 'niveau d enrichissement(s) du modele',NBENR1 188 : 189 : 190 : C********************************************************* 191 : c 192 : c PHASE 2 : PREPARATION DES OBJETS RESULTATS 193 : c 194 : C********************************************************* 195 : C 196 : C++++ RECHERCHE DES NOMS D'INCONNUES ET DES DUAUX 197 : c MOTINC et MODUA en dur pour l instant 198 : 199 : C++++ REMPLISSAGE DE DESCR de taille maxi 200 : c (on ajustera en enlevant si besoin) 201 : NLIGRP = LRE 202 : NLIGRD = LRE 203 : SEGINI DESCR 204 : IPDSCR = DESCR 205 : C 206 : KINC = 0 207 : C----->> BOUCLE SUR LES fonctions de Forme 208 : DO IENR=1,NBNIMAX 209 : C------->> BOUCLE SUR LES NOEUDS 210 : DO I=1,NBNN1 211 : C--------->> BOUCLE SUR LA DIMENSION 212 : DO KDIM=1,IDIM 213 : KINC = KINC + 1 214 : LISINC(KINC) = MOTINC((3*(IENR-1))+KDIM) 215 : LISDUA(KINC) = MOTDUA((3*(IENR-1))+KDIM) 216 : CYT LISINC(KINC) = MOTINC((IDIM*(IENR-1))+KDIM) 217 : CYT LISDUA(KINC) = MOTDUA((IDIM*(IENR-1))+KDIM) 218 : NOELEP(KINC) = I 219 : NOELED(KINC) = I 220 : ENDDO 221 : C--------------| 222 : ENDDO 223 : C-------------| 224 : ENDDO 225 : C------------| 226 : 227 : IF(KINC.NE.LRE) THEN 228 : WRITE(*,*) 'IL Y A UNE ERREUR DANS DESCR' 229 : WRITE(*,*) 'KINC , LRE = ',KINC,LRE 230 : RETURN 231 : ENDIF 232 : C 233 : SEGDES DESCR 234 : C 235 : C********************************************************* 236 : C INITIALISATIONS... 237 : C********************************************************* 238 : C 239 : c preparation des tables avec: 240 : 241 : if(NBENR1.ne.0) then 242 : do ienr=1,NBENR1 243 : c -le nombre d'inconnues de chaque sous-zone 244 : c determinee depuis le nombre de fonction de forme 245 : c ienr= 1: U+H(1+1=2), 2: U+H+F1(2+4=6), 3: U+H+F1+F2(6+4=10) 246 : nbniJ = 2 + ((ienr-1)*4) 247 : MLRE(1+ienr) = IDIM*NBNN1*nbniJ 248 : c -LOCIRI 249 : LOCIRI(5,1+ienr)= NIFOUR 250 : enddo 251 : endif 252 : C Tables + longues car 1er indice correspond au fontion de forme std 253 : MLRE(1) = IDIM*NBNN1*1 254 : LOCIRI(5,1)= NIFOUR 255 : 256 : if(NBENR1.lt.(NBENRMAX+1)) then 257 : do ienr=(NBENR1+1),(NBENRMAX) 258 : MLRE(1+ienr) = 0 259 : enddo 260 : endif 261 : c 262 : c ...DU SEGMENT WRK1 263 : NBENRMA2 = NBENRMAX 264 : NBBB = NBNN1 265 : SEGINI,WRK1 266 : 267 : c ...DU SEGMENT WRK2 268 : c NBNO = NBNI 269 : NBNO = LRE/IDIM 270 : c ici, pour le calcul de la masse, BGENE ne designe pas B mais N => lhook=2 271 : CTY LHOOK=2 272 : LHOOK=IDIM 273 : SEGINI,WRK2 274 : 275 : c ...DU SEGMENT MVELCH 276 : NV1 = NMATT 277 : SEGINI,MVELCH 278 : C 279 : c ...DU SEGMENT MRACC 280 : NBENRMA2 = NBENRMAX 281 : NBI = NBNN1 282 : segini,MRACC 283 : C 284 : C 285 : C********************************************************* 286 : C 287 : C>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> BOUCLE SUR LES ELEMENTS 288 : C 289 : DO 2000 J=1,NBEL1 290 : c write(ioimp,*) '========= EF',J,' /NBEL1 =========' 291 : C 292 : C 293 : C********************************************************* 294 : C POUR CHAQUE ELEMENT, ... 295 : C********************************************************* 296 : C 297 : C++++ ON RECUPERE LES COORDONNEES DES NOEUDS DE L ELEMENT 297 : CALL DOXE(XCOOR,IDIM,NBNN1,NUM,J,XE) 298 : C 299 : c 299 : c CALL ZEROI(TLREEL,NBENRMA2,NBI) 300 : C++++ NBENRJ = niveau d enrichissement de l element ++++ 301 : C =0 si EF std =1 si U+H =2 si U+H+F1 =3 si U+H+F1+F2 302 : NBENRJ=0 303 : if(NBENR1.ne.0) then 304 : do IENR=1,NBENR1 305 : MELVA1 = MCHAM1.IELVAL(IENR) 306 : segact,MELVA1 307 : do I=1,NBNN1 308 : mlree1 = MELVA1.IELCHE(I,J) 309 : C on en profite pour remplir MRACC table de raccourcis pour cet element 310 : TLREEL(IENR,I) = mlree1 311 : if(mlree1.ne.0) then 312 : NBENRJ = max(NBENRJ,IENR) 313 : C et on active la listreel 314 : segact,mlree1 315 : endif 316 : enddo 317 : segdes,MELVA1 318 : enddo 319 : endif 320 : if(NBENRMA2.gt.NBENR1) then 321 : do IENR2=(NBENR1+1),NBENRMA2 322 : do I=1,NBNN1 323 : TLREEL(IENR2,I) = 0 324 : enddo 325 : enddo 326 : endif 327 : C 328 : 329 : c++++ quelques indices pour dimensionner au plus juste 330 : c nbre total de ddl de l'élément considéré 331 : NLIGRD = MLRE(1+NBENRJ) 332 : NLIGRP = MLRE(1+NBENRJ) 333 : c nbre fonction de forme=((Ni_std+Ni_enrichi)*nbnoeud)=Ninconnues/idim 334 : NBNI = (MLRE(1+NBENRJ)) / IDIM 335 : 336 : c write(ioimp,*) 'EF',J,' NBENRJ=',NBENRJ, 337 : c &'donc',NLIGRD,' ddls et ',NBNI,' fonctions de forme' 338 : c 339 : 340 : c++++ Selon le niveau d enrichissement, 341 : c 342 : C + On copie cet element a la geo correspondante 343 : IPT1 = LOCIRI(1,1+NBENRJ) 344 : c si elle n existe pas, il faut la creer 345 : if(IPT1.eq.0) then 346 : NBNN =NBNN1 347 : NBELEM=1 348 : NBSOUS=0 349 : NBREF=0 350 : segini,IPT1 351 : IPT1.ITYPEL = ITYPEL 352 : LOCIRI(1,1+NBENRJ)=IPT1 353 : c si elle existe, il faut l agrandir 354 : else 355 : NBNN =NBNN1 356 : NBELEM = (IPT1.NUM(/2)) + 1 357 : NBSOUS=0 358 : NBREF=0 359 : segadj,IPT1 360 : endif 361 : c copie en cours 362 : do I=1,NBNN1 363 : IPT1.NUM(I,NBELEM) = NUM(I,J) 364 : enddo 365 : 366 : C + On crée le descr qui va bien si il n existe pas 367 : DES1 = LOCIRI(3,1+NBENRJ) 368 : if(DES1.eq.0) then 369 : segini,DES1=DESCR 370 : c write(ioimp,*)' on ne garde que les',MLRE(1+NBENRJ), 371 : c &' 1eres inconnues' 372 : c NLIGRP = MLRE(1+NBENRJ) fait + haut (ligne 313) 373 : c NLIGRD = MLRE(1+NBENRJ) 374 : segadj,DES1 375 : LOCIRI(3,1+NBENRJ) = DES1 376 : segdes,DES1 377 : endif 378 : 379 : C + On en profite pour créer IMATRI si il n existe pas 380 : xMATR1 = LOCIRI(4,1+NBENRJ) 381 : if(xMATR1.eq.0) then 382 : NELRIG = 1 383 : segini,xMATR1 384 : LOCIRI(4,1+NBENRJ) = xMATR1 385 : c si il existe, il faut l agrandir 386 : else 387 : NELRIG = xmatr1.RE(/3) + 1 388 : segadj,xMATR1 389 : endif 390 : 391 : 392 : C++++ ON MET A ZERO LA SOMME QUE L ON VA CALCULER 392 : c CALL ZERO(RINT,LRE,NLIGRP) 392 : * CALL ZERO(RINT,NLIGRD,NLIGRP) 393 : c |-> impossible car ZERO considere qu il sagit de la dimension de SHPWRK 393 : CALL ZEROX(RINT,NLIGRD,NLIGRP,LRE) 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 n est pas enrichi 398 : c car cela n'apprte rien de plus 399 : c On pourrait par ex. utiliser MINT2 = infell(12) qui contiendrait 400 : c le segment d'integration de l'EF std (QUA4 par ex.) 401 : if((NBENRJ.eq.0).and.(MINT2.ne.0)) then 402 : MINTE = MINT2 403 : NBPGAU= NGAU2 404 : else 405 : MINTE = MINT1 406 : NBPGAU= NGAU1 407 : endif 408 : C 409 : C********************************************************* 410 : C 411 : C>>>>>>>>>> BOUCLE SUR LES POINTS DE GAUSS 412 : C 413 : DO 2500 KGAU=1,NBPGAU 414 : C 415 : C********************************************************* 416 : C Initialisation à 0 417 : C********************************************************* 418 : 419 : c ZERO ne serait-il pas facultatif? 419 : * CALL ZERO(SHPWRK,6,NBNO) 419 : CALL ZERO(SHPWRK,6,NBNI) 420 : 421 : XGAU = 0. 422 : YGAU = 0. 423 : ZGAU = 0. 424 : C 425 : i6zz = 3 426 : IF (IDIM.EQ.3) i6zz = 4 427 : C 428 : c do ienr7=1,NBENRMAX 429 : do ienr7=1,NBENRJ 430 : do inod7=1,NBNN1 431 : c do i6=1,6 432 : CYT do i6=1,3 433 : do i6=1,i6zz 434 : LV7WRK(ienr7,1,i6,inod7) = 0 435 : LV7WRK(ienr7,2,i6,inod7) = 0 436 : enddo 437 : enddo 438 : enddo 439 : 440 : c********************************************************* 441 : c Calcul des fonction de forme std dans repere local 442 : c********************************************************* 443 : 444 : ccccc BOUCLE SUR LES NOEUDS ccccccccccccccccccccccccccccc 445 : c (et donc sur les Ni std) 446 : DO 2510 I=1,NBNN1 447 : 448 : C++++ Calcul des Ni std 449 : c (rappel: 1:Ni 2:Ni,qsi 3:Ni,eta avec i=1,4) 450 : SHPWRK(1,I) = SHPTOT(1,I,KGAU) 451 : SHPWRK(2,I) = SHPTOT(2,I,KGAU) 452 : SHPWRK(3,I) = SHPTOT(3,I,KGAU) 453 : IF (IDIM.EQ.3) SHPWRK(4,I) = SHPTOT(4,I,KGAU) 454 : 455 : C++++ Calcul des coordonnees dans le repere global 456 : XGAU = XGAU + (SHPWRK(1,I)*XE(1,I)) 457 : YGAU = YGAU + (SHPWRK(1,I)*XE(2,I)) 458 : IF (IDIM.EQ.3) ZGAU = ZGAU + (SHPWRK(1,I)*XE(3,I)) 459 : 460 : 2510 CONTINUE 461 : ccccc fin de BOUCLE SUR LES NOEUDS ccccccccccccccccccccccc 462 : 463 : 464 : 465 : c********************************************************* 466 : c Passage des fonctions de forme std dans repere global 467 : c********************************************************* 468 : c 469 : C++++ CALCUL DES Ni,x Ni,y (i=1,4) + CALCUL DE det(J) 469 : CALL JACOBI(XE,SHPWRK,IDIM,NBNN1,DJAC) 470 : c 471 : C PRISE EN COMPTE DU POIDS D'INTEGRATION 472 : DJAC = ABS(DJAC) * POIGAU(KGAU) 473 : 474 : c CAS DES CONTRAINTES PLANES 475 : C recuperation de l'epaisseur: DIM3 donnée facultative du materiau 476 : IF (IFOUR.EQ.-2) THEN 477 : MPTVAL=IVACAR 478 : IF (IVACAR.NE.0) THEN 479 : MELVAL=IVAL(1) 480 : IF (MELVAL.NE.0) THEN 481 : IGMN=MIN(KGAU,VELCHE(/1)) 482 : IBMN=MIN(J,VELCHE(/2)) 483 : DIM3=VELCHE(IGMN,IBMN) 484 : ELSE 485 : DIM3=1.D0 486 : ENDIF 487 : ENDIF 488 : DJAC=DJAC*DIM3 489 : MPTVAL=IVAMAT 490 : ENDIF 491 : 492 : 493 : c********************************************************* 494 : c Si on est pas du tout enrichi on peut sauter une grosse partie 495 : c********************************************************* 496 : if(NBENRJ.eq.0) goto 2999 497 : 498 : c********************************************************* 499 : c Calcul des level set + leurs derivees dans repere global 500 : c********************************************************* 501 : 502 : ccccc BOUCLE SUR LES enrichissements ccccccccccccccccccc 503 : do 2520 ienr=1,NBENRJ 504 : 505 : c MELVA1=MCHAM1.IELVAL(IENR) 506 : c segact,MELVA1 507 : 508 : ccccc BOUCLE SUR LES NOEUDS ccccccccccccccccccccccccccc 509 : DO 2521 I=1,NBNN1 510 : 511 : C++++ Le I eme noeud est-il ienr-enrichi? 512 : mlree1= TLREEL(IENR,I) 513 : if(mlree1.eq.0) goto 2521 514 : 515 : 516 : c Calcul du repere local de fissure(=PSI,PHI) 517 : c (rappel: 1,1:psi 1,2:phi 2,1 psi,x ...etc...) 518 : do 2522 inode=1,NBNN1 519 : c pour le H-enrichissement, on n a pas gardé PSI (inutile) 520 : if(ienr.ne.1) then 521 : c valeur de PSI au inode^ieme noeud 522 : xpsi1 = mlree1.PROG(inode) 523 : c qu on multiplie par la valeur de Ni^std au pt de G considéré 524 : LV7WRK(ienr,1,1,I)= LV7WRK(ienr,1,1,I) 525 : & + (SHPWRK(1,inode)*xpsi1) 526 : LV7WRK(ienr,1,2,I)= LV7WRK(ienr,1,2,I) 527 : & + (SHPWRK(2,inode)*xpsi1) 528 : LV7WRK(ienr,1,3,I)= LV7WRK(ienr,1,3,I) 529 : & + (SHPWRK(3,inode)*xpsi1) 530 : IF (IDIM.EQ.3) LV7WRK(ienr,1,4,I)= LV7WRK(ienr,1,4,I) 531 : & + (SHPWRK(4,inode)*xpsi1) 532 : c valeur de PHI au inode^ieme noeud 533 : xphi1 = mlree1.PROG(NBNN1+inode) 534 : else 535 : xphi1 = mlree1.PROG(inode) 536 : endif 537 : LV7WRK(ienr,2,1,I)= LV7WRK(ienr,2,1,I) 538 : & + (SHPWRK(1,inode)*xphi1) 539 : LV7WRK(ienr,2,2,I)= LV7WRK(ienr,2,2,I) 540 : & + (SHPWRK(2,inode)*xphi1) 541 : LV7WRK(ienr,2,3,I)= LV7WRK(ienr,2,3,I) 542 : & + (SHPWRK(3,inode)*xphi1) 543 : IF (IDIM.EQ.3) LV7WRK(ienr,2,4,I)= LV7WRK(ienr,2,4,I) 544 : & + (SHPWRK(4,inode)*xphi1) 545 : 2522 continue 546 : 547 : 2521 continue 548 : ccccc fin de BOUCLE SUR LES NOEUDS ccccccccccccccccccccccc 549 : 550 : 551 : 2520 CONTINUE 552 : ccccc fin de BOUCLE SUR LES enrichissements cccccccccccccccc 553 : 554 : c on a construit 555 : C LV7WRK(ienr, PSI/PHI, valeur/deriveeparqsi/pareta, iNOEUD) 556 : 557 : 558 : c********************************************************* 559 : c Ajout des fonctions de forme d enrichissement 560 : c + leurs derivees dans repere global 561 : c********************************************************* 561 : CYT CALL SHAPX(XGAU,YGAU,MELE,LV7WRK,NBENRMAX,NBENRJ, 561 : CALL SHAPX(XGAU,YGAU,ZGAU,MELE,LV7WRK,NBENRMAX,NBENRJ, 562 : &TLREEL,J,SHPWRK,IRET) 563 : 564 : 565 : c retour a la partie commune aux EF enrichi et non enrichi 566 : 2999 continue 567 : 568 : C********************************************************* 569 : C CALCUL DE [N] 570 : C********************************************************* 571 : c ZERO ne serait-il pas facultatif? 571 : c call ZERO(BGENE,LHOOK,NLIGRP) 572 : KB=1 573 : C boucle sur tous les Ni 574 : DO 3001 II=1,NBNI 575 : 576 : BGENE(1,KB) = SHPWRK(1,II) 577 : BGENE(2,KB+1) = SHPWRK(1,II) 578 : 579 : IF(IDIM.EQ.3) BGENE(3,KB+2) = SHPWRK(1,II) 580 : 581 : KB = KB + IDIM 582 : CTY KB = KB + 2 583 : 584 : 3001 CONTINUE 585 : C 586 : c if(J.eq.1.and.KGAU.eq.1) then 587 : c write(ioimp,*) 'BGENE(1,..)=',(BGENE(1,iou),iou=1,2*NBNI) 588 : c write(ioimp,*) 'BGENE(2,..)=',(BGENE(2,iou),iou=1,2*NBNI) 589 : c endif 590 : 591 : C********************************************************* 592 : C CALCUL DE rho 593 : C********************************************************* 594 : c 595 : C on remplit le segment MVELCH => valmat(1)=RHO 596 : IF (IVAL(1).NE.0) THEN 597 : MELVAL=IVAL(1) 598 : IBMN=MIN(J ,VELCHE(/2)) 599 : IGMN=MIN(KGAU,VELCHE(/1)) 600 : VALMAT(1)=VELCHE(IGMN,IBMN) 601 : ELSE 602 : VALMAT(1)=0.D0 603 : ENDIF 604 : 605 : C PRISE EN COMPTE DE RHO 606 : DJAC=DJAC*VALMAT(1) 607 : c if(J.eq.1) write(ioimp,*) 'DJAC=',DJAC 608 : c 609 : c 610 : C********************************************************* 611 : C CALCUL DE rho.Nt.N 612 : C********************************************************* 613 : 613 : CALL ZEROX(REL,NLIGRD,NLIGRP,LRE) 614 : 614 : CALL NTNST(BGENE,DJAC,LRE,LHOOK,REL) 615 : c pour gagner du temps cpu il faudrait qqchose du type: 615 : c CALL NTNSTX(BGENE,DJAC,DDHOOK,NLIGRP,2,LRE,REL) 616 : c if(J.eq.1) then 617 : c write(ioimp,*) 'REL',(REL(1,iou),iou=1,8) 618 : c write(ioimp,*) 'REL',(REL(2,iou),iou=1,8) 619 : c write(ioimp,*) 'REL',(REL(3,iou),iou=1,8) 620 : c write(ioimp,*) 'REL',(REL(4,iou),iou=1,8) 621 : c endif 622 : 623 : 624 : C********************************************************* 625 : C CUMUL DANS RINT(II,JJ) 626 : C********************************************************* 627 : DO 4201 II=1,NLIGRD 628 : DO 4202 JJ=1,NLIGRP 629 : RINT(II,JJ) = RINT(II,JJ) + REL(II,JJ) 630 : 4202 CONTINUE 631 : 4201 CONTINUE 632 : c 633 : c 634 : 2500 CONTINUE 635 : C FIN DE BOUCLE SUR LES POINTS DE GAUSS <<<<<<<<<<<<<< 636 : C 637 : 638 : 639 : C********************************************************* 640 : c write(ioimp,*) 'C STOCKAGE DANS XMATRI de RE(II,JJ)' 641 : C et de XMATRI dans IMATRI 642 : C********************************************************* 643 : C 644 : * SEGINI,XMATRI 645 : * IMATR1.IMATTT(NELRIG) = XMATRI 646 : c 647 : DO 6001 II=1,NLIGRD 648 : 649 : DO 6002 JJ=1,NLIGRP 650 : 651 : xmatr1.RE(II,JJ,nelrig) = RINT(II,JJ) 652 : 653 : c si NON enrichi, on met tout a 0 sauf la diago 654 : c if(mlree1.eq.0) RE(II,JJ) = 0. inutile 655 : 656 : 6002 CONTINUE 657 : c if(J.eq.1) write(ioimp,*) 'RE_',II,'=',(RE(II,jou),jou=1,NLIGRP) 658 : 659 : 6001 CONTINUE 660 : c 661 : c quelques desactivations 662 : CYT SEGDES,XMATRI 663 : do IENR=1,NBENR1 664 : do I=1,NBNN1 665 : mlree1 = TLREEL(IENR,I) 666 : if(mlree1.ne.0) segdes,mlree1 667 : enddo 668 : enddo 669 : c 670 : c 671 : c 672 : 2000 CONTINUE 673 : C FIN DE BOUCLE SUR LES ELEMENTS <<<<<<<<<<<<<<<<<<<<< 674 : c 675 : c 676 : 677 : C********************************************************* 678 : C SUPPRESSION ET DESACTIVATION DE SEGMENTS 679 : C********************************************************* 680 : C 681 : SEGSUP,WRK1,WRK2,MVELCH 682 : SEGSUP,MRACC 683 : 684 : segdes,MELEME 685 : segdes,MINT1 686 : if(MINT2.ne.0) segdes,MINT2 687 : c desactivation des maillages et segment imatri 688 : do ienr=1,(1+NBENR1) 689 : IPT1 = LOCIRI(1,ienr) 690 : if(IPT1.ne.0) segdes,IPT1 691 : xMATR1 = LOCIRI(4,ienr) 692 : if(xMATR1.ne.0) segdes,xMATR1 693 : enddo 694 : 695 : c SEGDES dans RIGI1 = IMODEL 696 : C 697 : RETURN 698 : END 699 : 700 : 701 :
© Cast3M 2003 - All rights reserved.
Disclaimer