1 : C XPOST2 SOURCE BP208322 13/11/13 21:15:27 7865 2 : SUBROUTINE XPOST2 3 : c 4 : c Construit 2 Chpoints avec des ddl "physiques" (UX UY) 5 : c sur le maillage de la fissure (en phi=0+ et phi=0-) 6 : c en RECOmbinant les ddl Xfem (UX AX B1X C1X ...) 7 : c -> Utile pour le tracé de l ouverture de fissure 8 : c 9 : C SYNTAXE : 10 : c chpo2 chpo3 = XFEM 'FISS' geofissure chpo1 MODEL_XFEM 11 : C 12 : c 13 : IMPLICIT INTEGER(I-N) 14 : IMPLICIT REAL*8 (A-H,O-Z) 15 : C 16 : PARAMETER (NBRMAX=5) 17 : 18 : C 19 : C SEGMENTS INCLUDE 20 : -INC CCOPTIO 21 : -INC CCREEL 22 : -INC SMCOORD 23 : -INC SMELEME 24 : -INC SMCHPOI 25 : -INC SMCHAML 26 : -INC SMMODEL 27 : -INC SMLREEL 28 : c 29 : POINTEUR MCHEX1.MCHELM 30 : C 31 : SEGMENT IDEJVU(NBPT2) 32 : c 33 : segment wrk4 34 : real*8 xel(3,nbnn),shpp(6,nbnn),qsi(3),XPT2(3) 35 : endsegment 36 : C 37 : C_____________________________________________________________ 38 : C INITIALISATION DES INCONNUES obligatoires et facultatives 39 : PARAMETER (NOBL=3,NFAC=27) 40 : CHARACTER*4 DDLOBL(NOBL),DDLFAC(NFAC),MODDL,MODDL2 41 : DATA DDLOBL/'UX ','UY ','UZ '/ 42 : DATA DDLFAC/'AX ','AY ','AZ ', 43 : >'B1X ','B1Y ','B1Z ', 44 : >'C1X ','C1Y ','C1Z ', 45 : >'D1X ','D1Y ','D1Z ', 46 : >'E1X ','E1Y ','E1Z ', 47 : >'B2X ','B2Y ','B2Z ', 48 : >'C2X ','C2Y ','C2Z ', 49 : >'D2X ','D2Y ','D2Z ', 50 : >'E2X ','E2Y ','E2Z '/ 51 : INTEGER TENR(NFAC),TNI(NFAC),TF2O(NFAC) 52 : c ddlfac correspond a quel enrichissement? 53 : DATA TENR/1,1,1,2,2,2,2,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3,3,3/ 54 : c ddlfac correspond a quel fonction d'enrichissement? 55 : DATA TNI/1,1,1,1,1,1,2,2,2,3,3,3,4,4,4,1,1,1,2,2,2,3,3,3,4,4,4/ 56 : c ddlfac correspond a quel ddlobl? 57 : DATA TF2O/1,2,3,1,2,3,1,2,3,1,2,3,1,2,3,1,2,3,1,2,3,1,2,3,1,2,3/ 58 : c tables pour mise en concordance des ddl 59 : INTEGER TOBL(NOBL),TFAC(NOBL,4,NBRMAX) 60 : INTEGER TIFAC(NOBL,4,NBRMAX) 61 : c fonctions de forme 62 : REAL*8 SHPWRK(4),SHPWR2(4) 63 : 64 : C_____________________________________________________________ 65 : C 66 : NBERR1=0 67 : 68 : C LECTURE 68 : CALL LIROBJ('MAILLAGE',IPT2in,1,IRETOU) 68 : CALL LIROBJ('CHPOINT',MCHPOI,1,IRETOU) 68 : CALL LIROBJ('MMODEL',MMODEL,1,IRETOU) 69 : IF(IERR.NE.0) RETURN 70 : 71 : SEGACT,MCOORD 72 : 73 : C_____________________________________________________________ 74 : C TRANSFO EN MCHELM DU CHPOINT avec ddl xfem 75 : * SEGACT,MCHPOI 76 : * NSOUP0 = IPCHP(/1) 76 : call CHAME1(0,MMODEL,MCHPOI,' ',IPCHEL,1) 77 : 78 : C_____________________________________________________________ 79 : C CHANGEMENT DE SUPPORT POUR LE CHPOINT DE SORTIE 80 : ipt2 = IPT2in 80 : call change(ipt2,1) 81 : c write(6,*) 'xpost2: retour de change (ipt2 doit etre actif)',ipt2 82 : 83 : C_____________________________________________________________ 84 : C ACTIVATION DU MODELE 85 : SEGACT,MMODEL 86 : NSOUS = KMODEL(/1) 87 : 88 : C_____________________________________________________________ 89 : c write(6,*) 'C INITIALISATION du mchpo2' 90 : C 91 : segact,MCHPOI 92 : NAT=JATTRI(/1) 93 : NSOUPO=1 94 : SEGINI,MCHPO2=MCHPOI 95 : segadj,MCHPO2 96 : SEGINI,MCHPO3=MCHPO2 97 : c on sait que l'on a seulement besoin de UX et UY en 2D 98 : c (+léger que dans xpost1.eso) 99 : NC=IDIM 100 : segini,MSOUP2,MSOUP3 101 : MCHPO2.IPCHP(1) = MSOUP2 102 : MCHPO3.IPCHP(1) = MSOUP3 103 : do ic2=1,IDIM 104 : MSOUP2.NOCOMP(ic2) = DDLOBL(ic2) 105 : MSOUP3.NOCOMP(ic2) = DDLOBL(ic2) 106 : MSOUP2.NOHARM(ic2) = IFOUR 107 : MSOUP3.NOHARM(ic2) = IFOUR 108 : enddo 109 : MSOUP2.IGEOC = IPT2 110 : MSOUP3.IGEOC = IPT2 111 : * segact,IPT2 112 : NBNN2 = IPT2.NUM(/1) 113 : NBELE2 = IPT2.NUM(/2) 114 : NBPT2 = NBELE2 115 : N = NBPT2 116 : C segini,MPOVA2,MPOVA3 117 : C bp: on ajoute IDEJVU pour ne caluler le champ de depl qu'une seule fois si le point arive sur le bord d un element 118 : segini,MPOVA2,MPOVA3,IDEJVU 119 : MSOUP2.IPOVAL = MPOVA2 120 : MSOUP3.IPOVAL = MPOVA3 121 : segdes,MSOUP2,MCHPO2,MCHPOI 122 : segdes,MSOUP3,MCHPO3 123 : 124 : C CHAMP DE DEPLACEMENT X-FEM 125 : MCHELM = IPCHEL 126 : SEGACT,MCHELM 127 : 128 : 129 : C_____________________________________________________________ 130 : C>>>>>>>>>>>>>>>>>>>>>>> BOUCLE SUR LES ZONES DU Modele 131 : DO 1000 ISOUS=1,NSOUS 132 : c write(6,*) '____________________________________' 133 : c write(6,*) 'xpost2: ZONE',ISOUS,' / ',NSOUS 134 : 135 : IMODEL = KMODEL(ISOUS) 136 : SEGACT,IMODEL 137 : 138 : C_____________________________________ 139 : C RECHERCHE DU MCHAM1 issu du MCHEX1 D ENRICHISSEMENT 140 : MCHAM1=0 141 : NBENR2=0 142 : NOBMOD=IVAMOD(/1) 143 : IF (NOBMOD.NE.0) THEN 144 : c boucle sur les objets contenus dans ivamod 145 : DO 1002 iobmo1=1,NOBMOD 146 : c write(6,*) 'objet=',iobmo1,'/',NOBMOD 147 : if((TYMODE(iobmo1)).ne.'MCHAML') goto 1002 148 : MCHEX1 = IVAMOD(iobmo1) 149 : segact,MCHEX1 150 : if((MCHEX1.TITCHE).ne.'ENRICHIS') goto 1003 151 : c REM. IMPORTANTE: a ce jour 1 seule zone XFEM, donc on prend ICHAML(isouX=1) 152 : c mais on ecrit ci dessous un prototype de la recherche de la zone du MCHEX1 153 : c associée a la zone du sous modele en cours via le pointeur du maillage 154 : nsouX = MCHEX1.ICHAML(/1) 155 : do isouX=1,nsouX 156 : if(IMAMOD.eq.IMACHE(isouX)) goto 1004 157 : enddo 158 : goto 1003 159 : 1004 continue 160 : MCHAM1 = MCHEX1.ICHAML(isouX) 161 : segact,MCHAM1 162 : NBENR2 = MCHAM1.IELVAL(/1) 163 : do IENR2=1,NBENR2 164 : MELVA1=MCHAM1.IELVAL(IENR2) 165 : segact,MELVA1 166 : enddo 167 : 1003 continue 168 : segdes,MCHEX1 169 : 1002 CONTINUE 170 : ENDIF 171 : c write(6,*) 'IMAMOD,IMACHE=',IMAMOD,IMACHE(isouX) 172 : c write(6,*) 'isouX/n,MCHAM1,NBENR2=',isouX,nsouX,MCHAM1,NBENR2 173 : 174 : C_____________________________________ 175 : C ACTIVATION DU SOUS CHPOINT devenu MCHELM ET DES COMPOSANTES MELVAL 176 : MCHAML = ICHAML(ISOUS) 177 : SEGACT,MCHAML 178 : N2= IELVAL(/1) 179 : DO I2=1,N2 180 : MELVAL = IELVAL(I2) 181 : SEGACT,MELVAL 182 : ENDDO 183 : 184 : c nbre de composante, de points 185 : NCOMP = N2 186 : 187 : 188 : C_____________________________________ 189 : c quelles sont les composantes obligatoires (=physiques) et ou sont elles? 190 : c on en deduit NC 191 : NC=0 192 : DO 1011 IOBL=1,NOBL 193 : MODDL2 = DDLOBL(IOBL) 194 : DO ICOMP=1,NCOMP 195 : MODDL = NOMCHE(ICOMP) 196 : * on a trouvé cette comp obl dans le chpoint d entree 197 : IF (MODDL.eq.MODDL2) THEN 198 : NC=NC+1 199 : TOBL(NC) = ICOMP 200 : GOTO 1011 201 : ENDIF 202 : ENDDO 203 : 1011 CONTINUE 204 : IF (NC.lt.NOBL) THEN 205 : DO IOBL=(NC+1),NOBL 206 : TOBL(IOBL) = 0 207 : ENDDO 208 : ENDIF 209 : 210 : c ...et les facultatives(enrichissement)? 211 : do i1=1,NOBL 212 : do i2=1,4 213 : do i3=1,NBRMAX 214 : TFAC(i1,i2,i3) = 0 215 : enddo 216 : enddo 217 : enddo 218 : NF=0 219 : IFAC = 0 220 : DO 1012 IFAC=1,NFAC 221 : MODDL2 = DDLFAC(IFAC) 222 : DO ICOMP=1,NCOMP 223 : MODDL = NOMCHE(ICOMP) 224 : * on a trouvé une comp fac qui va etre ajouté a la comp obl dans le chpoint de sortie 225 : IF (MODDL.eq.MODDL2) THEN 226 : NF=NF+1 227 : IOBL=TF2O(IFAC) 228 : INI =TNI (IFAC) 229 : IENR=TENR(IFAC) 230 : TFAC(IOBL,INI,IENR) = ICOMP 231 : TIFAC(IOBL,INI,IENR) = IFAC 232 : GOTO 1012 233 : ENDIF 234 : ENDDO 235 : * on n a pas trouvé la composante facultative 236 : IOBL=TF2O(IFAC) 237 : INI =TNI (IFAC) 238 : IENR=TENR(IFAC) 239 : TIFAC(IOBL,INI,IENR) = IFAC 240 : 1012 CONTINUE 241 : c write(*,*) 'TOBL et TFAC remplis' 242 : 243 : c recup du maillage 244 : MELEME = IMAMOD 245 : SEGACT,MELEME 246 : NBNN0=NUM(/1) 247 : NBELE0=NUM(/2) 248 : MELGEO = ITYPEL 249 : 250 : c initialisation du segment de travail 251 : NBNN = NBNN0 252 : segini,wrk4 253 : 254 : 255 : C_____________________________________ 256 : C>>>>>>>>>>>>> BOUCLE SUR LES ELEMENTS DU MODELE 257 : DO 2000 J0=1,NBELE0 258 : 259 : c XEL = coor des noeuds de l'ef J0 de MELEME 259 : call DOXE(XCOOR,idim,NBNN0,NUM,J0,XEL) 260 : 261 : 262 : C_____________________________________ 263 : C>>>>>>>>>>>>> BOUCLE SUR LES NOEUD DU CHP2 264 : DO 2002 J2=1,NBPT2 265 : 266 : * point J2 IDEJaVU dans un autre element => on saute 267 : if(IDEJVU(J2).eq.1) goto 2002 268 : 269 : C inoeu2 = IPT2.NUM(1,J2) 270 : c XPT2 = coor du pt J2 de IPT2 270 : call DOXE(XCOOR,idim,1,IPT2.NUM,J2,XPT2) 271 : 272 : c calcul iteratif des coord QSI du pt XPT2 dans le repere local de l'ef defini par XEL 272 : call qsijs(XEL,MELGEO,NBNN,idim,XPT2,SHPP,QSI,iret) 273 : if(iret.ne.0) goto 2002 274 : do i0=1,NBNN0 275 : * if (SHPP(1,i0).lt.0.) goto 2002 276 : * a revoir car pb lorsqu'on tombe "pile" sur un bord 277 : if(SHPP(1,i0).lt.(-1.E-7)) goto 2002 278 : enddo 279 : 280 : * On a trouvé que le point J2 appartenait a l'EF J0 281 : IDEJVU(J2) = 1 282 : c write(6,*) '_______________________________' 283 : c write(*,*) 'ELEM ',J0,'/',NBELE0,' point ',J2,'/',NBPT2 284 : c write(6,*) 'QSI= ',(QSI(iou),iou=1,3) 285 : c write(6,*) 'SHPP=',(SHPP(1,iou),iou=1,4) 286 : 287 : C______________________ 288 : C On commence par interpoler la partie relative aux ddls (UX,UY) 289 : DO IC2=1,NC 290 : MPOVA2.VPOCHA(J2,IC2) = 0.d0 291 : ICOMP = TOBL(IC2) 292 : MELVAL = IELVAL(ICOMP) 293 : DO I0=1,NBNN0 294 : MPOVA2.VPOCHA(J2,IC2) = MPOVA2.VPOCHA(J2,IC2) 295 : $ + ( SHPP(1,I0) * VELCHE(I0,J0) ) 296 : ENDDO 297 : MPOVA3.VPOCHA(J2,IC2) = MPOVA2.VPOCHA(J2,IC2) 298 : ENDDO 299 : 300 : C______________________ 301 : C le noeud inoeu (I0^ieme noeud de l ef J0) est il IENR2-enrichi? 302 : if(NBENR2.eq.0) goto 2002 303 : DO 3000 I0=1,NBNN0 304 : inoeu = NUM(I0,J0) 305 : 306 : DO 3001 IENR2=1,NBENR2 307 : 308 : MELVA1=MCHAM1.IELVAL(IENR2) 309 : MLREEL=MELVA1.IELCHE(I0,J0) 310 : 311 : C si ce noeud n est pas IENR2-enrichi on ne fait rien 312 : IF(MLREEL.eq.0) GOTO 3001 313 : 314 : C si ce noeud est enrichi,... 315 : c write(6,*) 'I0,inoeu,IENR2,MLREEL=', 316 : c $ I0,inoeu,IENR2,MLREEL 317 : 318 : c ...on calcule les fonctions d enrichissement : 319 : 320 : c------------pour IENR=1, fonction H, ddl AX et AY 321 : IF (IENR2.eq.1) THEN 322 : NBNI = 1 323 : SHPWRK(1) = 1.d0 324 : SHPWR2(1) = -1.d0 325 : ENDIF 326 : c------------fin du cas IENR=1, fonction H 327 : 328 : c------------pour IENR>1, 4 fonctions F 329 : IF (IENR2.ge.2) THEN 330 : NBNI = 4 331 : PHIX = 0.d0 332 : PSIX = 0.d0 333 : SEGACT,MLREEL 334 : c WRITE(6,*) 'XPSI1=',(PROG(IOU),IOU=1,4) 335 : do iii0 = 1,NBNN0 336 : XSHPP = SHPP(1,iii0) 337 : XPSI1 = PROG(iii0) 338 : PSIX = PSIX + (XSHPP * XPSI1) 339 : enddo 340 : SEGDES,MLREEL 341 : * write(*,*) 'PHIX,PSIX=',PHIX,PSIX 342 : RX = ( (PSIX**2) + (PHIX**2) )**0.5 343 : * IF(RX.LT.XTOL) THEN 344 : IF (PSIX.ge.0.) THEN 345 : THETAX = 0. 346 : THETA2 = 0. 347 : ELSE 348 : * THETAX = HX * ((XPI/2.) - (ATAN(PSIX/ABS(PHIX)))) 349 : THETAX = XPI 350 : THETA2 = -1.*XPI 351 : ENDIF 352 : RX05= (RX**0.5) 353 : c write(6,*) 'PSIX,PHIX,RX05=',PSIX,PHIX,RX05 354 : SIN1T = SIN(THETAX) 355 : C COS1T = COS(THETAX) 356 : SIN05T = SIN(THETAX/2.) 357 : COS05T = COS(THETAX/2.) 358 : c if (IENR2.lt.2) then 359 : c SHPWRK(1) = 0.d0 360 : c else 361 : SHPWRK(1) = (RX05) * SIN05T 362 : c endif 363 : SHPWRK(2) = (RX05) * COS05T 364 : SHPWRK(3) = (RX05) * SIN05T * SIN1T 365 : SHPWRK(4) = (RX05) * COS05T * SIN1T 366 : SIN12 = SIN(THETA2) 367 : C COS12 = COS(THETA2) 368 : SIN052 = SIN(THETA2/2.) 369 : COS052 = COS(THETA2/2.) 370 : c if (IENR2.lt.2) then 371 : c SHPWR2(1) = 0.d0 372 : c else 373 : SHPWR2(1) = (RX05) * SIN052 374 : c endif 375 : SHPWR2(2) = (RX05) * COS052 376 : SHPWR2(3) = (RX05) * SIN052 * SIN12 377 : SHPWR2(4) = (RX05) * COS052 * SIN12 378 : ENDIF 379 : c------------fin du cas IENR>1, fonction F 380 : 381 : c on ajoute les fonctions d enrichissement 382 : DO 3900 IC2=1,NC 383 : DO 3900 INI=1,NBNI 384 : ICOMP = TFAC(IC2,INI,IENR2) 385 : c cas ou on a pas trouvé cette composante dans cette zone du 386 : c chpoint solution => AVERTISSEMENT puis on saute 387 : if (ICOMP.eq.0) then 388 : NBERR1=NBERR1+1 389 : if (IIMPI.ge.1) then 390 : write(IOIMP,991) DDLFAC(icomp),inoeu 391 : 991 format(2X,'ABSENCE DE LA COMPOSANTE ',A4,' AU NOEUD ', 392 : $ I6,' DANS LE CHPOINT FOURNI a XFEM FISS') 393 : endif 394 : goto 3900 395 : endif 396 : MELVAL = IELVAL(ICOMP) 397 : MPOVA2.VPOCHA(J2,IC2) = MPOVA2.VPOCHA(J2,IC2) 398 : $ + ( (SHPP(1,I0) * SHPWRK(INI)) * VELCHE(I0,J0) ) 399 : MPOVA3.VPOCHA(J2,IC2) = MPOVA3.VPOCHA(J2,IC2) 400 : $ + ( (SHPP(1,I0) * SHPWR2(INI)) * VELCHE(I0,J0) ) 401 : 3900 CONTINUE 402 : 403 : 404 : 3001 CONTINUE 405 : C<<<<<<<<< FIN DE BOUCLE SUR LES enrichissements 406 : 3000 CONTINUE 407 : C<<<<<<<<< FIN DE BOUCLE SUR LES noeuds 408 : 409 : 410 : 411 : 2002 CONTINUE 412 : C<<<<<<<<<<<<<<< FIN DE BOUCLE SUR LES noeud du chp2 413 : 2000 CONTINUE 414 : C<<<<<<<<<<<<<<< FIN DE BOUCLE SUR LES elements 415 : 416 : IF (MCHAM1.ne.0) THEN 417 : do IENR2=1,NBENR2 418 : MELVA1=MCHAM1.IELVAL(IENR2) 419 : segdes,MELVA1 420 : enddo 421 : segdes,MCHAM1 422 : ENDIF 423 : DO I2=1,N2 424 : MELVAL = IELVAL(I2) 425 : SEGDES,MELVAL 426 : ENDDO 427 : SEGDES,MCHAML 428 : SEGDES,MELEME,IMODEL 429 : segsup,wrk4 430 : 431 : 1000 CONTINUE 432 : C<<<<<<<<<<<<<<<<<<<<<<<<< FIN DE BOUCLE SUR LES ZONES 433 : 434 : SEGDES,MMODEL,MCHELM 435 : SEGDES,MPOVA2,MPOVA3 436 : SEGDES,IPT2 437 : SEGSUP,IDEJVU 438 : 439 : c --cas ou on a une ou des erreurs-- 440 : IF (NBERR1.gt.0) THEN 441 : write(IOIMP,*) 'OPERATEUR XFEM FISS : ABSENCE DANS LE CHPOINT ', 442 : & 'DEPLACEMENT DE CERTAINES INCONNUES ATTENDUES PAR LE MODELE' 443 : ENDIF 444 : 445 : C_____________________________________________________________ 446 : C 447 : C ECRITURE 447 : CALL ECROBJ('CHPOINT',MCHPO2) 447 : CALL ECROBJ('CHPOINT',MCHPO3) 448 : 449 : 450 : RETURN 451 : END 452 : 453 : 454 : 455 : 456 : 457 : 458 : 459 : 460 :
© Cast3M 2003 - All rights reserved.
Disclaimer