1 : C XPOST1 SOURCE BP208322 13/11/13 21:15:25 7865 2 : c 3 : SUBROUTINE XPOST1(IPCHP1,IPMOD1,IPCHP2) 4 : c 5 : c Construit un Chpoint avec des ddl "physiques" (UX UY) en 6 : c RECOmbinant les ddl Xfem (UX AX B1X C1X ...) 7 : c -> Utile pour le tracé de déformées 8 : c 9 : C SYNTAXE : 10 : c chpo2 = XFEM 'RECO' chpo1 MODEL_XFEM 11 : C 12 : c 13 : IMPLICIT INTEGER(I-N) 14 : IMPLICIT REAL*8 (A-H,O-Z) 15 : C 16 : PARAMETER (XTOL1=1.d-7) 17 : PARAMETER (NBRMAX=5) 18 : 19 : C 20 : C SEGMENTS INCLUDE 21 : -INC CCREEL 22 : -INC CCOPTIO 23 : -INC SMCOORD 24 : -INC SMELEME 25 : -INC SMCHPOI 26 : -INC SMCHAML 27 : -INC SMMODEL 28 : -INC SMLREEL 29 : c 30 : POINTEUR MCHEX1.MCHELM,MPHI1.MPOVAL,MPSI1.MPOVAL 31 : C 32 : 33 : c SEGMENT MRACC 34 : c INTEGER TLENR(NBENR1,NBPTS) 35 : c ENDSEGMENT 36 : SEGMENT MRACC 37 : INTEGER TLENR(NBENR1,NBPTS) 38 : INTEGER IPOCHA(NBPTS,NCOMP1) 39 : ENDSEGMENT 40 : c SEGMENT MRACC 41 : c INTEGER NOD2PT(NBPTS) 42 : c INTEGER TLENR(NBENR1,NBPT1) 43 : c INTEGER IPOCHA(NCOMP1,NBPT1) 44 : c ENDSEGMENT 45 : C 46 : C 47 : C_____________________________________________________________ 48 : C INITIALISATION DES INCONNUES obligatoires et facultatives 49 : PARAMETER (NOBL=4,NFAC=27) 50 : CHARACTER*4 DDLOBL(NOBL),DDLFAC(NFAC),MODDL,MODDL2 51 : DATA DDLOBL/'UX ','UY ','UZ ','LX '/ 52 : DATA DDLFAC/'AX ','AY ','AZ ', 53 : >'B1X ','B1Y ','B1Z ', 54 : >'C1X ','C1Y ','C1Z ', 55 : >'D1X ','D1Y ','D1Z ', 56 : >'E1X ','E1Y ','E1Z ', 57 : >'B2X ','B2Y ','B2Z ', 58 : >'C2X ','C2Y ','C2Z ', 59 : >'D2X ','D2Y ','D2Z ', 60 : >'E2X ','E2Y ','E2Z '/ 61 : INTEGER TENR(NFAC),TNI(NFAC),TF2O(NFAC) 62 : c ddlfac correspond a quel enrichissement? 63 : 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/ 64 : c ddlfac correspond a quel fonction d'enrichissement? 65 : 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/ 66 : c ddlfac correspond a quel ddlobl? 67 : 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/ 68 : c tables pour mise en concordance des ddl 69 : INTEGER TOBL(NOBL),TFAC(NOBL,4,NBRMAX) 70 : INTEGER TIFAC(NOBL,4,NBRMAX) 71 : c fonctions de forme 72 : REAL*8 SHPWRK(4) 73 : 74 : C_____________________________________________________________ 75 : C ACTIVATION ... 76 : 77 : C ...DU MODELE 78 : MMODEL=IPMOD1 79 : SEGACT,MMODEL 80 : NSOUS = KMODEL(/1) 81 : 82 : C ...DU CHPOINT avec ddl xfem 83 : MCHPOI=IPCHP1 84 : SEGACT,MCHPOI 85 : NSOUP0 = IPCHP(/1) 86 : 87 : C_____________________________________________________________ 88 : C TRAVAIL PRELIMINAIRE SUR L ENRICHISSEMENT' 89 : C on les transforme en 2 chpoints psi et phi 90 : c avec le type d enrichissement pour composante 91 : 92 : SEGACT,MCOORD 93 : NBPTS = (XCOOR(/1)) / (idim+1) 94 : NBENR1= NBRMAX 95 : NCOMP1=NFAC 96 : SEGINI,MRACC 97 : N=NBPTS 98 : NC=NBRMAX 99 : segini,MPHI1,MPSI1 100 : NBPTS = 0 101 : NBENR2= 0 102 : 103 : c ----boucle sur les zones du modele 104 : DO 0001 ISOUS=1,NSOUS 105 : c write(6,*) 'zone du modele sISOUS=',ISOUS,'/',NSOUS 106 : IMODEL = KMODEL(ISOUS) 107 : SEGACT,IMODEL 108 : NOBMOD=IVAMOD(/1) 109 : IF(NOBMOD.NE.0) THEN 110 : 111 : c -------boucle sur les objet ivamod jusqu'a trouver un mchelm d'enrichissement 112 : DO 0002 iobmo1=1,NOBMOD 113 : if((TYMODE(iobmo1)).ne.'MCHAML') goto 0002 114 : MCHEX1 = IVAMOD(iobmo1) 115 : segact,MCHEX1 116 : if((MCHEX1.TITCHE).ne.'ENRICHIS') goto 0003 117 : 118 : * on a trouvé une zone enrichie => on crée une table d'enrichissment nodal 119 : N1 = MCHEX1.ICHAML(/1) 120 : c ---------boucle sur les maillages elementaire du chaml d'enrichissement 121 : do 0010 i1=1,N1 122 : if((MCHEX1.INFCHE(i1,4)).ne.0) goto 0010 123 : MELEME = MCHEX1.IMACHE(i1) 124 : MCHAM1 = MCHEX1.ICHAML(i1) 125 : segact,MCHAM1,MELEME 126 : NBNN = NUM(/1) 127 : NBELEM = NUM(/2) 128 : N2=MCHAM1.IELVAL(/1) 129 : c ------------boucle sur les enrichissements possibles 130 : do i2=1,N2 131 : MELVAL=MCHAM1.IELVAL(i2) 132 : segact,MELVAL 133 : do j=1,NBELEM 134 : do i=1,NBNN 135 : inoeu = NUM(i,j) 136 : NBPTS = max(NBPTS,inoeu) 137 : MLREEL = IELCHE(i,j) 138 : c write(6,*) 'i,j,inoeu,mlreel',i,j,inoeu,MLREEL,NBPTS 139 : if (MLREEL.ne.0) then 140 : TLENR(i2,inoeu) = MLREEL 141 : NBENR2 = max(NBENR2,i2) 142 : segact,MLREEL 143 : if (i2.eq.1) then 144 : MPHI1.VPOCHA(inoeu,i2) = PROG(i) 145 : else 146 : MPSI1.VPOCHA(inoeu,i2) = PROG(i) 147 : MPHI1.VPOCHA(inoeu,i2) = PROG(NBNN+i) 148 : endif 149 : c avec MPSI1 et MPHI1, on peut fermer MLREEL 150 : segdes,MLREEL 151 : endif 152 : enddo 153 : enddo 154 : segdes,MELVAL 155 : enddo 156 : c ------------fin de boucle sur les enrichissements possibles 157 : segdes,MCHAM1,MELEME 158 : 0010 continue 159 : c ---------fin de boucle sur les maillages elementaire du chaml d'enrichissement 160 : 161 : 0003 CONTINUE 162 : segdes,MCHEX1 163 : 0002 CONTINUE 164 : c -------fin de boucle sur les objet ivamod 165 : 166 : ENDIF 167 : SEGDES,IMODEL 168 : 0001 CONTINUE 169 : c ----fin de boucle sur les zones du modele 170 : 171 : N=NBPTS 172 : * NC=NBENR1 173 : NC=NBENR2 174 : segadj,MPHI1,MPSI1 175 : 176 : SEGDES,MMODEL 177 : 178 : C_____________________________________________________________ 179 : C INITIALISATION 180 : c 181 : c ...du mchpo2 182 : NAT=JATTRI(/1) 183 : NSOUPO=NSOUP0 184 : SEGINI,MCHPO2=MCHPOI 185 : C NPOCHA=0 186 : 187 : 188 : C_____________________________________________________________ 189 : C>>>>>>>>>>>>>>>>>>>>>>> BOUCLE SUR LES ZONES DU CHPOINT 190 : DO 1000 ISOUP=1,NSOUP0 191 : c write(*,*) 'xpost1: ZONE ',ISOUP,' / ',NSOUP0 192 : 193 : C_____________________________________ 194 : C ACTIVATION DU SOUS CHPOINT 195 : MSOUPO = IPCHP(ISOUP) 196 : SEGACT,MSOUPO 197 : MELEME = IGEOC 198 : MPOVAL = IPOVAL 199 : SEGACT,MELEME,MPOVAL 200 : 201 : c nbre de composante, de points 202 : NCOMP = NOCOMP(/2) 203 : NBPT0 = NUM(/2) 204 : NC = NCOMP 205 : 206 : C INITIALISATION du MSOUP2 (avec NC =celui de départ) 207 : SEGINI,MSOUP2=MSOUPO 208 : MCHPO2.IPCHP(ISOUP)=MSOUP2 209 : 210 : c quelles sont les composantes obligatoires (=physiques) et ou sont elles? 211 : c on en deduit NC 212 : NC=0 213 : DO 1001 IOBL=1,NOBL 214 : MODDL2 = DDLOBL(IOBL) 215 : DO ICOMP=1,NCOMP 216 : MODDL = NOCOMP(ICOMP) 217 : * on a trouvé cette comp obl dans le chpoint d entree 218 : IF(MODDL.eq.MODDL2) THEN 219 : NC=NC+1 220 : MSOUP2.NOCOMP(NC) = MODDL2 221 : MSOUP2.NOCONS(NC) = NOCONS(ICOMP) 222 : TOBL(NC) = ICOMP 223 : GOTO 1001 224 : ENDIF 225 : ENDDO 226 : 1001 CONTINUE 227 : IF(NC.lt.NOBL) THEN 228 : DO IOBL=(NC+1),NOBL 229 : TOBL(IOBL) = 0 230 : ENDDO 231 : ENDIF 232 : 233 : c ...et les facultatives(enrichissement)? 234 : do i1=1,NOBL 235 : do i2=1,4 236 : do i3=1,NBRMAX 237 : TFAC(i1,i2,i3) = 0 238 : enddo 239 : enddo 240 : enddo 241 : NF=0 242 : IFAC = 0 243 : DO 1002 IFAC=1,NFAC 244 : MODDL2 = DDLFAC(IFAC) 245 : DO ICOMP=1,NCOMP 246 : MODDL = NOCOMP(ICOMP) 247 : * on a trouvé une comp fac qui va etre ajouté a la comp obl dans le chpoint de sortie 248 : IF(MODDL.eq.MODDL2) THEN 249 : NF=NF+1 250 : IOBL=TF2O(IFAC) 251 : INI =TNI (IFAC) 252 : IENR=TENR(IFAC) 253 : TFAC(IOBL,INI,IENR) = ICOMP 254 : TIFAC(IOBL,INI,IENR) = IFAC 255 : GOTO 1002 256 : ENDIF 257 : ENDDO 258 : * on n a pas trouvé la composante facultative 259 : IOBL=TF2O(IFAC) 260 : INI =TNI (IFAC) 261 : IENR=TENR(IFAC) 262 : TIFAC(IOBL,INI,IENR) = IFAC 263 : 1002 CONTINUE 264 : 265 : 266 : C INITIALISATION du MSOUP2 (avec NC ajusté) et du MPOVA2 267 : segadj,MSOUP2 268 : N = VPOCHA(/1) 269 : SEGINI,MPOVA2 270 : MSOUP2.IPOVAL = MPOVA2 271 : 272 : SEGDES,MSOUPO,MSOUP2 273 : c actifs a ce stade : MCHPO2,MCHPOI,MELEME,MPOVAL,MPOVA2 274 : 275 : C_____________________________________ 276 : C>>>>>>>>>>>>> BOUCLE SUR LES POINTS 277 : DO 2000 J=1,NBPT0 278 : c write(*,*) 'point ',J,' / ',NBPT0,NC 279 : 280 : C______________________ 281 : C On commence par recopier les valeurs obligatoires (UX,UY) 282 : DO IC2=1,NC 283 : ICOMP = TOBL(IC2) 284 : MPOVA2.VPOCHA(J,IC2) = VPOCHA(J,ICOMP) 285 : ENDDO 286 : 287 : C______________________ 288 : C ce noeud est il IENR2-enrichi? 289 : inoeu = NUM(1,J) 290 : 291 : DO 3000 IENR2=1,NBENR2 292 : 293 : MLREEL = TLENR(IENR2,inoeu) 294 : 295 : C si ce noeud n est pas enrichi on ne fait rien 296 : IF(MLREEL.eq.0) GOTO 3000 297 : 298 : C si ce noeud est enrichi, 299 : c on recupere les infos relatives a l enrichissement 300 : 301 : c on calcule les fonctions d enrichissement 302 : 303 : c------------pour IENR=1, fonction H, ddl AX et AY 304 : IF (IENR2.eq.1) THEN 305 : PHIX = MPHI1.VPOCHA(inoeu,IENR2) 306 : NBNI = 1 307 : SHPWRK(1) = SIGN(1.D0,PHIX) 308 : if (abs(phix).lt.xtol1) then 309 : SHPWRK(1) = 0.d0 310 : else 311 : SHPWRK(1) = SIGN(1.D0,PHIX) 312 : endif 313 : ENDIF 314 : c------------fin du cas IENR=1, fonction H 315 : 316 : c------------pour IENR>1, 4 fonctions F 317 : IF (IENR2.ge.2) THEN 318 : PHIX = MPHI1.VPOCHA(inoeu,IENR2) 319 : PSIX = MPSI1.VPOCHA(inoeu,IENR2) 320 : if (ABS(phix).lt.xtol1) then 321 : HX = 0.d0 322 : else 323 : HX = SIGN(1.D0,PHIX) 324 : endif 325 : RX = ( (PSIX**2.) + (PHIX**2.) )**0.5 326 : IF (RX.LT.XTOL1) THEN 327 : THETAX = 0.d0 328 : ELSE 329 : THETAX = HX * ((XPI/2.) - (ATAN(PSIX/ABS(PHIX)))) 330 : ENDIF 331 : SIN1T = SIN(THETAX) 332 : C COS1T = COS(THETAX) 333 : SIN05T = SIN(THETAX/2.) 334 : COS05T = COS(THETAX/2.) 335 : NBNI = 4 336 : SHPWRK(1) = (RX**0.5) * SIN05T 337 : SHPWRK(2) = (RX**0.5) * COS05T 338 : SHPWRK(3) = (RX**0.5) * SIN05T * SIN1T 339 : SHPWRK(4) = (RX**0.5) * COS05T * SIN1T 340 : ENDIF 341 : c------------fin du cas IENR>1, fonction F 342 : 343 : c on ajoute les fonctions d enrichissement 344 : DO 3900 IC2=1,NC 345 : DO 3900 INI=1,NBNI 346 : ICOMP = TFAC(IC2,INI,IENR2) 347 : IFAC = TIFAC(IC2,INI,IENR2) 348 : c 1ere fois qu on voit le besoin de (composante + noeud) 349 : if(IPOCHA(inoeu,IFAC).eq.0) IPOCHA(inoeu,IFAC)=1 350 : c cas ou on a pas trouvé cette composante dans cette zone du 351 : c chpoint solution => on saute simplement 352 : if(ICOMP.eq.0) goto 3900 353 : MPOVA2.VPOCHA(J,IC2) = MPOVA2.VPOCHA(J,IC2) 354 : $ + ( SHPWRK(INI) * VPOCHA(J,ICOMP) ) 355 : IPOCHA(inoeu,IFAC) = IPOCHA(inoeu,IFAC)+1 356 : 3900 CONTINUE 357 : 358 : 3000 CONTINUE 359 : C<<<<<<<<< FIN DE BOUCLE SUR LES enrichissements 360 : 361 : 362 : 363 : 2000 CONTINUE 364 : C<<<<<<<<<<<<<<< FIN DE BOUCLE SUR LES POINTS 365 : 366 : SEGDES,MPOVA2,MPOVAL,MELEME 367 : 368 : 1000 CONTINUE 369 : C<<<<<<<<<<<<<<<<<<<<<<<<< FIN DE BOUCLE SUR LES ZONES 370 : 371 : SEGDES,MCHPOI,MCHPO2 372 : segsup,MPHI1,MPSI1 373 : 374 : 375 : C_____________________________________________________________ 376 : C PETIT AVERTISSEMENT S'IL MANQUE DES ENRICHISSEMENTS 377 : C 378 : C --initialisation du nombre d erreur sur le noms de composantes 379 : NBERR1=0 380 : 381 : c --recherche des erreurs 382 : DO inoeu=1,NBPTS 383 : DO icomp1=1,NCOMP1 384 : c -cas ou on a pas besoin de cette composante en ce noeud : 385 : c (IPOCHA(inoeu,icomp1).EQ.0) 386 : c -Cas ou on a trouvé un noeud enrichi sans la composante associée 387 : c dans le chpoint => avertissement : 388 : IF (IPOCHA(inoeu,icomp1).EQ.1) THEN 389 : NBERR1=NBERR1+1 390 : if (IIMPI.ge.1) then 391 : write(IOIMP,991) DDLFAC(icomp1),inoeu 392 : 991 format(2X,'ABSENCE DE LA COMPOSANTE ',A4,' AU NOEUD ',I6, 393 : $ ' DANS LE CHPOINT FOURNI a XFEM RECO') 394 : endif 395 : ENDIF 396 : ENDDO 397 : ENDDO 398 : c --cas ou on a une ou des erreurs-- 399 : IF (NBERR1.gt.0) THEN 400 : write(IOIMP,*) 'OPERATEUR XFEM RECO : ABSENCE DANS LE CHPOINT ', 401 : & 'DEPLACEMENT DE CERTAINES INCONNUES ATTENDUES PAR LE MODELE' 402 : ENDIF 403 : 404 : segsup,MRACC 405 : 406 : 407 : C_____________________________________________________________ 408 : C ON RETOURNE LE CHPOINT 409 : 410 : IPCHP2=MCHPO2 411 : 412 : 413 : RETURN 414 : END 415 : 416 : 417 : 418 : 419 : 420 : 421 : 422 : 423 : 424 : 425 : 426 :
© Cast3M 2003 - All rights reserved.
Disclaimer