1 : C TRIELE SOURCE BP208322 12/12/20 21:15:26 7632 2 : C======================================================= 3 : C CREATION : bp, le 18.10.2005 4 : C MODIFS : -modif de la projection : bp, le 07/02/2006 5 : c -refonte de l'operateur pour integration dans Cast3M (aout2008) 6 : c -mise à 0. de Phi et Psi moins aleatoire (septembre 2009) 7 : c -syntaxe 2 (bp, juin 2010) 8 : c -ajout option 'DESENRICHISSEMENT' (bp, avril 2012) 9 : c 10 : C======================================================= 11 : C FONCTION : 12 : C TRIE LES ELEMENTS SELON LES ISOZEROS DES CHPOINTS LEVEL SET 13 : C 14 : c 15 : C======================================================= 16 : c SYNTAXE 1 : 17 : c chelX_n+1 | | = TRIE mod_n psi_n+1 phi_n+1 | | ; 18 : c | | | 'SAUT' | 19 : c | relaX | | 'DESENR' | 20 : c 21 : c (mod_n est actualisé --> mod_n+1) 22 : C 23 : C E: phi_n+1 = MCHPOI distance au plan de fissure a linstant n+1 24 : C E: psi_n+1 = MCHPOI distance au front de fissure a linstant n+1 25 : C E: mod_n = MMODEL a linstant n 26 : C (contient eventuellement une zone avec des element X-FEM 27 : C et IVAMOD pointe alors vers chelX_n) 28 : C E: 'SAUT' = mot cle optionnel permettant de limiter le niveau 29 : C d'enrichissement a H. La pointe doit etre sur un bord 30 : C d'element 31 : C E: 'DESENR' = mot cle optionnel permettant de 'desenrichir' le modele 32 : C S: mod_n+1 = MMODEL a linstant n+1 (mise a jour de IVAMOD) 33 : C S: chelX_n+1 = MCHELM d'enrichissement a linstant n+1 34 : C (construit d'apres les iso-zero de phi et psi) 35 : C S: relaX : relation associée au desenrichissement 36 : c (mise a zero des 'pointes devenues obsolètes') 37 : C 38 : C======================================================= 39 : c SYNTAXE 2 : 40 : c TRIE mod_? chelX_n+1 41 : c mod_n+1 42 : c 43 : C E: mod_? = MMODEL contenant ou pas des enrichissements 44 : C E: chelX_n+1 = MCHELM d'enrichissement a l instant n+1 45 : C S: mod_n+1 = MMODEL pointant vers chelX_n+1 (inscrit dans IVAMOD) 46 : C 47 : C 48 : C======================================================= 49 : C 50 : SUBROUTINE TRIELE 51 : c 52 : C======================================================= 53 : c 54 : IMPLICIT INTEGER(I-N) 55 : IMPLICIT REAL*8 (A-H,O-Z) 56 : 57 : LOGICAL LOG1 58 : CHARACTER*4 nomc1,nomc2,char1 59 : PARAMETER (NMOT0=2) 60 : CHARACTER*4 MOTCLE(NMOT0) 61 : DATA MOTCLE/'SAUT','DESE'/ 62 : 63 : PARAMETER (NDDL1=12) 64 : CHARACTER*4 MOTPR1(NDDL1),MOTDU1(NDDL1) 65 : DATA MOTPR1/'B1X ','B1Y ','B1Z ','C1X ','C1Y ','C1Z ', 66 : >'D1X ','D1Y ','D1Z ','E1X ','E1Y ','E1Z '/ 67 : DATA MOTDU1/'FB1X','FB1Y','FB1Z','FC1X','FC1Y','FC1Z', 68 : >'FD1X','FD1Y','FD1Z','FE1X','FE1Y','FE1Z'/ 69 : CHARACTER*4 MOTPR2(NDDL1),MOTDU2(NDDL1) 70 : DATA MOTPR2/'B2X ','B2Y ','B2Z ','C2X ','C2Y ','C2Z ', 71 : >'D2X ','D2Y ','D2Z ','E2X ','E2Y ','E2Z '/ 72 : DATA MOTDU2/'FB2X','FB2Y','FB2Z','FC2X','FC2Y','FC2Z', 73 : >'FD2X','FD2Y','FD2Z','FE2X','FE2Y','FE2Z'/ 74 : 75 : 76 : C SEGMENTS INCLUDE 77 : -INC CCOPTIO 78 : -INC SMCOORD 79 : -INC SMELEME 80 : -INC SMCHPOI 81 : -INC SMMODEL 82 : -INC SMCHAML 83 : -INC SMLREEL 84 : -INC SMINTE 85 : -INC SMRIGID 86 : C 87 : C Segment (type LISTENTI) contenant les informations sur un element 88 : SEGMENT INFO 89 : INTEGER INFELL(JG) 90 : ENDSEGMENT 91 : c 92 : SEGMENT ITYPND(NI) 93 : SEGMENT NDPT1(NI) 94 : SEGMENT NDPT2(NI) 95 : SEGMENT NUMND(NI) 96 : SEGMENT KKGAU(NI) 97 : SEGMENT IDEJVU(NBPTI) 98 : 99 : c segment pour tableau de travail de reel 100 : SEGMENT XPSI(NI) 101 : SEGMENT XPHI(NI) 102 : 103 : c autrespointeurs 104 : POINTEUR MCHEX1.MCHELM,MCHEX2.MCHELM 105 : POINTEUR IPTX1.MELEME 106 : POINTEUR MINT1.MINTE,MINT2.MINTE,MINT0.MINTE 107 : 108 : 109 : 110 : C======================================================= 111 : c 112 : IDIM1 = IDIM + 1 113 : segact,MCOORD*mod 114 : NBPTS = XCOOR(/1)/(IDIM1) 115 : 116 : C======================================================= 117 : c WRITE(*,*) 'LECTURE' 118 : C======================================================= 119 : 120 : C++++ le modele a l'instant n 120 : CALL LIROBJ('MMODEL',MMODE1,1,IRETOU) 121 : if (IERR.ne.0) then 122 : write(*,*) 'Pb de lecture du modele' 123 : return 124 : endif 125 : SEGACT,MMODE1 126 : NBSMOD = MMODE1.KMODEL(/1) 127 : 128 : c++++ L'option pour le niveau d'enrichissement max souhaité 129 : *as 2010_01_20 : selon l'option le niveau d'enrichissement est limité à 1 ou non: 130 : lsaut=100 130 : c CALL LIRCHA(char1,0,iretou) 131 : c if (iretou.ne.0) then 132 : c if (char1.ne.'SAUT') then 133 : c MOTERR(1:4)=char1 133 : c CALL ERREUR(7) 134 : c RETURN 135 : c else 136 : c lsaut=1000 137 : c endif 138 : c endif 139 : *fin as 2010_01_20 139 : CALL LIRMOT(MOTCLE,NMOT0,IMOT0,0) 140 : if(IMOT0.eq.1) lsaut=1000 141 : 142 : c++++ lecture d un MCHEX1 d'enrichissement a brancher 142 : CALL LIROBJ('MCHAML',MCHEX1,0,ISYNTAX) 143 : 144 : if (ISYNTAX.ne.0) then 145 : 146 : segact,MCHEX1 147 : if ((MCHEX1.TITCHE).ne.'ENRICHIS') then 148 : write(*,*) 'Le MCHELM fourni doit avoir ENRICHIS pour titre' 149 : segdes,MCHEX1 150 : return 151 : endif 152 : segdes,MCHEX1 153 : 154 : else 155 : 156 : C++++ lecture des Level Set 156 : CALL LIROBJ('CHPOINT',MCHPO1,1,IRETOU) 156 : CALL LIROBJ('CHPOINT',MCHPO2,1,IRETOU) 157 : if (IERR.ne.0) then 158 : write(*,*) 'Besoin des chpoints Level Sets' 159 : return 160 : endif 161 : C qq verif de base sur les chpoints 162 : segact MCHPO1,MCHPO2 163 : nsoup1 = MCHPO1.IPCHP(/1) 164 : nsoup2 = MCHPO2.IPCHP(/1) 165 : if ((nsoup1.ne.1).or.(nsoup2.ne.1)) then 166 : write(*,*) 'Chpoint avec plus d 1 zone: Non traite' 167 : return 168 : endif 169 : MSOUP1 = MCHPO1.IPCHP(1) 170 : MSOUP2 = MCHPO2.IPCHP(1) 171 : segdes MCHPO1,MCHPO2 172 : C composante des chpoints: 1=PSI 2=PHI 173 : segact MSOUP1,MSOUP2 174 : nc1 = MSOUP1.NOCOMP(/2) 175 : nc2 = MSOUP2.NOCOMP(/2) 176 : if ((nc1.ne.1).or.(nc2.ne.1)) then 177 : write(*,*) 'Chpoint avec plus d 1 composante: Non traite' 178 : return 179 : endif 180 : nomc1 = MSOUP1.NOCOMP(1) 181 : nomc2 = MSOUP2.NOCOMP(1) 182 : if (nomc1.ne.'PSI ') then 183 : write(*,*) 'Chpoint 1 doit avoir pour composante PSI ' 184 : return 185 : endif 186 : if (nomc2.ne.'PHI ') then 187 : write(*,*) 'Chpoint 2 doit avoir pour composante PHI ' 188 : return 189 : endif 190 : C recuperation des maillages et valeurs 191 : IPT1 = MSOUP1.IGEOC 192 : IPT2 = MSOUP2.IGEOC 193 : c if (IPT1.ne.IPT2) then 194 : c write(*,*) 'Chpoints doivent avoir le meme support maillage' 195 : c return 196 : c bp (07/07/2011) : trop contraignant on introduit ndpt2 et on 197 : c va faire des test + tard 198 : c endif 199 : MPOVA1 = MSOUP1.IPOVAL 200 : MPOVA2 = MSOUP2.IPOVAL 201 : segdes MSOUP1,MSOUP2 202 : C on crée MPOVA3 et MPOVA4 = versions "simplifiées" de MPOVA1 et MPOVA2 203 : SEGINI,MPOVA3=MPOVA1 204 : SEGINI,MPOVA4=MPOVA2 205 : 206 : endif 207 : 208 : 209 : C======================================================= 210 : C BOUCLE SUR LES SOUS MODELES 211 : C======================================================= 212 : DO 100 ISMOD1 = 1,NBSMOD 213 : c 214 : CcccccC ouverture du sous modele 215 : IMODEL = MMODE1.KMODEL(ISMOD1) 216 : segact IMODEL*mod 217 : 218 : c si EF XFEM, on travaille 219 : * if(NEFMOD.ne.263) goto 100 220 : *as 2010_01_22 : ajout de l'element 3D 221 : c if (NEFMOD.eq.263.or.NEFMOD.eq.264) then 222 : if(NEFMOD.ne.263.and.NEFMOD.ne.264) goto 100 223 : 224 : NOBMO1 = IVAMOD(/1) 225 : 226 : ccccccC branchement d 1 eventuel MCHEX1 d'enrichissement fourni 227 : if (ISYNTAX.ne.0) then 228 : NOBMO1 = IVAMOD(/1) 229 : c NOBMOD = NOBMO1 + 1 230 : if(NOBMO1.ne.0) write(6,*)'Attention: Les objets contenus', 231 : & ' dans le IVAMOD du modele seront oubliés !!!' 232 : NOBMOD=1 233 : cbp: attention: on ecrase tout pour l instant 234 : MN3 = INFMOD(/1) 235 : NFOR = FORMOD(/2) 236 : NMAT = MATMOD(/2) 237 : segadj,IMODEL 238 : TYMODE(NOBMOD)='MCHAML ' 239 : IVAMOD(NOBMOD) = MCHEX1 240 : segdes,IMODEL 241 : goto 100 242 : endif 243 : 244 : 245 : ccccccC recup du MCHEX1 d'enrichissement pre-existant 246 : if (NOBMO1.ne.0) then 247 : do iobmo1=1,NOBMO1 248 : if ((TYMODE(iobmo1)).eq.'MCHAML ') then 249 : MCHEX1 = IVAMOD(iobmo1) 250 : segact,MCHEX1 251 : if ((MCHEX1.TITCHE).eq.'ENRICHIS') then 252 : isouX=1 253 : MCHAM1 = MCHEX1.ICHAML(isouX) 254 : cbp: on suppose donc l'exitence d'autant de MCHEX1 (avec 1 seul MCHAM1) 255 : c que de zones elementaires, au lieu de 1 MCHEX1 avec autant de MCHAM1 256 : c que de zones elementaires... il faudrait corriger ca tres vite... 257 : segact,MCHAM1 258 : goto 101 259 : endif 260 : segdes,MCHEX1 261 : endif 262 : enddo 263 : endif 264 : 265 : ccccccC creation d un MCHEX1 d'enrichissement vierge 266 : write(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' 267 : write(*,*) 'Le modele ne contient pas le MCHELM ENRICHIS' 268 : write(*,*) '=> On initialise a 0 le MCHELM ENRICHIS' 269 : write(*,*) '!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' 270 : L1 = 8 271 : N1 = 1 272 : N3 = 6 273 : SEGINI,MCHEX1 274 : MCHEX1.TITCHE = 'ENRICHIS' 275 : MCHEX1.IFOCHE = IFOUR 276 : MCHEX1.IMACHE(1)= IMAMOD 277 : MCHEX1.CONCHE(1)= IMODEL.CONMOD 278 : MCHEX1.INFCHE(1,2) = 0 279 : MCHEX1.INFCHE(1,3) = NIFOUR 280 : MCHEX1.INFCHE(1,6) = 1 281 : *as 2010_01_20: 282 : * N2 = 2 283 : N2 = 1 284 : SEGINI,MCHAM1 285 : MCHEX1.ICHAML(1) = MCHAM1 286 : MELEME=IMAMOD 287 : SEGACT,MELEME 288 : N1PTEL = 0 289 : N1EL = 0 290 : N2PTEL = NUM(/1) 291 : N2EL = NUM(/2) 292 : SEGDES,MELEME 293 : MCHAM1.NOMCHE(1)= 'H ' 294 : *as 2010_01_20: 295 : * MCHAM1.NOMCHE(2)= 'F1 ' 296 : DO I2=1,N2 297 : MCHAM1.TYPCHE(I2)= 'POINTEURLISTREEL' 298 : SEGINI,MELVA1 299 : MCHAM1.IELVAL(I2) = MELVA1 300 : SEGDES,MELVA1 301 : ENDDO 302 : c on ajoute une petite place dans IVAMOD 303 : NOBMOD=NOBMO1+1 304 : MN3=INFMOD(/1) 305 : NFOR=FORMOD(/2) 306 : NMAT=MATMOD(/2) 307 : segadj,IMODEL 308 : iobmo1=NOBMOD 309 : 310 : ccccccC point de rassemblement une fois obtenu MCHEX1 et MCHAM1 311 : 101 continue 312 : c qq verif de base 313 : NCHAM1 = MCHEX1.ICHAML(/1) 314 : if (NCHAM1.ne.1) then 315 : write(*,*) 'MCHELM doit avoir 1 seule zone' 316 : return 317 : endif 318 : 319 : c recup du maillage et du MCHAML associé a cette zone 320 : MELEME = MCHEX1.IMACHE(1) 321 : if (MELEME.ne.IMAMOD) then 322 : write(*,*) 'MCHAML et ss-Model doivent avoir le meme support' 323 : return 324 : endif 325 : segact,MELEME 326 : MCHAM1 = MCHEX1.ICHAML(1) 327 : 328 : c recup de l enrichissement précédent (ITIP1) et calcul suivant (ITIP2) 329 : ITIP1=MCHEX1.INFCHE(1,2) 330 : if(ITIP1.lt.0.or.ITIP1.gt.2) then 331 : write(ioimp,*) 'Pb dans MCHEX1.INFCHE(1,2)' 331 : call erreur(21) 332 : endif 333 : if(IMOT0.eq.1) then 334 : ITIP2=0 335 : else 336 : if(ITIP1.eq.1) then 337 : ITIP2=2 338 : else 339 : c cas ou ITIP1 = 0 (initial) ou 2 340 : ITIP2=1 341 : endif 342 : endif 343 : 344 : 345 : C======================================================= 346 : c write(*,*) 'CREATION D UN JOLI MCHEX2 identique au MCHEX1' 347 : c et eventuellement d une matrice de relation pour désenrichir 348 : C======================================================= 349 : 350 : SEGINI,MCHEX2=MCHEX1 351 : IVAMOD(iobmo1)=MCHEX2 352 : TYMODE(iobmo1)='MCHAML ' 353 : MELE = NEFMOD 354 : c write(*,*) 'on branche MCHEX2=',MCHEX2,' dans IVAMOD(',iobmo1 355 : SEGINI,MCHAM2=MCHAM1 356 : MCHEX2.ICHAML(1)=MCHAM2 357 : MCHEX2.INFCHE(1,2)=ITIP2 358 : segdes,MCHEX1,MCHEX2 359 : N2ENR = MCHAM2.IELVAL(/1) 360 : do i2enr=1,N2ENR 361 : MELVA1 = MCHAM1.IELVAL(i2enr) 362 : SEGINI,MELVA2=MELVA1 363 : MCHAM2.IELVAL(i2enr) = MELVA2 364 : SEGDES,MELVA2 365 : enddo 366 : 367 : c --- DESENRICHISSEMENT --- 368 : IPRIG2=0 369 : if(IMOT0.eq.2) then 370 : 371 : c -On "oublie" les très vieux enrichissements devenus obsolètes 372 : c et qu on va réutiliser 373 : if(MCHAM2.IELVAL(/1).ge.ITIP2+1) then 374 : MELVA2=MCHAM2.IELVAL(ITIP2+1) 375 : segact,MELVA2*mod 376 : do j=1,MELVA2.IELCHE(/2) 377 : do i=1,MELVA2.IELCHE(/1) 378 : MELVA2.IELCHE(i,j)=0 379 : enddo 380 : enddo 381 : segdes,MELVA2 382 : endif 383 : 384 : c -On met a 0 (avec une relation) les vieux enrichissements 385 : if(MCHAM2.IELVAL(/1).ge.ITIP1+1) then 386 : 387 : c On commence par creer un maillage prototype de la rigidite (IPT6) 388 : c en reperant les noeuds a traiter via idejvu 389 : NBPTI=NBPTS 390 : segini,idejvu 391 : NBSOUS=0 392 : NBREF=0 393 : NBNN=2 394 : NBELEM=1000 395 : segini,IPT6 396 : IPT6.ITYPEL=22 397 : c on noublie pas de créer des noeuds pour les LX 398 : NBPTS0=NBPTS 399 : iilx=NBPTS0 400 : NBPTS=NBPTS+NBELEM 401 : segadj,MCOORD 402 : c on travaille sur le melva2 de la pointe précédente 403 : MELVA2=MCHAM2.IELVAL(ITIP1+1) 404 : segact,MELVA2*mod 405 : ii2=0 406 : do 120 j=1,MELVA2.IELCHE(/2) 407 : do 120 i=1,MELVA2.IELCHE(/1) 408 : i120 = NUM(i,j) 409 : c rem: commme l enrichissement est nodal et que l on travaille sur 1 seul 410 : c enrichissement a annuler on se permet les lignes ci dessous 411 : if(idejvu(i120).gt.0) goto 120 412 : idejvu(i120)=1 413 : if (MELVA2.IELCHE(i,j).ne.0) then 414 : ii2=ii2+1 415 : if(ii2.gt.NBELEM) then 416 : NBELEM=NBELEM + 1000 417 : NBPTS=NBPTS + 1000 418 : segadj,IPT6,MCOORD 419 : endif 420 : c on crée un noeud pour le multiplicateur de lagrange LX 421 : c iilx = NBPTS0+ii2 422 : iilx = iilx+1 423 : XCOOR((iilx-1)*(IDIM1)+1) = XCOOR((i120-1)*(IDIM1)+1) 424 : XCOOR((iilx-1)*(IDIM1)+2) = XCOOR((i120-1)*(IDIM1)+2) 425 : if(idim.eq.3) 426 : & XCOOR((iilx-1)*(IDIM1)+3) = XCOOR((i120-1)*(IDIM1)+3) 427 : XCOOR(iilx*(IDIM1)) = XCOOR(i120*(IDIM1)) 428 : IPT6.NUM(1,ii2)=iilx 429 : IPT6.NUM(2,ii2)=i120 430 : endif 431 : 120 continue 432 : c on ajuste IPT6 et on supprime 433 : NBELEM=ii2 434 : segadj,IPT6 435 : segdes,MELVA2 436 : segsup,idejvu 437 : 438 : c on finit d ajuster MCOORD 439 : NRIGEL=idim*4 440 : NBPTS = NBPTS0+(ii2*NRIGEL) 441 : segadj,MCOORD 442 : 443 : c On crée la relation qui permet de désenrichir 444 : segini,MRIGID 445 : IPRIG2=MRIGID 446 : IFORIG = IFOMOD 447 : MTYMAT ='RIGIDITE' 448 : NLIGRD=2 449 : NLIGRP=2 450 : NELRIG=ii2 451 : jb=0 452 : jdim=0 453 : do iriri=1,NRIGEL 454 : COERIG(iriri)=1.d0 455 : segini,DESCR,XMATRI 456 : c IRIGEL(1,iriri)=IPT6 457 : c On duplique le maillage prototype de la rigidite 458 : if(iriri.eq.1) then 459 : IPT7=IPT6 460 : else 461 : segini,IPT7=IPT6 462 : c on crée un noeud pour le multiplicateur de lagrange LX 463 : c pour chaque element de ipt7 464 : do ii2=1,NBELEM 465 : iilx = iilx+1 466 : c iilx6=IPT6.NUM(1,ii2) on prend ipt7 actif et pas ipt6 467 : iilx6=IPT7.NUM(1,ii2) 468 : XCOOR((iilx-1)*(IDIM1)+1) = XCOOR((iilx6-1)*(IDIM1)+1) 469 : XCOOR((iilx-1)*(IDIM1)+2) = XCOOR((iilx6-1)*(IDIM1)+2) 470 : if(idim.eq.3) 471 : & XCOOR((iilx-1)*(IDIM1)+3) = XCOOR((iilx6-1)*(IDIM1)+3) 472 : XCOOR(iilx*(IDIM1)) = XCOOR(iilx6*(IDIM1)) 473 : IPT7.NUM(1,ii2)=iilx 474 : enddo 475 : endif 476 : segdes,IPT7 477 : IRIGEL(1,iriri)=IPT7 478 : IRIGEL(2,iriri)=0 479 : IRIGEL(3,iriri)=DESCR 480 : IRIGEL(4,iriri)=XMATRI 481 : IRIGEL(5,iriri)=NIFOUR 482 : IRIGEL(6,iriri)=0 483 : IRIGEL(7,iriri)=0 484 : IRIGEL(8,iriri)=0 485 : c remplissage du DESCR 486 : LISINC(1)='LX' 487 : LISDUA(1)='FLX' 488 : jb=jb+1 489 : jdim=jdim+1 490 : if(ITIP1.eq.1) then 491 : LISINC(2)=MOTPR1(jb) 492 : LISDUA(2)=MOTDU1(jb) 493 : elseif(ITIP1.eq.2) then 494 : LISINC(2)=MOTPR2(jb) 495 : LISDUA(2)=MOTDU2(jb) 496 : endif 497 : c petit ajustement pour "sauter" les composantes en Z 498 : if(jdim.eq.idim) then 499 : jdim=0 500 : if(idim.eq.2) jb=jb+1 501 : endif 502 : do jinc=1,NLIGRP 503 : NOELEP(jinc)=jinc 504 : NOELED(jinc)=jinc 505 : enddo 506 : c remplissage de XMATRI 507 : do jel=1,NELRIG 508 : RE(1,1,jel)=0.d0 509 : RE(1,2,jel)=1.d0 510 : RE(2,1,jel)=1.d0 511 : RE(2,2,jel)=0.d0 512 : enddo 513 : segdes,DESCR,XMATRI 514 : enddo 515 : segdes,MRIGID 516 : 139 continue 517 : endif 518 : 519 : endif 520 : 521 : segdes,MCHAM1,MCHAM2 522 : 523 : 524 : C======================================================= 525 : c BOUCLE SUR LES POINTS DES CHPOINTS 526 : C======================================================= 527 : 528 : c -chpoint psi 529 : segact,IPT1 530 : NBNO1 = IPT1.NUM(/2) 531 : C on crée la table qui depuis num de noeud donne point dans chpoint 532 : NI = xcoor(/1)/(IDIM1) 533 : SEGINI,NDPT1 534 : cccccc Vrai Debut de boucle 535 : DO I1=1,NBNO1 536 : C on recupere le numero de noeud 537 : INUM1 = IPT1.NUM(1,I1) 538 : C on remplit la table 539 : NDPT1(INUM1) = I1 540 : ENDDO 541 : 542 : c -chpoint phi 543 : if (IPT1.eq.IPT2) then 544 : NDPT2=NDPT1 545 : segdes,IPT1 546 : else 547 : segact,IPT2 548 : NBNO2 = IPT2.NUM(/2) 549 : if (NBNO1.ne.NBNO2) then 550 : write(*,*) 'Les maillages support des chpoints ', 551 : & 'doivent avoir le meme nombre de points!' 552 : return 553 : endif 554 : SEGINI,NDPT2 555 : DO I2=1,NBNO2 556 : C on recupere le numero de noeud 557 : INUM2 = IPT2.NUM(1,I2) 558 : if (NDPT1(INUM2).eq.0) then 559 : write(*,*) 'Les maillages support des chpoints ', 560 : & 'doivent avoir les memes points!' 561 : return 562 : endif 563 : C on remplit la table 564 : NDPT2(INUM2) = I2 565 : ENDDO 566 : segdes,IPT1,IPT2 567 : endif 568 : C=======FIN DE BOUCLE SUR LES POINTS DES CHPOINTS======= 569 : 570 : 571 : C======================================================= 572 : c 1ere BOUCLE SUR LES ELEMENTS 573 : c POUR LA "SIMPLIFICATION" 574 : C======================================================= 575 : C Ouverture du maillage et des chpoints level set 576 : c SEGACT MELEME,MPOVA1,MPOVA2 577 : SEGACT MPOVA1,MPOVA2 578 : NBNN1 = NUM(/1) 579 : NBEL1 = NUM(/2) 580 : c on recupere les coord dans l element parents des noeuds 580 : call elquoi(MELE,0,1,IPINF,IMODEL) 581 : INFO = IPINF 582 : MINT0 = INFELL(11) 583 : segact,MINT0 584 : c on recupere les coord dans l element parents des points de gauss 584 : call elquoi(MELE,0,3,IPINF,IMODEL) 585 : INFO = IPINF 586 : MINT1 = INFELL(11) 587 : segact,MINT1 588 : c segment integration ef std correspondant -> inutile pour l instant 589 : c MINT2 = INFELL(12) 590 : c if(MINT2.ne.0) segact,MINT2 591 : segdes,IMODEL 592 : c recherche du num du pt de gauss KKGAU le plus proche du noeud inode1 593 : NGAU1 = INFELL(6) 594 : NI = NBNN1 595 : segini,KKGAU 596 : do inode1=1,NBNN1 597 : kgau1 = 1 598 : DXi1 = ((MINT1.QSIGAU(kgau1)) - (MINT0.QSIGAU(inode1)))**2 599 : & + ((MINT1.ETAGAU(kgau1)) - (MINT0.ETAGAU(inode1)))**2 600 : & + ((MINT1.DZEGAU(kgau1)) - (MINT0.DZEGAU(inode1)))**2 601 : do kgau2=2,NGAU1 602 : DXi2 = ((MINT1.QSIGAU(kgau2)) - (MINT0.QSIGAU(inode1)))**2 603 : & + ((MINT1.ETAGAU(kgau2)) - (MINT0.ETAGAU(inode1)))**2 604 : & + ((MINT1.DZEGAU(kgau2)) - (MINT0.DZEGAU(inode1)))**2 605 : if (DXi2.lt.DXi1) then 606 : DXi1 = DXi2 607 : kgau1 = kgau2 608 : endif 609 : enddo 610 : KKGAU(inode1) = kgau1 611 : enddo 612 : c write(6,*) 'KKGAU' 613 : c write(6,*) (KKGAU(iou),iou=1,4) 614 : C Initialisation des tableaux de sortie 615 : NI = xcoor(/1)/(IDIM1) 616 : SEGINI,ITYPND 617 : NI = NBEL1 618 : c SEGINI,TYPEF 619 : C Initialisation des tableaux de travail 620 : NI = NBNN1 621 : segini,NUMND,XPSI,XPHI 622 : 623 : 624 : c++++++ Debut de boucle sur les elements 625 : DO 200 J=1,NBEL1 626 : 627 : C++++++++ Mini-Boucle(S) sur les noeuds de l element 628 : DO K=1,NBNN1 629 : C numero INUM du Kieme noeud du Jieme element 630 : INUM = NUM(K,J) 631 : C numero du point correspondant des chpoints 632 : I1 = NDPT1(INUM) 633 : I2 = NDPT2(INUM) 634 : c on verifie que le chpoint y est bien defini sinon on saute tout! 635 : if(I1.eq.0) goto 200 636 : C valeurs de psi et phi au Kieme noeud 637 : XPSI(K) = MPOVA1.VPOCHA(I1,1) 638 : XPHI(K) = MPOVA2.VPOCHA(I2,1) 639 : ENDDO 640 : C++++++++ Fin de la Mini-Boucle sur les noeuds de l element 641 : 642 : c c distance pour permettre la simplification 643 : c c bp: a supprimer tres bientot car inutile et aleatoire ... 644 : c (bp: supprimé le 13./07/2011) 645 : c EPSTOL = 1.D-3 646 : c INUM1 = NUM(1,J) 647 : c INUM2 = NUM(2,J) 648 : c XOM1 = ( ( XCOOR( (INUM2-1)*(IDIM1) + 1 ) ) 649 : c & - ( XCOOR( (INUM1-1)*(IDIM1) + 1 ) ) )**2 650 : c XOM2 = ( ( XCOOR( (INUM2-1)*(IDIM1) + 2 ) ) 651 : c & - ( XCOOR( (INUM1-1)*(IDIM1) + 2 ) ) )**2 652 : c if (idim.eq.3) then 653 : c XOM3 = ( ( XCOOR( (INUM2-1)*(IDIM1) + 3 ) ) 654 : c & - ( XCOOR( (INUM1-1)*(IDIM1) + 3 ) ) )**2 655 : c 656 : c XTOL = sqrt(XOM1+XOM2+XOM3) 657 : c else 658 : c XTOL = sqrt(XOM1+XOM2) 659 : c endif 660 : c XTOL = EPSTOL * XTOL 661 : c if(XTOL.le.1.D-20) 662 : c & write(6,*) '!!! attention distance inter noeud faible =',XTOL 663 : 664 : C++++++++ Mini-Boucle(S) sur les noeuds de l element 665 : DO K=1,NBNN1 666 : C on calcule phi au pts de gauss le + proche du noeud considéré 667 : kgau1 = KKGAU(K) 668 : XPSIGAU = 0.d0 669 : XPHIGAU = 0.d0 670 : do ii1=1,NBNN1 671 : XPSIGAU = XPSIGAU 672 : & + ( (MINT1.SHPTOT(1,ii1,kgau1)) * (XPSI(ii1)) ) 673 : XPHIGAU = XPHIGAU 674 : & + ( (MINT1.SHPTOT(1,ii1,kgau1)) * (XPHI(ii1)) ) 675 : enddo 676 : C numero INUM du Kieme noeud du Jieme element 677 : INUM = NUM(K,J) 678 : C numero du point correspondant des chpoints 679 : I1 = NDPT1(INUM) 680 : I2 = NDPT2(INUM) 681 : c C simplification 682 : c IF (((XPSI(K))*XPSIGAU).LE.0.d0) THEN 683 : c c XPSI(K) = 0.d0 684 : c MPOVA3.VPOCHA(I1,1) = 0.d0 685 : c ELSEIF (ABS(XPSI(K)).LE.XTOL) THEN 686 : c XPSI(K) = 0.d0 687 : c MPOVA3.VPOCHA(I1,1) = 0.d0 688 : c ENDIF 689 : c IF (((XPHI(K))*XPHIGAU).LE.0.d0) THEN 690 : c c XPHI(K) = 0.d0 691 : c MPOVA4.VPOCHA(I2,1) = 0.d0 692 : c ELSEIF (ABS(XPHI(K)).LE.XTOL) THEN 693 : c c XPHI(K) = 0.d0 694 : c MPOVA4.VPOCHA(I2,1) = 0.d0 695 : c ENDIF 696 : C simplification (bp: 13./07/2011) 697 : IF(((XPSI(K))*XPSIGAU).LE.0.d0) 698 : & MPOVA3.VPOCHA(I1,1) = 0.d0 699 : IF(((XPHI(K))*XPHIGAU).LE.0.d0) 700 : & MPOVA4.VPOCHA(I2,1) = 0.d0 701 : ENDDO 702 : C++++++++ Fin de la Mini-Boucle sur les noeuds de l element 703 : 704 : 200 CONTINUE 705 : C=======FIN DE la 1ere BOUCLE SUR LES ELEMENTS========== 706 : 707 : 708 : 709 : C======================================================= 710 : c 2EME BOUCLE SUR LES ELEMENTS 711 : C POUR LE TRI DES NOEUDS ("TEST EN PAGAILLE") 712 : C======================================================= 713 : 714 : c++++++ Debut de boucle sur les elements 715 : DO 222 J=1,NBEL1 716 : 717 : C++++++++ Mini-Boucle(S) sur les noeuds de l element 718 : DO K=1,NBNN1 719 : C numero INUM du Kieme noeud du Jieme element 720 : INUM = NUM(K,J) 721 : C numero du point correspondant des chpoints 722 : I1 = NDPT1(INUM) 723 : I2 = NDPT2(INUM) 724 : c on verifie que le chpoint y est bien defini sinon on saute tout! 725 : if(I1.eq.0) goto 222 726 : C numero INUM du Kieme noeud de l'element en cours (le Jieme) 727 : NUMND(K) = INUM 728 : C valeurs de psi et phi au Kieme noeud 729 : XPSIK = MPOVA3.VPOCHA(I1,1) 730 : XPHIK = MPOVA4.VPOCHA(I1,1) 731 : XPSI(K) = XPSIK 732 : XPHI(K) = XPHIK 733 : C Calcul de valeurs utiles POUR LE TRI 734 : if (K.eq.1) then 735 : PSIMAX = XPSIK 736 : PSIMIN = XPSIK 737 : PSIMIA = ABS(XPSIK) 738 : PHIMAX = XPHIK 739 : PHIMIN = XPHIK 740 : PHIMIA = ABS(XPHIK) 741 : else 742 : PSIMAX = MAX(PSIMAX,XPSIK) 743 : PSIMIN = MIN(PSIMIN,XPSIK) 744 : PSIMIA = MIN(PSIMIA,ABS(XPSIK)) 745 : PHIMAX = MAX(PHIMAX,XPHIK) 746 : PHIMIN = MIN(PHIMIN,XPHIK) 747 : PHIMIA = MIN(PHIMIA,ABS(XPHIK)) 748 : endif 749 : 750 : ENDDO 751 : C++++++++ Fin de la Mini-Boucle sur les noeuds de l element 752 : c if (JWRITE.eq.1) then 753 : c write(6,*) 'EF ',J,' psi,phi_node apres simplif=' 754 : c write(6,*) (XPSI(iou),iou=1,4) 755 : c write(6,*) (XPHI(iou),iou=1,4) 756 : c endif 757 : PSIPRD = PSIMAX * PSIMIN 758 : PHIPRD = PHIMAX * PHIMIN 759 : c 760 : C 761 : C++++++++ Tests en pagaille pour trier 762 : 763 : c LOG1 = Tous les noeuds de lelement sont situés en arriere du front? 764 : LOG1 = (PSIMAX.LE.0.d0) 765 : c ...noeuds en arriere du front + element coupé par la fissure 766 : IF (LOG1.AND.(PHIPRD.LT.0.d0)) THEN 767 : C => tous les noeud doivent etre H-enrichi (=10) 768 : DO K=1,NBNN1 769 : INUM = NUMND(K) 770 : ITYPND(INUM) = MAX(10,ITYPND(INUM)) 771 : ENDDO 772 : C 773 : c ...noeuds en arriere du front + element tangent avec la fissure 774 : ELSEIF (LOG1.AND.(PHIMIA.EQ.0.d0)) THEN 775 : C => seuls les noeud tangents doivent etre H-enrichi (=10) 776 : DO K=1,NBNN1 777 : IF (XPHI(K).EQ.0.d0) THEN 778 : INUM = NUMND(K) 779 : ITYPND(INUM) = MAX(10,ITYPND(INUM)) 780 : ENDIF 781 : ENDDO 782 : C 783 : c Element contenant la pointe de fissure 784 : ELSEIF ((PSIPRD.LT.0.d0).AND.(PHIPRD.LT.0.d0)) THEN 785 : C => tous les noeud doivent etre F-enrichi (=100) 786 : DO K=1,NBNN1 787 : INUM = NUMND(K) 788 : ITYPND(INUM) = 100 789 : ENDDO 790 : C 791 : c Element dont un bord contient la pointe de fissure 792 : ELSEIF ((PSIPRD.LT.0.d0).AND.(PHIMIA.EQ.0.d0)) THEN 793 : C => tous les noeud doivent etre F-enrichi (=100) 794 : DO K=1,NBNN1 795 : INUM = NUMND(K) 796 : ITYPND(INUM) = 100 797 : ENDDO 798 : C 799 : c La pointe arrive sur un bord de l element 800 : ELSEIF ((PSIMIA.EQ.0.d0).AND.(PHIPRD.LT.0.d0)) THEN 801 : C => choix: seuls les noeud en avant du front sont F-enrichi (=100) 802 : DO K=1,NBNN1 803 : INUM = NUMND(K) 804 : if (XPSI(K).GE.0.d0) ITYPND(INUM) = 100 805 : ENDDO 806 : C 807 : c Element dont un coin contient la pointe de fissure 808 : ELSEIF ((PSIMIA.EQ.0.d0).AND.(PHIMIA.EQ.0.d0)) THEN 809 : C => choix: seuls les noeud en avant du front sont F-enrichi (=100) 810 : DO K=1,NBNN1 811 : INUM = NUMND(K) 812 : if (XPSI(K).GE.0.d0) ITYPND(INUM) = 100 813 : ENDDO 814 : C 815 : c Les autres elements ne sont pas directement coupés ni tangents avec la fissure 816 : ELSE 817 : C => on met a 1 pour dire qu on les a vu 818 : DO K=1,NBNN1 819 : INUM = NUMND(K) 820 : ITYPND(INUM) = MAX(1,ITYPND(INUM)) 821 : ENDDO 822 : C 823 : ENDIF 824 : C++++++++ Fin des Tests en pagaille pour trier 825 : 826 : c DO K=1,NBNN1 827 : c INUM = NUMND(K) 828 : c TYPEF.VINT(J) = TYPEF.VINT(J) + ITYPND(INUM) 829 : c if(j.eq.660) write(*,*) 'K',K,' typnd',ITYPND(INUM) 830 : c if(j.eq.661) write(*,*) 'K',K,' typnd',ITYPND(INUM) 831 : c ENDDO 832 : 833 : 222 CONTINUE 834 : C=======FIN DE la 2eme BOUCLE SUR LES ELEMENTS========== 835 : segsup,NUMND,XPSI,XPHI,KKGAU 836 : segdes,MINT0,MINT1 837 : 838 : 839 : 840 : C======================================================= 841 : c 3EME BOUCLE SUR LES ELEMENTS 842 : C POUR ATTRIBUER L'ENRICHISSEMENT 843 : C======================================================= 844 : 845 : segact,MCHAM2 846 : * on evite d'ouvrir les MELVA2 pour chaque element 847 : DO I2ENR=1,N2ENR 848 : MELVA2 = MCHAM2.IELVAL(I2ENR) 849 : if(MELVA2.ne.0) segact,MELVA2*mod 850 : ENDDO 851 : segact,MPOVA1,MPOVA2 852 : 853 : c++++++ Debut de boucle sur les elements 854 : DO 300 J=1,NBEL1 855 : 856 : c pour eviter de creer bcp de listreels on utilise l astuce : 857 : mphi =0 858 : mpsiph=0 859 : 860 : c++++++++ Boucle sur les noeuds de l element j 861 : DO 310 K=1,NBNN1 862 : 863 : C Kieme noeud du Jieme element a pour numero INUM 864 : INUM = NUM(K,J) 865 : 866 : c on saute les noeuds non enrichis 867 : if((ITYPND(INUM)).le.1) goto 310 868 : 869 : Ccccccccccc noeud H-enrichi 870 : if ((ITYPND(INUM)).eq.10) then 871 : MELVA2 = MCHAM2.IELVAL(1) 872 : c segact,MELVA2*mod 873 : MLREE1 = MELVA2.IELCHE(K,J) 874 : 875 : c on H-enrichit ssi ce noeud n etait pas deja H-enrichi 876 : if (MLREE1.eq.0) then 877 : if (mphi.ne.0) then 878 : MELVA2.IELCHE(K,J) = mphi 879 : else 880 : JG = NBNN1 881 : segini,MLREE2 882 : MELVA2.IELCHE(K,J) = MLREE2 883 : c on remplit la listreel des valeurs de PHI aux noeuds de l element 884 : do inode=1,NBNN1 885 : INUM2 = NUM(inode,J) 886 : I2 = NDPT2(INUM2) 887 : MLREE2.PROG(inode) = MPOVA2.VPOCHA(I2,1) 888 : enddo 889 : segdes,MLREE2 890 : mphi=MLREE2 891 : endif 892 : endif 893 : c segdes,MELVA2 894 : endif 895 : Ccccccccccc fin des noeud H-enrichi 896 : 897 : 898 : Ccccccccccc noeud F-enrichi 899 : * if((ITYPND(INUM)).eq.100) then 900 : *as 2010_01_20: lsaut =1000 si option 'SAUT' , =100 sinon 901 : if ((ITYPND(INUM)).eq.lsaut) then 902 : 903 : c si F1 pas prevu on le rajoute 904 : if (N2ENR.lt.2) then 905 : N2ENR = 2 906 : N2 = N2ENR 907 : segadj,MCHAM2 908 : MCHAM2.NOMCHE(N2ENR) = 'F1' 909 : MCHAM2.TYPCHE(N2ENR) = 'POINTEURLISTREEL' 910 : N1PTEL=0 911 : N1EL =0 912 : N2PTEL=NBNN1 913 : N2EL =NBEL1 914 : c segdes,MELVA2 915 : segini,MELVA2 916 : MCHAM2.IELVAL(N2ENR) = MELVA2 917 : else 918 : MELVA2 = MCHAM2.IELVAL(2) 919 : c segact,MELVA2*mod 920 : endif 921 : 922 : c - on F1-enrichit si ce noeud n etait pas deja F1-enrichi 923 : c OU si on est dans une stratégie de désenrichissement 924 : c et que c'est le tour de F1 925 : MLREE1 = MELVA2.IELCHE(K,J) 926 : c if (MLREE1.eq.0) then 927 : if ((MLREE1.eq.0.and.IMOT0.ne.2).OR 928 : & .(ITIP2.eq.1.and.IMOT0.eq.2)) then 929 : 930 : if (mpsiph.ne.0) then 931 : MELVA2.IELCHE(K,J) = mpsiph 932 : else 933 : NBLV7=2 934 : JG = NBLV7*NBNN1 935 : segini,MLREE2 936 : MELVA2.IELCHE(K,J) = MLREE2 937 : c segdes,MELVA2 938 : c on remplit la listreel des valeurs de PSI et PHI 939 : do inode=1,NBNN1 940 : INUM2 = NUM(inode,J) 941 : I1 = NDPT1(INUM2) 942 : I2 = NDPT2(INUM2) 943 : MLREE2.PROG(inode) = MPOVA1.VPOCHA(I1,1) 944 : MLREE2.PROG(NBNN1+inode) = MPOVA2.VPOCHA(I2,1) 945 : enddo 946 : segdes,MLREE2 947 : mpsiph=MLREE2 948 : endif 949 : 950 : c - si noeud deja F1-enrichi, alors on essaie de le F2-enrichir 951 : c OU si on est dans une stratégie de désenrichissement 952 : c et que c'est le tour de F2 953 : else 954 : 955 : c si F2 pas prevu on le rajoute 956 : if (N2ENR.lt.3) then 957 : N2ENR = 3 958 : N2 = N2ENR 959 : segadj,MCHAM2 960 : MCHAM2.NOMCHE(N2ENR) = 'F2' 961 : MCHAM2.TYPCHE(N2ENR) = 'POINTEURLISTREEL' 962 : N1PTEL=0 963 : N1EL =0 964 : N2PTEL=NBNN1 965 : N2EL =NBEL1 966 : c segdes,MELVA2 967 : segini,MELVA2 968 : MCHAM2.IELVAL(N2ENR) = MELVA2 969 : else 970 : MELVA2 = MCHAM2.IELVAL(3) 971 : segact,MELVA2*mod 972 : endif 973 : c +-- on F2-enrichit si ce noeud n etait pas deja F1-enrichi 974 : MLREE1 = MELVA2.IELCHE(K,J) 975 : if (MLREE1.eq.0) then 976 : if (mpsiph.ne.0) then 977 : MELVA2.IELCHE(K,J) = mpsiph 978 : else 979 : NBLV7=2 980 : JG = NBLV7*NBNN1 981 : segini,MLREE2 982 : MELVA2.IELCHE(K,J) = MLREE2 983 : c segdes,MELVA2 984 : c on remplit la listreel des valeurs de PSI et PHI 985 : do inode=1,NBNN1 986 : INUM2 = NUM(inode,J) 987 : I1 = NDPT1(INUM2) 988 : I2 = NDPT2(INUM2) 989 : MLREE2.PROG(inode) = MPOVA1.VPOCHA(I1,1) 990 : MLREE2.PROG(NBNN1+inode) = MPOVA2.VPOCHA(I2,1) 991 : enddo 992 : segdes,MLREE2 993 : mpsiph=MLREE2 994 : endif 995 : 996 : c +-- si on est la, c est que le noeud est deja F1 et F2 enrichi 997 : c => pour l instant on genere une erreur, 998 : c mais pour le futur ajouter la possiblite d avoir N2ENR >3 999 : else 1000 : write(*,*) 'Trop d enrichissement (> 2 pointes) au noeud',INUM 1001 : write(*,*) K,'ieme noeud du',J,'ieme element' 1002 : return 1003 : endif 1004 : c +-- fin de la distinction noeud deja F2-enrichi ou pas 1005 : 1006 : endif 1007 : c - fin de la distinction noeud deja F1-enrichi ou pas 1008 : 1009 : endif 1010 : Ccccccccccc fin du cas noeud F-enrichi 1011 : 1012 : 310 CONTINUE 1013 : c++++++++ Fin de Boucle sur les noeuds de l element j 1014 : 1015 : 1016 : 300 CONTINUE 1017 : C=======FIN DE la 3eme BOUCLE SUR LES ELEMENTS=========== 1018 : 1019 : 1020 : DO I2ENR=1,N2ENR 1021 : MELVA2 = MCHAM2.IELVAL(I2ENR) 1022 : if(MELVA2.ne.0) segdes,MELVA2 1023 : ENDDO 1024 : SEGDES,MELEME,MPOVA1,MPOVA2,MCHAM2 1025 : SEGSUP,MPOVA3,MPOVA4 1026 : if(NDPT1.ne.NDPT2) SEGSUP,NDPT2 1027 : SEGSUP,NDPT1,ITYPND 1028 : 1029 : *as 2010_01_22 :ajout de lelement 3D 1030 : c endif 1031 : 100 CONTINUE 1032 : C=====FIN DE BOUCLE SUR LES SOUS MODELES================ 1033 : c 1034 : SEGDES,MMODE1 1035 : C 1036 : 1037 : if(ISYNTAX.eq.0) then 1037 : if(IPRIG2.ne.0) CALL ECROBJ('RIGIDITE',IPRIG2) 1037 : CALL ECROBJ('MCHAML',MCHEX2) 1038 : endif 1039 : C 1040 : 1041 : RETURN 1042 : END 1043 : 1044 : C 1045 : 1046 : 1047 : 1048 : 1049 : 1050 : 1051 : 1052 : 1053 : 1054 : 1055 : 1056 :
© Cast3M 2003 - All rights reserved.
Disclaimer