1 : C ACCRO2 SOURCE BP208322 13/12/20 21:15:01 7890 2 : SUBROUTINE ACCRO2 3 : C======================================================================== 4 : C Cree la matrice de liaison entre le champ u et w 5 : C avec une formulation faible a 3 champs : u w et t 6 : C FAIBLE => On utilise les fonctions de forme 7 : c 8 : C Option PENA et STAB pour créer les matrices de penalite Kww 9 : C et de stabilisation Ktt (utile si LATIN utilisée) 10 : c 11 : C En sortie, on a Kfaibl = [Ktt Ktu Ktw ; Kut 0 0 ; Kwt 0 Kww] 12 : C [0 Ktu Ktw ; sym. ] = \int t * (u-w) ds 13 : c 14 : c ajout des enrichissements XFEM 15 : c 16 : C Creation : Benoit Trolle (SNCF-Lamcos), avril 2012 17 : C Integration dans Cast3m : BP, decembre 2012 18 : C======================================================================== 19 : 20 : IMPLICIT INTEGER(I-N) 21 : IMPLICIT REAL*8 (A-H,O-Z) 22 : 23 : -INC CCREEL 24 : -INC CCOPTIO 25 : -INC SMCOORD 26 : -INC SMELEME 27 : -INC SMRIGID 28 : -INC SMLMOTS 29 : -INC CCHAMP 30 : -INC CCGEOME 31 : -INC SMMODEL 32 : -INC SMINTE 33 : -INC SMCHAML 34 : -INC SMLREEL 35 : 36 : 37 : C Segment contenant les informations sur un element 38 : SEGMENT INFO 39 : INTEGER INFELL(JG) 40 : ENDSEGMENT 41 : 42 : POINTEUR INFO1.INFO,INFO2.INFO 43 : POINTEUR MCHEX1.MCHELM 44 : external shape 45 : DATA EPSI/1.D-9/ 46 : DIMENSION ICOR(6),IMEL(6),imtt(10) 47 : DIMENSION QSI(3),XPO(3) 48 : 49 : C INITIALISATION DES INCONNUES obligatoires et facultatives 50 : PARAMETER (NOBL=3,NFAC=9,NMOCLE=3) 51 : CHARACTER*4 DDLOBL(NOBL),DDLFAC(NFAC),MODDL,MODDL2 52 : CHARACTER*4 DUAOBL(NOBL),DUAFAC(NFAC),MOCLE(NMOCLE) 53 : CHARACTER*4 DDLOB2(NOBL),DDLFA2(NOBL) 54 : CHARACTER*4 DUAOB2(NOBL),DUAFA2(NOBL) 55 : 56 : DATA DDLOBL/'UX ','UY ','UZ '/ 57 : DATA DDLFAC/'AX ','AY ','AZ ', 58 : >'B1X ','B1Y ','B1Z ', 59 : >'B2X ','B2Y ','B2Z '/ 60 : DATA DUAOBL/'FX ','FY ','FZ '/ 61 : DATA DUAFAC/'FAX ','FAY ','FAZ ', 62 : >'FB1X','FB1Y','FB1Z', 63 : >'FB2X','FB2Y','FB2Z'/ 64 : 65 : c nom de composantes speciales 66 : DATA DDLOB2/'UX ','UY ','UZ '/ 67 : DATA DDLFA2/'AX ','AY ','AZ '/ 68 : DATA DUAOB2/'FX ','FY ','FZ '/ 69 : DATA DUAFA2/'FAX ','FAY ','FAZ '/ 70 : 71 : DATA MOCLE /'PENA','STAB','NOMU'/ 72 : 73 : c----------------------------------------------------- 74 : c Segment de travail 75 : c----------------------------------------------------- 76 : SEGMENT MTRAV 77 : REAL*8 SHPP2(6,NBNN2) 78 : REAL*8 XE2(3,NBNN2),XEL2(3,NBNN2),NGENE2(1,NBNN2) 79 : REAL*8 XE3(3,NBNN2) 80 : REAL*8 SHPP1(6,NBNN1) 81 : REAL*8 XE1(3,NBNN1),XEL1(3,NBNN1) 82 : REAL*8 NGENE1(1,NBNN1) 83 : REAL*8 REL(nbnn2,nbnn2),REL1(nbnn2,nbnn1) 84 : REAL*8 RINT(nbnn2,nbnn2) 85 : REAL*8 RINT0(nbnn2,ngau2*nbnn1),RINT1(nbnn2,ngau2*nbnn1) 86 : REAL*8 RINT2(nbnn2,ngau2*nbnn1),RINT3(nbnn2,ngau2*nbnn1) 87 : INTEGER ITYMAT(2,4*nbnn1*ngau2) 88 : INTEGER ITAMYT(ngau2*NBNN1,4) 89 : INTEGER PRELEM(3*nbnn1*ngau2) 90 : INTEGER ElemPG(ngau2) 91 : ENDSEGMENT 92 : 93 : c tableaux comptant le nbre d EF de chaque ddl 94 : PARAMETER(NDDLMAX=6) 95 : INTEGER NELDDL(NDDLMAX) 96 : 97 : if(iimpi.ge.2) then 98 : write(ioimp,*) '-----------------------------------------------' 99 : write(ioimp,*) ' ENTREE dans ACCRO2' 100 : write(ioimp,*) '-----------------------------------------------' 101 : endif 102 : 103 : SEGACT MCOORD 104 : 105 : c----------------------------------------------------- 106 : c LECTURE DES MOT CLES PENAlite et STABilistation 107 : c----------------------------------------------------- 108 : RLATIN=0.d0 109 : XISTAB=0.d0 110 : INOMU=0 111 : 700 ICLE=0 111 : CALL LIRMOT(MOCLE,NMOCLE,ICLE,0) 112 : IF (ICLE.EQ.1) GOTO 701 113 : IF (ICLE.EQ.2) GOTO 702 114 : IF (ICLE.EQ.3) THEN 115 : INOMU=1 116 : GOTO 700 117 : ENDIF 118 : GOTO 799 119 : 120 : c LECTURE DE LA DIRECTION DE RECHERCHE (FLOTTANT) 120 : 701 CALL LIRREE(RLATIN,1,IRETOU) 121 : IF (IERR.NE.0) RETURN 122 : if(iimpi.ge.2) write(ioimp,*) 'On a lu la raideur :',RLATIN 123 : GOTO 700 124 : 125 : c LECTURE DE LA STABILISATION (FLOTTANT) 125 : 702 CALL LIRREE(XISTAB,1,IRETOU) 126 : IF (IERR.NE.0) RETURN 127 : if(iimpi.ge.2) write(ioimp,*) 'On a lu la STABILISATION :',XISTAB 128 : GOTO 700 129 : 130 : 799 CONTINUE 131 : 132 : 133 : c----------------------------------------------------- 134 : c RECUPERATION DU MAILLAGE MASSIF 135 : c----------------------------------------------------- 135 : CALL LIROBJ ('MMODEL',IPMODL,1,IRETOU) 136 : IF(IERR.NE.0) RETURN 137 : MMODE1=IPMODL 138 : segact,MMODE1 139 : c Récupération du nombre de zones du modèle 140 : N1 = MMODE1.KMODEL(/1) 141 : if(N1.gt.1) write(ioimp,*) 'attention 1 seule zone a ce jour!' 142 : IMODE1 = MMODE1.KMODEL(1) 143 : segdes,MMODE1 144 : C Récupération du maillage et du numéro d'élément du modèle 145 : segact,IMODE1 146 : nele1 = IMODE1.NEFMOD 147 : IPT1 = IMODE1.IMAMOD 148 : SEGACT IPT1 149 : C Récupération du numéro d'élément du maillage, du nombre de noeuds et d'éléments 150 : iele1 = IPT1.itypel 151 : nbnn1 = IPT1.num(/1) 152 : nbel1 = IPT1.num(/2) 153 : c récupération des caractéristique EF IPT1 154 : cbp : on essaie d'abord avec INFELE 155 : mele1 = IMODE1.IMAMOD 156 : MINTE1 = IMODE1.INFELE(11) 157 : cbp : si echec, on retente avec elquoi 158 : if(mele1.eq.0.or.MINTE1.eq.0) then 158 : call elquoi(nele1,0,3,IPTR1,IMODE1) 159 : INFO = IPTR1 160 : segact,INFO 161 : mele1 = INFELL(1) 162 : MINTE1 = INFELL(11) 163 : segdes,INFO 164 : endif 165 : segdes,IMODE1 166 : cbp : si toujours pb, erreur... 167 : if(mele1.eq.0.or.MINTE1.eq.0) then 168 : write(ioimp,*) 'donnees relatives a l element fini inconnues' 168 : call erreur(591) 169 : endif 170 : 171 : 172 : c----------------------------------------------------- 173 : c RECUPERATION DU MAILLAGE INTERFACE 174 : c----------------------------------------------------- 175 : IPMAI2 = 0 175 : CALL LIROBJ ('MMODEL',IPMODL,0,IRETOU) 176 : IF(IERR.NE.0) RETURN 177 : 178 : C -DANS LE CAS OU ON A UN MODELE EN ENTREE 179 : IF (IRETOU.EQ.1)THEN 180 : MMODE2=IPMODL 181 : segact,MMODE2 182 : N2 = MMODE2.KMODEL(/1) 183 : if(N2.gt.1) write(ioimp,*) 'attention 1 seule zone a ce jour!' 184 : IMODE2 = MMODE2.KMODEL(1) 185 : segdes,MMODE2 186 : segact,IMODE2 187 : nele2 = IMODE2.NEFMOD 188 : IPT2 = IMODE2.IMAMOD 189 : SEGACT IPT2 190 : c pour l'instant on dit que nele = iele (marche pour iele entre 2 et 26, voir bdata.eso) 191 : iele2 = IPT2.itypel 192 : nbnn2 = IPT2.num(/1) 193 : nbel2 = IPT2.num(/2) 194 : 195 : c recuperation des caracteritiques de l'element 196 : cbp : on essaie d'abord avec INFELE 197 : MINTE2 = IMODE2.INFELE(11) 198 : ngau2 = IMODE2.INFELE(6) 199 : cbp : si echec, on retente avec elquoi 200 : if(MINTE2.eq.0.or.ngau2.eq.0) then 201 : c write(ioimp,*) 'appel elquoi',iele2,nele2 201 : call elquoi(nele2,0,3,IPTR2,IMODE2) 202 : INFO = IPTR2 203 : segact,INFO 204 : MINTE2 = INFELL(11) 205 : ngau2 = INFELL(6) 206 : segdes,INFO 207 : endif 208 : segdes,IMODE2 209 : cbp : si toujours pb, erreur... 210 : if(mele1.eq.0.or.MINTE1.eq.0) then 211 : write(ioimp,*) 'donnees relatives a l element fini INTERFACE' 212 : & ,'inconnues' 212 : call erreur(591) 213 : endif 214 : 215 : C -DANS LE CAS OU ON A UN MAILLAGE EN ENTREE 216 : ELSE 216 : CALL LIROBJ('MAILLAGE',IPMAI2,1,IRETOU) 217 : IF(IERR.NE.0) RETURN 218 : IPT2 = IPMAI2 219 : SEGACT IPT2 220 : c pour l'instant on dit que nele = iele (marche pour iele entre 2 et 26, voir bdata.eso) 221 : iele2 = IPT2.itypel 222 : nele2 = iele2 223 : if (nele2.lt.2.or.nele2.gt.26) then 224 : write(ioimp,*)'element geometrique different de l element fini' 224 : call erreur(16) 225 : endif 226 : nbnn2 = IPT2.num(/1) 227 : nbel2 = IPT2.num(/2) 228 : ngau2 = -3 229 : c SEG2 230 : if (nele2.EQ.2) ngau2 = 2 231 : c SEG3 232 : if (nele2.EQ.3) ngau2 = 3 233 : c TRI3 234 : if (nele2.EQ.4) ngau2 = 1 235 : c TRI6 236 : if (nele2.EQ.6) ngau2 = 4 237 : c QUA4 238 : if (nele2.EQ.8) ngau2 = 4 239 : c QUA8 240 : if (nele2.EQ.10) ngau2 = 9 241 : if(ngau2.eq.-3) then 242 : write(ioimp,*)'nombre de points d integration inconnu pour', 243 : & 'ce type d element geometrique' 243 : call erreur(16) 244 : endif 245 : ENDIF 246 : 246 : call RESHPT(ngau2,nbnn2,iele2,nele2,0,IPTR2,IRET) 247 : MINTE2 = IPTR2 248 : segact,MINTE1,MINTE2 249 : 250 : 251 : c----------------------------------------------------- 252 : c RECHERCHE DU MCHAML ISSU MCHEX1 D ENRICHISSEMENT 253 : c----------------------------------------------------- 254 : MCHAM1=0 255 : NBENR2=0 256 : segact,IMODE1 257 : NOBMOD = IMODE1.IVAMOD(/1) 258 : IF (NOBMOD.NE.0) THEN 259 : DO 1002 iobmo1=1,NOBMOD 260 : if((IMODE1.TYMODE(iobmo1)).ne.'MCHAML') goto 1002 261 : MCHEX1 = IMODE1.IVAMOD(iobmo1) 262 : segact,MCHEX1 263 : if((MCHEX1.TITCHE).ne.'ENRICHIS') goto 1003 264 : MCHAM1 = MCHEX1.ICHAML(1) 265 : segact,MCHAM1 266 : NBENR2 = MCHAM1.IELVAL(/1) 267 : do ienr2=1,NBENR2 268 : MELVA1=MCHAM1.IELVAL(IENR2) 269 : if(MELVA1.ne.0) segact,MELVA1 270 : enddo 271 : 1003 continue 272 : segdes,MCHEX1 273 : 1002 CONTINUE 274 : ENDIF 275 : 276 : c------------------------------------- 277 : c INITIALISATION DES OBJETS DE TRAVAIL : MTRAV et ITYMAT 278 : c------------------------------------- 279 : 280 : segini,MTRAV 281 : 282 : *bp : NTYMAT = (U ou H ou HB1 ou HB1B2) * dim * nbre de configs possibles 283 : * NTYMAT = nbre de types de matrices 284 : * nbre de configs varie entre : 285 : * -on peut avoir tous les pt de Gauss dans 1 seul elements de structure 286 : * -jusqu'à chaque pt de Gauss dans un element distinct 287 : ity=0 288 : c boucle sur les enrichissements possibles U, H , HB1, HB1B2 289 : if(iimpi.ge.3) write(ioimp,*) 'ity , ino1 , ienr1' 290 : do ienr1=1,4 291 : c boucle sur les nombre de noeuds "structure" 292 : do ino1=(nbnn1),(ngau2*nbnn1) 293 : ity=ity+1 294 : ITYMAT(1,ity)=ino1 295 : ITYMAT(2,ity)=ienr1 296 : ITAMYT(ino1,ienr1)=ity 297 : if(iimpi.ge.3) write(ioimp,*) ity,ino1,ienr1 298 : enddo 299 : enddo 300 : NTYMAT = 4*(nbnn1*(ngau2-1)+1) 301 : if(iimpi.ge.3) write(ioimp,*) 'ity*idim=',ity,NTYMAT 302 : 303 : 304 : c-------------------------------------------------------------------- 305 : c INITIALISATION DU SEGMENT MRIGID 306 : c-------------------------------------------------------------------- 307 : NRIGEL = NTYMAT*idim 308 : segini,MRIGID 309 : IFORIG = IFOUR 310 : MTYMAT ='RIGIDITE' 311 : c -on prepare le meleme 312 : NBSOUS = 0 313 : NBREF = 0 314 : 315 : ityty=0 316 : c on initialise la taille matrice en fonction du type de matrice 317 : do ity=1,NTYMAT 318 : do iidim=1,idim 319 : 320 : ityty=ityty+1 321 : COERIG(ityty) = 1.D0 322 : 323 : * dim de la matrice RE elementaire = nbnoeud*TX + nbnoeud*AX + nbnoeud*AX 324 : nbno1 = ITYMAT(1,ity) 325 : nenr1 = ITYMAT(2,ity) 326 : if(nenr1.le.2) then 327 : NLIGRP = nbnn2 + nbnn2 + nbno1 328 : else 329 : NLIGRP = nbnn2 + nbnn2 + ((nenr1-1)*nbno1) 330 : endif 331 : NLIGRD = NLIGRP 332 : 333 : c -creation du MELEME 334 : NBNN = nbnn2 + nbno1 335 : NBELEM=0 336 : SEGINI,MELEME 337 : ITYPEL=28 338 : IRIGEL(1,ityty) = MELEME 339 : 340 : c -remplissage du DESCR 341 : SEGINI,DESCR 342 : IRIGEL(3,ityty) = DESCR 343 : iddl=0 344 : c remplissage des ddl LX de la fissure (ici appelé FX) 345 : c on iverse donc inconnue et duale []*FX=UX 346 : do ino2=1,nbnn2 347 : iddl=iddl+1 348 : if (nenr1.eq.1) then 349 : LISINC(iddl)=DUAOB2(iidim) 350 : LISDUA(iddl)=DDLOB2(iidim) 351 : else 352 : LISINC(iddl)=DUAFA2(iidim) 353 : LISDUA(iddl)=DDLFA2(iidim) 354 : endif 355 : NOELEP(iddl)=ino2 356 : NOELED(iddl)=ino2 357 : enddo 358 : c remplissage des ddl UX de la fissure (ici appelé WX) []*WX=TX 359 : do ino2=1,nbnn2 360 : iddl=iddl+1 361 : if (nenr1.eq.1) then 362 : LISINC(iddl)=DDLOB2(iidim) 363 : LISDUA(iddl)=DUAOB2(iidim) 364 : else 365 : LISINC(iddl)=DDLFA2(iidim) 366 : LISDUA(iddl)=DUAFA2(iidim) 367 : endif 368 : NOELEP(iddl)=ino2 369 : NOELED(iddl)=ino2 370 : enddo 371 : c remplissage des ddl de la structure 372 : if (nenr1.eq.1) then 373 : do ino1=1,nbno1 374 : iddl=iddl+1 375 : LISINC(iddl)=DDLOBL(iidim) 376 : LISDUA(iddl)=DUAOBL(iidim) 377 : NOELEP(iddl)=nbnn2+ino1 378 : NOELED(iddl)=nbnn2+ino1 379 : enddo 380 : else 381 : do ini1=1,(nenr1-1) 382 : do ino1=1,nbno1 383 : iddl=iddl+1 384 : LISINC(iddl)=DDLFAC(iidim+(3*(ini1-1))) 385 : LISDUA(iddl)=DUAFAC(iidim+(3*(ini1-1))) 386 : NOELEP(iddl)=nbnn2+ino1 387 : NOELED(iddl)=nbnn2+ino1 388 : enddo 389 : enddo 390 : endif 391 : if(iimpi.ge.3) write(ioimp,*) ityty,(LISINC(iou),iou=1,NLIGRP) 392 : SEGDES,DESCR 393 : 394 : c -initialisation du XMATRI 395 : NELRIG=0 396 : SEGINI,XMATRI 397 : IRIGEL(4,ityty) = XMATRI 398 : IRIGEL(5,ityty) = NIFOUR 399 : IRIGEL(6,ityty) = 0 400 : IRIGEL(7,ityty) = 0 401 : IRIGEL(8,ityty) = 0 402 : 403 : enddo 404 : enddo 405 : 406 : 407 : c---------------------------------------------------------------------- 408 : c 1. RECHERCHE DES ELEMENTS DE STRUCTURE CONTENANT DES POINTS DE GAUSS 409 : c DES ELEMENTS DE LA FISSURE 410 : c 2. REMPLISSAGE DU MELEME + MRIGID en FONCTION 411 : c---------------------------------------------------------------------- 412 : 413 : C 414 : c==== Boucle sur les elements de fissure ============================== 415 : DO 1100 iem2=1,nbel2 416 : 416 : call doxe(xcoor,idim,nbnn2,ipt2.num,iem2,xe2) 417 : nbenrj = 0 418 : nbnno2 = 0 418 : call ZERO(RINT,nbnn2,nbnn2) 418 : call ZERO(RINT0,nbnn2,ngau2*nbnn1) 418 : call ZERO(RINT1,nbnn2,ngau2*nbnn1) 418 : call ZERO(RINT2,nbnn2,ngau2*nbnn1) 418 : call ZERO(RINT3,nbnn2,ngau2*nbnn1) 419 : 420 : c======= Boucle sur les pt de G de fissure ============================ 421 : DO 1132 igau2=1,ngau2 422 : izok = 0 423 : c récupération des coordonnees du point de gauss dans le repère global 424 : XPO(1) = 0.D0 425 : XPO(2) = 0.D0 426 : XPO(3) = 0.D0 427 : DO ino=1,nbnn2 428 : XPO(1) = XPO(1) + (MINTE2.SHPTOT(1,ino,IGAU2)*xe2(1,ino)) 429 : XPO(2) = XPO(2) + (MINTE2.SHPTOT(1,ino,IGAU2)*xe2(2,ino)) 430 : IF(IDIM.EQ.3) then 431 : XPO(3) = XPO(3) + (MINTE2.SHPTOT(1,ino,IGAU2)*xe2(3,ino)) 432 : ENDIF 433 : ENDDO 434 : 435 : c if(iimpi.ge.5) write(ioimp,*) 'Fissure: element',iem2, 436 : c & ' point de Gauss',igau2,' x=',XPO(1),XPO(2) 437 : 438 : c---------- Boucle sur les elements de structure ---------------------- 439 : DO 1131 iem1=1,nbel1 440 : 441 : c si pas d'enrichissement, on travaille sur tous les elements 442 : if(MCHAM1.eq.0) goto 1133 443 : c on saute les elements non enrichi car a priori ne contiennent pas la fissure 444 : do ienr2=1,NBENR2 445 : MELVA1=MCHAM1.IELVAL(IENR2) 446 : if(MELVA1.ne.0) then 447 : do inode1=1,nbnn1 448 : if(MELVA1.IELCHE(inode1,iem1).ne.0) goto 1133 449 : enddo 450 : endif 451 : enddo 452 : goto 1131 453 : 1133 continue 454 : 455 : c recuperation des coordonnées des noeuds de IPT1 : xel1 (dans le repère x,y,z) 455 : call doxe(xcoor,idim,nbnn1,ipt1.num,iem1,xel1) 456 : 457 : c calcul des fonctions de formes de IPT1 au pt de Gauss de IPT2 457 : call QSIJS(xel1,iele1,nbnn1,idim,XPO,SHPP1,qsi,IRET) 458 : c if(iimpi.ge.5) write(ioimp,*) 'Structure: element',iem1, 459 : c & ' noeuds:',(ipt1.num(iou,iem1),iou=1,nbnn1) 460 : c if(iimpi.ge.5) write(ioimp,*) 461 : c & ' N_i=',(SHPP1(1,iou),iou=1,nbnn1) 462 : 463 : 464 : c test pour savoir si PG est dans EF de IPT1 465 : DO 1130 ino1=1,NBNN1 466 : if (SHPP1(1,ino1).LT.-1.01D-7) then 467 : go to 1131 468 : endif 469 : 1130 continue 470 : izok = 1 471 : c ON a trouvé l'élément de structure ! 472 : c write(ioimp,*) 'ON a trouve l element de structure !' 473 : c DETECTION DU TYPE D'ENRICHISSEMENT MAX DE CET ELEMENT = nbenrj 474 : DO 3001 IENR2=1,NBENR2 475 : MELVA1=MCHAM1.IELVAL(IENR2) 476 : IF(MELVA1.eq.0) GOTO 3001 477 : DO 3002 ino1=1,nbnn1 478 : MLREEL = MELVA1.IELCHE(ino1,iem1) 479 : c Test pour savoir si le noeud est enrichi 480 : IF(MLREEL.eq.0) GOTO 3002 481 : nbenrj=max(nbenrj,IENR2) 482 : 3002 continue 483 : 3001 continue 484 : 485 : if(iimpi.ge.3) write(ioimp,*) 'EF fissure ',iem2, 486 : & ' ptdeG ',igau2,' -> EF MASSIF ',iem1,' nbenrj=',nbenrj 487 : 488 : c Préparation du Remplissage du MELEME 489 : do 3003 ino1=1,nbnn1 490 : iino1=IPT1.NUM(ino1,iem1) 491 : if(nbnno2.gt.0) then 492 : do iino2=1,nbnno2 493 : if(PRELEM(iino2).eq.iino1) goto 3003 494 : enddo 495 : endif 496 : c nouveau noeud : on l'ajoute 497 : nbnno2=nbnno2+1 498 : PRELEM(nbnno2)=iino1 499 : 3003 continue 500 : if(iimpi.ge.3) then 501 : write(ioimp,*) ' PRELEM=',( PRELEM(iou),iou=1,nbnno2) 502 : write(ioimp,*) 'EF fissure ',iem2,' ptdeG ',igau2 503 : endif 504 : 505 : C CALCUL DE REL ET REL1 506 : 507 : c matrices des fonctions de forme de l'element associé à IPT2 508 : 509 : do ino2=1,NBNN2 510 : xe3(1,ino2) = sqrt( xe2(1,ino2)*xe2(1,ino2) 511 : 1 + xe2(2,ino2)*xe2(2,ino2) ) 512 : NGENE2(1,ino2)=MINTE2.SHPTOT(1,ino2,IGAU2) 513 : do ipp2=1,6 514 : SHPP2(ipp2,ino2)=MINTE2.SHPTOT(ipp2,ino2,IGAU2) 515 : enddo 516 : if(iimpi.ge.4) then 517 : write(ioimp,*) 'SHPP2',(SHPP2(iou,ino2),iou=1,6) 518 : write(ioimp,*) 'xe2',(xe2(iou,ino2),iou=1,2) 519 : endif 520 : enddo 521 : c calcul du jacobien 522 : if(idim.eq.2) then 523 : c calcul du jacobien en 2D pour un seg2 524 : dXdQsi=0.D0 525 : dYdQsi=0.D0 526 : DO i=1,NBNN2 527 : dXdQsi=dXdQsi+SHPP2(2,i)*XE2(1,i) 528 : dYdQsi=dYdQsi+SHPP2(2,i)*XE2(2,i) 529 : ENDDO 530 : c Dans le cas des seg2, poids des pts de gauss = 1 531 : DJAC2=SQRT(dXdQsi*dXdQsi+dYdQsi*dYdQsi) 532 : else 533 : dXdQsi=0.D0 534 : dYdQsi=0.D0 535 : dZdQsi=0.D0 536 : dXdEta=0.D0 537 : dYdEta=0.D0 538 : dZdEta=0.D0 539 : DO i=1,NBNN2 540 : dXdQsi=dXdQsi+SHPP2(2,i)*XE2(1,i) 541 : dXdEta=dXdEta+SHPP2(3,i)*XE2(1,i) 542 : dYdQsi=dYdQsi+SHPP2(2,i)*XE2(2,i) 543 : dYdEta=dYdEta+SHPP2(3,i)*XE2(2,i) 544 : dZdQsi=dZdQsi+SHPP2(2,i)*XE2(3,i) 545 : dZdEta=dZdEta+SHPP2(3,i)*XE2(3,i) 546 : ENDDO 547 : z = (dXdQsi*dYdEta-dXdEta*dYdQsi) 548 : x = (dYdQsi*dZdEta-dYdEta*dZdQsi) 549 : y = (dZdQsi*dXdEta-dZdEta*dXdQsi) 550 : c Dans le cas des tri3, poids du pts de gauss = 0.5 551 : DJAC2 = sqrt(x*x+y*y+z*z)*0.5 552 : endif 553 : 554 : c matrices des fonctions de forme de l'element associé à IPT1 555 : do ino1=1,NBNN1 556 : NGENE1(1,ino1)=SHPP1(1,ino1) 557 : enddo 558 : if(iimpi.ge.4) then 559 : write(ioimp,*) 'DJAC2=',DJAC2 560 : write(ioimp,*) 'NGENE1=',(NGENE1(1,iou),iou=1,NBNN1) 561 : write(ioimp,*) 'NGENE2=',(NGENE2(1,iou),iou=1,NBNN2) 562 : endif 563 : 564 : c calcul des matrices elementaires = produits des fonctions de forme 565 : DO II=1,nbnn2 566 : XGENE2 = NGENE2(1,II)*DJAC2 567 : c calcul de K_ww, Kwt et Ktt 568 : DO JJ=1,nbnn2 569 : jjj = JJ 570 : C CC = Jacobien * transposee(N)(II,i) * N(i,JJ) (somme sur i) 571 : CC = XGENE2*NGENE2(1,jjj) 572 : C REL est une matrice symetrique 573 : REL(II,JJ) = CC 574 : REL(JJ,II) = REL(II,JJ) 575 : ENDDO 576 : c calcul de K_tu 577 : DO JJ=1,nbnn1 578 : C CC = Jacobien * transposee(N)(II,i) * N(i,JJ) (somme sur i) 579 : CC=XGENE2*NGENE1(1,JJ) 580 : c write(ioimp,*) ' CC = ', CC 581 : c REL1(II,JJ)=-1.D0*CC 582 : REL1(II,JJ)=CC 583 : c REL(JJdec,II)=REL(II,JJdec) 584 : ENDDO 585 : if(iimpi.ge.4) write(ioimp,*) 586 : 1 'REL1(',II,',:)=',(REL1(II,iou),iou=1,NBNN1) 587 : ENDDO 588 : 589 : C CUMUL DANS RINT ET RINT1 590 : C CUMUL DANS RINT(II,JJ) pour Kww, Kwt, Ktw et Ktt 591 : c (4 fois la meme matrices = N^{w T} * N^{w} = N^{w T} * N^{t} = ...) 592 : DO II=1,nbnn2 593 : DO JJ=1,nbnn2 594 : RINT(II,JJ) = RINT(II,JJ) + REL(II,JJ) 595 : ENDDO 596 : ENDDO 597 : 598 : C CUMUL DANS RINT1(II,JJ) pour Kut ( = N^{t T} * N^{u}) 599 : C INTEGRATION POSSIBLEMENT INCOMPATIBLE POUR KUT ! 600 : c on peut prendre le 1er enrichissement (ou le standard) car on a 601 : c toujours les memes noeuds au debut (sans tenir compte des B1X) 602 : 603 : c Stockage de l'element iem1 pour les differents PG 604 : elemPG(igau2) = iem1 605 : 606 : DO 4205 JJ=1,nbnn1 607 : iino1=IPT1.NUM(JJ,iem1) 608 : 609 : do JJ1 = 1,nbnno2 610 : if(PRELEM(JJ1).eq.iino1) goto 4207 611 : enddo 612 : write(IOIMP,*) 'On n a pas trouvé le noeud ',iino1 612 : call ERREUR(5) 613 : 4207 continue 614 : itoto = 0 615 : 616 : c UX 617 : DO II=1,nbnn2 618 : RINT0(II,JJ1) = RINT0(II,JJ1) + REL1(II,JJ) 619 : ENDDO 620 : c AX 621 : DO II=1,nbnn2 622 : MELVA1 = MCHAM1.IELVAL(1) 623 : MLREEL = MELVA1.IELCHE(JJ,iem1) 624 : if(MLREEL.ne.0) then 625 : RINT1(II,JJ1) = RINT1(II,JJ1) + REL1(II,JJ) 626 : c write(ioimp,*) ' RINT1(',II,JJ1,') = ', RINT1(II,JJ1) 627 : endif 628 : ENDDO 629 : 630 : 631 : c B1X 632 : IF(nbenrj.ge.2) THEN 633 : MELVA1 = MCHAM1.IELVAL(2) 634 : MLREEL = MELVA1.IELCHE(JJ,iem1) 635 : if(MLREEL.ne.0) then 636 : PSIX = 0.d0 637 : SEGACT,MLREEL 638 : do iii0 = 1,nbnn1 639 : PSIX = PSIX + (SHPP1(1,iii0) * PROG(iii0)) 640 : enddo 641 : SEGDES,MLREEL 642 : RX05= SQRT(ABS(PSIX)) 643 : c write(ioimp,*) 'RX05 = ', RX05 644 : do II=1,nbnn2 645 : RINT2(II,JJ1) = RINT2(II,JJ1) + RX05*REL1(II,JJ) 646 : enddo 647 : endif 648 : ENDIF 649 : 650 : 651 : c B2X 652 : IF(nbenrj.ge.3) THEN 653 : MELVA1 = MCHAM1.IELVAL(3) 654 : MLREEL = MELVA1.IELCHE(JJ,iem1) 655 : if(MLREEL.ne.0) then 656 : PSIX = 0.d0 657 : SEGACT,MLREEL 658 : do iii0 = 1,nbnn1 659 : PSIX = PSIX + (SHPP1(1,iii0) * PROG(iii0)) 660 : enddo 661 : SEGDES,MLREEL 662 : RX05= SQRT(ABS(PSIX)) 663 : do II=1,nbnn2 664 : RINT3(II,JJ1) = RINT3(II,JJ1) + RX05*REL1(II,JJ) 665 : enddo 666 : endif 667 : ENDIF 668 : 669 : 4205 CONTINUE 670 : 671 : if(iimpi.ge.4) then 672 : DO II=1,nbnn2 673 : write(ioimp,*) 'RINT1(',II,',:)=',(RINT1(II,iou),iou=1,nbnno2) 674 : ENDDO 675 : endif 676 : 677 : 678 : c on a trouvé le iem1 élément et on a fait le travail : 679 : c on passe au point de gauss suivant 680 : goto 1132 681 : 682 : 1131 CONTINUE 683 : c---------- fin de la Boucle sur les elements de structure ---------------------- 684 : if(iimpi.ge.2) write(ioimp,*) 685 : 1 'Fin de la boucle sur les elements de structure' 686 : if(izok.eq.0) write(ioimp,*) 'Attention le pt de Gauss ' 687 : & ,igau2,'de l element ',iem2,' est hors support' 688 : 689 : 690 : 1132 CONTINUE 691 : c======= fin de la Boucle sur les pt de G de fissure ============================ 692 : 693 : c --- Remplissage du MELEME et DU XMATRI --- 694 : DO ityp1=1,2 695 : 696 : c on commence par le MELEME des ddls Obligatoire (U) 697 : izone = 1 698 : c on poursuit avec le MELEME des ddls Sauts (H ou HB1 ou HB1B2) 699 : if(ityp1.eq.2) izone=nbenrj+1 700 : ity = ITAMYT(nbnno2,izone) 701 : 702 : do iidim=1,idim 703 : ityty = ((ity-1) * idim) + iidim 704 : if(iimpi.ge.4) write(ioimp,*) ity,nbnno2,izone,' -> ',ityty 705 : 706 : c --- Remplissage du MELEME --- 707 : MELEME = IRIGEL(1,ityty) 708 : NBNN = NUM(/1) 709 : NBELEM = NUM(/2) + 1 710 : if(iimpi.ge.4) then 711 : write(ioimp,*) 'on remplit le type de matrice ',ityty, 712 : 1 ' NBNN,NBELEM=',NBNN,NBELEM 713 : endif 714 : segadj,MELEME 715 : c element de fissure 716 : do ino2=1,nbnn2 717 : NUM(ino2,NBELEM) = IPT2.NUM(ino2,iem2) 718 : c NUM(ino2+nbnn2,NBELEM) = IPT2.NUM(ino2,iem2) 719 : enddo 720 : c element de structure 721 : do iino2=1,nbnno2 722 : c NUM((2*nbnn2)+iino2,NBELEM) = PRELEM(iino2) 723 : NUM(nbnn2+iino2,NBELEM) = PRELEM(iino2) 724 : enddo 725 : 726 : c --- Remplissage du XMATRI.RE --- 727 : XMATRI = IRIGEL(4,ityty) 728 : NLIGRD = RE(/1) 729 : NLIGRP = RE(/2) 730 : NELRIG = RE(/3) + 1 731 : segadj,XMATRI 732 : 733 : DO II=1,nbnn2 734 : 735 : c -Ktt 736 : c on choisit de ne garder que la partie diagonale 737 : JJ=II 738 : XMATRI.RE(II,JJ,NBELEM) = RINT(II,JJ)*XISTAB 739 : 740 : c -Kwt 741 : DO JJ=1,nbnn2 742 : XMATRI.RE(II,JJ+nbnn2,NBELEM) = RINT(II,JJ) 743 : XMATRI.RE(JJ+nbnn2,II,NBELEM) = RINT(II,JJ) 744 : c write(ioimp,*) 'XMATRI.RE(',II,JJ+nbnn2,NBELEM,') =', RINT(II,JJ) 745 : ENDDO 746 : 747 : C -Kww 748 : c on choisit de ne garder que la partie diagonale 749 : JJ=II 750 : XMATRI.RE(nbnn2+II,nbnn2+JJ,NBELEM) = RINT(II,JJ)*RLATIN 751 : 752 : C -Kut (on a deja fait le travail de positionnement avant) 753 : if(ityp1.eq.1) then 754 : c UX 755 : JJdec = 2*nbnn2 756 : DO JJ=1,nbnno2 757 : XMATRI.RE(II,JJdec+JJ,NBELEM) = -1.*RINT0(II,JJ) 758 : XMATRI.RE(JJdec+JJ,II,NBELEM) = -1.*RINT0(II,JJ) 759 : c write (6,* ) 'XMATRI.RE(',II,JJdec+JJ,NBELEM,') =', XMA 760 : c 1TRI.RE(II,JJdec+JJ,NBELEM) 761 : ENDDO 762 : else 763 : c AX 764 : JJdec = 2*nbnn2 765 : DO JJ=1,nbnno2 766 : XMATRI.RE(II,JJdec+JJ,NBELEM) = -1.*RINT1(II,JJ) 767 : XMATRI.RE(JJdec+JJ,II,NBELEM) = -1.*RINT1(II,JJ) 768 : c write (6,* ) 'XMATRI.RE(',II,JJdec+JJ,NBELEM,') =', XMA 769 : c 1TRI.RE(II,JJdec+JJ,NBELEM) 770 : ENDDO 771 : 772 : c B1X 773 : IF(nbenrj.ge.2) THEN 774 : JJdec=JJdec+nbnno2 775 : DO JJ=1,nbnno2 776 : XMATRI.RE(II,JJdec+JJ,NBELEM) = -1.*RINT2(II,JJ) 777 : XMATRI.RE(JJdec+JJ,II,NBELEM) = -1.*RINT2(II,JJ) 778 : 779 : c Cas ou les inconnues aux noeuds ne sont pas présentes 780 : c mis à 0 de ces inconnues en mettant un 1 sur la diagonale 781 : 782 : c Début de la boucle pour tester les noeuds de structure d'ancien point de gauss 783 : cbp DO 5002 iigau2 = 1,igau2 784 : DO 5002 iigau2 = 1,ngau2 785 : iiem1 = elemPG(iigau2) 786 : DO JJ1=1,nbnn1 787 : iino1=IPT1.NUM(JJ1,iiem1) 788 : if(PRELEM(JJ).eq.iino1) then 789 : MELVA1 = MCHAM1.IELVAL(2) 790 : MLREEL = MELVA1.IELCHE(JJ1,iiem1) 791 : if(MLREEL.ne.0) goto 5002 792 : XMATRI.RE(JJdec+JJ,JJdec+JJ,NBELEM) = 1 793 : c write(ioimp,*) 'XMATRI.RE(',JJd ec+JJ,JJdec+JJ,NBEL 794 : c 1EM,') =', XMATRI.RE(JJdec+JJ,JJdec+JJ,NBELEM) 795 : endif 796 : ENDDO 797 : 5002 CONTINUE 798 : ENDDO 799 : 800 : 801 : 802 : ENDIF 803 : c B2X 804 : IF(nbenrj.ge.3) THEN 805 : JJdec=JJdec+nbnno2 806 : DO JJ=1,nbnno2 807 : XMATRI.RE(II,JJdec+JJ,NBELEM) = -1.*RINT3(II,JJ) 808 : XMATRI.RE(JJdec+JJ,II,NBELEM) = -1.*RINT3(II,JJ) 809 : ENDDO 810 : ENDIF 811 : endif 812 : 813 : ENDDO 814 : 815 : enddo 816 : 817 : 818 : ENDDO 819 : 820 : 821 : 822 : 1100 CONTINUE 823 : c==== fin de la Boucle sur les elements de fissure ============================== 824 : c---------------------------------------------------------------------- 825 : 826 : 827 : c------------------------------------ 828 : c AJUSTEMENT AVANT DE QUITTER 829 : c------------------------------------ 830 : 831 : if(iimpi.ge.2) write(ioimp,*) 'AJUSTEMENT AVANT DE QUITTER' 832 : 833 : 834 : C BOUCLE SUR LES SOUS RIGIDITÉS 835 : ityok = 0 836 : DO 2000 ityty=1,(idim*NTYMAT) 837 : 838 : MELEME = IRIGEL(1,ityty) 839 : DESCR = IRIGEL(3,ityty) 840 : XMATRI = IRIGEL(4,ityty) 841 : NBELEM = NUM(/2) 842 : if(NBELEM.ne.0) then 843 : ityok = ityok + 1 844 : IRIGEL(1,ityok) = MELEME 845 : IRIGEL(3,ityok) = DESCR 846 : IRIGEL(4,ityok) = XMATRI 847 : segdes,MELEME,XMATRI 848 : else 849 : segsup,MELEME,DESCR,XMATRI 850 : endif 851 : 2000 continue 852 : NRIGEL = ityok 853 : segadj,MRIGID 854 : 855 : c------------------------------- 856 : c MENAGE AVANT DE QUITTER 857 : c------------------------------- 858 : 859 : segsup,MTRAV 860 : segdes,MINTE1,MINTE2 861 : segdes,IPT1,IPT2 862 : segdes,MCHAM1 863 : 864 : 865 : c------------------------------- 866 : c ON REMPLACE LES FX PAR DES LX 867 : c------------------------------- 868 : IF(INOMU.eq.0) THEN 869 : JGN=4 870 : JGM=2*idim 871 : segini,MLMOT1 872 : do iidim=1,idim 873 : MLMOT1.MOTS(iidim) = DUAOBL(iidim) 874 : MLMOT1.MOTS(iidim+idim) = DUAFAC(iidim) 875 : enddo 876 : ILMOT1 = MLMOT1 877 : IPRIG1 = MRIGID 878 : if(iimpi.ge.2) write(ioimp,*)' APPEL A FX2LX ' 878 : CALL FX2LX(IPRIG1,ILMOT1,IPRIG2) 879 : MRIGID = IPRIG2 880 : SEGSUP,MLMOT1 881 : ENDIF 882 : 883 : c write(ioimp,*) 'desactivation du MRIGID',MRIGID 884 : SEGDES,MRIGID 885 : 886 : if(iimpi.ge.3) write(ioimp,*) 'ecriture du MRIGID',MRIGID 886 : CALL ECROBJ('RIGIDITE',MRIGID) 887 : 888 : RETURN 889 : 890 : END 891 : 892 : 893 : 894 : 895 : 896 : 897 : 898 :
© Cast3M 2003 - All rights reserved.
Disclaimer