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