1 : C TRJMA2 SOURCE PV 08/09/11 21:17:12 6150 2 : SUBROUTINE TRJMA2(XARI,XDEP,UELEM,DTINT, 3 : * IDIM,ICONT,ITYEL,IART,INOEU) 4 : C 5 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 6 : C CALCUL DES TRAJECTOIRES ET 7 : C INTERSECTION TRAJECTOIRE-ELEMENT ( calcul analytique pour une 8 : C formulation EFMH) 9 : C 10 : C ENTREES 11 : C XDEP POSITION INITIALE 12 : C UELEM FLUX AUX FACES 13 : C IDIM DIMENSION DE L ESPACE 14 : C ITYEL TYPE DE L ELEMENT 15 : C 16 : C SORTIES 17 : C XARI POSITION A L INTERSECTION 18 : C DTINT TEMPS ECOULE JUSQU A LA SORTIE 19 : C ICONT NUMERO DE LA PORTE DE SORTIE 20 : C IART DIFFERENT DE 0 ON EST SORTI DE L ELEMENT 3D PAR UNE ARETE 21 : C INOEU DIFFERENT DE 0 ON EST SORTI DE L ELEMENT NOEUD 22 : C 23 : C TYPES D ELEMENTS CONSIDERES 24 : C TRI3 QUA4 CUB8 PRI6 TET 25 : C 4 8 14 16 23 26 : C 27 : C Auteur Patrick Meyniel 28 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 29 : C 30 : IMPLICIT INTEGER(I-N) 31 : IMPLICIT REAL*8 (A-H,O-Z) 32 : LOGICAL TEST 33 : C 34 : DIMENSION IJ(2),ZI(2),IFAT3(3),YI(3) 35 : DIMENSION IFAQ4(4),IJKL(3,3),IJT4(2,3),IFAT4(4),IFAC6(4),IFAC8(6) 36 : DIMENSION XARI(3),XDEP(3),UELEM(6),X(3),Y(3) 37 : C 38 : TEST(DTINT,DTI)=(DTINT*DTI*(DTINT-DTI).GE.0.D0) 39 : DATA IJKL/1,2,3,2,1,3,3,2,1/ 40 : DATA IJ /2,1/ 41 : DATA IJT4 /2,3, 1,3, 1,2/ 42 : DATA IFAT3 /3,1,2/ 43 : DATA IFAQ4 /1,3,4,2/ 44 : DATA IFAC6 /6,4,1,2/ 45 : DATA IFAC8 /6,4,3,5,1,2/ 46 : DATA IFAT4 /4,1,2,3/ 47 : DATA YI/5.D-1,5.D-1,-1.D0/ 48 : DATA ZI/-1.D0,1.D0/ 49 : 50 : C 51 : ICONT=0 52 : IART=0 53 : INOEU=0 54 : ITEST=0 55 : EPS=1.D-14 56 : DTINT=1.D30 57 : DO 50 I=1,IDIM 58 : XARI(I)=XDEP(I) 59 : 50 CONTINUE 60 : 61 : C 62 : C*** TRIANGLES TRI3 * 63 : C * * 64 : C FACES 3 * * 2 65 : C * * 66 : C ****** 67 : C 1 68 : IF(ITYEL.EQ.4)THEN 69 : IDIM=2 70 : C 71 : C A PARTIR DES FLUX ON RECHERCHE LE POINT DE SORTIE DE 72 : C LA TRAJECTOIRE DANS L'ELEMENT 73 : C 74 : C ON TESTE L'INTERSECTION AVEC LA FACE 2 ET TOUTES LES 75 : C FACES AYANT UN FLUX POSITIF 76 : C FACE 2 77 : C 78 : IF(UELEM(2).GT.0.D0)THEN 79 : DTINT=(XDEP(2)+XDEP(1)-1)/(UELEM(1)+UELEM(3)) 80 : XARI(1)=XDEP(1)-(UELEM(3)*DTINT) 81 : XARI(2)=XDEP(2)-(UELEM(1)*DTINT) 82 : ICONT=2 83 : ENDIF 84 : C 85 : C FACE 1 ET 3 86 : C 87 : DO 1 I=1,IDIM 88 : C 89 : IF(UELEM(IFAT3(I)).GT.0.D0)THEN 90 : DTI=XDEP(I)/UELEM(IFAT3(I)) 91 : C 92 : C ON NE GARDE QUE LE PLUS PETIT TEMPS POSITIF 93 : C 94 : IF(TEST(DTINT,DTI))THEN 95 : DTINT=DTI 96 : C ON GARDE LES COORD DU POINT ET LA FACE DE SORTIE 97 : C CORRESPONDANTS AU TEMPS LE PLUS PETIT 98 : XARI(IJ(I))=XDEP(IJ(I))-(UELEM(IFAT3(IJ(I))) 99 : * *DTINT) 100 : XARI(I)=0.D0 101 : ICONT=IFAT3(I) 102 : ENDIF 103 : ENDIF 104 : 1 CONTINUE 105 : C 106 : C ON CHERCHE SI ON A AFFAIRE A UN NOEUD DU MAILLAGE 107 : C 108 : 109 : DO 30 I=1,IDIM 110 : IF(ABS(XARI(I)).LT.EPS)ITEST=ITEST+1 111 : IF(ABS(1.D0-XARI(I)).LT.EPS)ITEST=ITEST+1 112 : 30 CONTINUE 113 : C 114 : C 115 : C 116 : C*** QUADRANGLE QUA4 QUA9 3 117 : C ***** 118 : C * * 119 : C FACES 4 * * 2 120 : C ***** 121 : C 1 122 : ELSEIF(ITYEL.EQ.8)THEN 123 : IDIM=2 124 : Q=MAX(ABS(UELEM(1)),ABS(UELEM(2)),ABS(UELEM(3)),ABS(UELEM(4))) 125 : C SI LES FLUX OPPOSES NE SE COMPENSENT PAS 126 : C 127 : X(1)=UELEM(3)+UELEM(1) 128 : X(2)=UELEM(2)+UELEM(4) 129 : Y(1)=UELEM(3)-UELEM(1) 130 : Y(2)=UELEM(2)-UELEM(4) 131 : C 132 : IF(ABS(X(1))/Q.GT.1D-10)THEN 133 : C ON TESTE TOUTES LES FACES DE SORTIE QUI ONT UN FLUX POSITIF 134 : 135 : DO 2 I=1,IDIM 136 : DO 3 I1=1,IDIM 137 : IF(UELEM(IFAQ4(2*I+I1-2)).GT.0.D0)THEN 138 : VAR1=ZI(I1)+(Y(I)/X(I)) 139 : VAR2=XDEP(IJ(I))+(Y(I)/X(I)) 140 : C ON ELIMINE LES CAS OÙ LE LOG EST NEG OU NON DEFINI 141 : IF((VAR1/VAR2).LE.0.D0)GOTO 7 142 : DTI=(4.D0/X(I))*LOG(VAR1/VAR2) 143 : IF(ABS(LOG(VAR1/VAR2)).LT.1.D-15)DTI=0.D0 144 : C ON GARDE LE PLUS PETIT PAS DE TEMPS 145 : IF(TEST(DTINT,DTI))THEN 146 : DTINT=DTI 147 : I2=I 148 : I3=I1 149 : ENDIF 150 : ENDIF 151 : 7 CONTINUE 152 : 3 CONTINUE 153 : 2 CONTINUE 154 : C 155 : C ON EXHIBE LES COORDONNEES CORRESPONDANT AU PLUS PETIT TEMPS 156 : XARI(I2)=(XDEP(I2)*EXP(25.D-2*X(IJ(I2))*DTINT))- 157 : * ((Y(IJ(I2))/X(IJ(I2)))*(1.-EXP(25.D-2*X(IJ(I2))*DTINT))) 158 : C 159 : C SINON (LES FLUX OPPOSES SE COMPENSENT) 160 : C 161 : ELSE 162 : C 163 : C 164 : C ON TESTE TOUTES LES FACES AYANT UN FLUX POSITIF 165 : C 166 : DO 4 I=1,IDIM 167 : DO 5 I1=1,IDIM 168 : IF(UELEM(IFAQ4(2*I+I1-2)).GT.0.D0)THEN 169 : C 170 : DTI=(ZI(I1)-XDEP(IJ(I)))/(25.D-2*Y(I)) 171 : C ON GARDE LE PLUS PETIT TEMPS POSITIF 172 : IF(TEST(DTINT,DTI))THEN 173 : DTINT=DTI 174 : I2=I 175 : I3=I1 176 : ENDIF 177 : ENDIF 178 : 5 CONTINUE 179 : 4 CONTINUE 180 : C 181 : C ON GARDE LES COORDONNEES CORRESPONDANT AU PLUS PETIT TEMPS 182 : XARI(I2)=XDEP(I2)+25.D-2*Y(IJ(I2))*DTINT 183 : ENDIF 184 : XARI(IJ(I2))=ZI(I3) 185 : C ON EXHIBE LE NUMERO DE LA FACE DE SORTIE 186 : ICONT=IFAQ4((2*I2)+I3-2) 187 : C 188 : C ON CHERCHE SI ON A AFFAIRE A UN NOEUD DU MAILLAGE 189 : 190 : DO 35 I=1,IDIM 191 : IF(ABS(1.D0-XARI(I)).LT.EPS)ITEST=ITEST+1 192 : IF(ABS(1.D0+XARI(I)).LT.EPS)ITEST=ITEST+1 193 : 35 CONTINUE 194 : C 195 : C 196 : C*** CUBE CUB8 197 : ELSEIF(ITYEL.EQ.14)THEN 198 : IDIM=3 199 : Q=MAX(ABS(UELEM(1)),ABS(UELEM(2)),ABS(UELEM(3)), 200 : * ABS(UELEM(4)),ABS(UELEM(5)),ABS(UELEM(6))) 201 : X(1)=125.D-3*(UELEM(4)+UELEM(6)) 202 : X(2)=125.D-3*(UELEM(3)+UELEM(5)) 203 : X(3)=125.D-3*(UELEM(2)+UELEM(1)) 204 : C 205 : Y(1)=125.D-3*(UELEM(6)-UELEM(4)) 206 : Y(2)=125.D-3*(UELEM(3)-UELEM(5)) 207 : Y(3)=125.D-3*(UELEM(1)-UELEM(2)) 208 : C 209 : C SI LES FLUX NE SECOMPENSENT PAS DEUX A DEUX 210 : C 211 : IF(ABS(X(1))/Q.GT.1.D-11)THEN 212 : IF(ABS(X(2))/Q.GT.1.D-11)THEN 213 : IF(ABS(X(3))/Q.GT.1.D-11)THEN 214 : C 215 : C ON TESTE TOUTES LES FACES DE SORTIE 216 : C 217 : DO 17 I=1,IDIM 218 : DO 18 I1=1,2 219 : C ON TESTE TOUTES LES FACES QUI ONT UN FLUX POSITIF 220 : IF(UELEM(IFAC8(2*I+I1-2)).GT.0.D0)THEN 221 : VAR1=ZI(I1)-(Y(I)/X(I)) 222 : VAR2=XDEP(I)-(Y(I)/X(I)) 223 : C ON ELIMINE LES CAS OÙ LE LOG EST NEG OU NON DEFINI 224 : IF((VAR1/VAR2).LE.0.D0)GOTO 6 225 : DTI=(1.D0/X(I))*LOG(VAR1/VAR2) 226 : IF(ABS(LOG(VAR1/VAR2)).LT.1.D-15)DTI=0.D0 227 : C ON GARDE LE PLUS PETIT TEMPS POSITIF 228 : IF(TEST(DTINT,DTI))THEN 229 : DTINT=DTI 230 : I2=I 231 : I3=I1 232 : I4=IJT4(1,I) 233 : I5=IJT4(2,I) 234 : ENDIF 235 : ENDIF 236 : 6 CONTINUE 237 : 18 CONTINUE 238 : 17 CONTINUE 239 : C ON GARDE LES COORDONNEES CORRESPONDANT AU TEMPS LE PLUS PETIT 240 : XARI(I2)=ZI(I3) 241 : XARI(I4)=(XDEP(I4)*EXP(X(I4)*DTINT))+ 242 : * ((Y(I4)/X(I4))*(1-EXP(X(I4)*DTINT))) 243 : XARI(I5)=(XDEP(I5)*EXP(X(I5)*DTINT))+ 244 : * ((Y(I5)/X(I5))*(1-EXP(X(I5)*DTINT))) 245 : ICONT=IFAC8((2*I2)+I3-2) 246 : 247 : ELSE 248 : C ON TESTE TOUTES LES FACES MAIS L'EQUATION SUIVANT Z DIFFERE 249 : C (Q1+Q2=O) 250 : DO 21 I=1,2 251 : DO 19 I1=1,2 252 : IF(UELEM(IFAC8(2*I+I1-2)).GT.0.D0)THEN 253 : VAR1=ZI(I1)-(Y(I)/X(I)) 254 : VAR2=XDEP(I)-(Y(I)/X(I)) 255 : C ON ELIMINE LES CAS OÙ LE LOG EST NEG OU NON DEFINI 256 : IF((VAR1/VAR2).LE.0.D0)GOTO 9 257 : DTI=(1.D0/X(I))*LOG(VAR1/VAR2) 258 : IF(ABS(LOG(VAR1/VAR2)).LT.1.D-15)DTI=0.D0 259 : c ON GARDE LE PLUS PETIT TEMPS POSITIF 260 : IF(TEST(DTINT,DTI))THEN 261 : DTINT=DTI 262 : I2=I 263 : I3=I1 264 : I4=IJT4(1,I) 265 : ENDIF 266 : ENDIF 267 : 9 CONTINUE 268 : 19 CONTINUE 269 : 21 CONTINUE 270 : C ON GARDE LES COORDONNEES ET LA FACE CORRESPONDANT AU TEMPS 271 : XARI(I2)=ZI(I3) 272 : XARI(I4)=(XDEP(I4)*EXP(X(I4)*DTINT))+ 273 : * ((Y(I4)/X(I4))*(1-EXP(X(I4)*DTINT))) 274 : XARI(3)=XDEP(3)-(Y(3)*DTINT) 275 : ICONT=IFAC8((2*I2)+I3-2) 276 : C ON TESTE MAINTENANT LES FACES Z=1 ET Z=-1 277 : IF(UELEM(2).GT.0.D0)THEN 278 : DTI=(XDEP(3)-1.D0)/Y(3) 279 : INT=2 280 : ELSE 281 : DTI=(XDEP(3)+1.D0)/Y(3) 282 : INT=1 283 : ENDIF 284 : C ON TESTE LE TEMPS 285 : IF(TEST(DTINT,DTI))THEN 286 : DTINT=DTI 287 : XARI(1)=(XDEP(1)*EXP(X(1)*DTINT))+ 288 : * ((Y(1)/X(1))*(1-EXP(X(1)*DTINT))) 289 : XARI(2)=(XDEP(2)*EXP(X(2)*DTINT))+ 290 : * ((Y(2)/X(2))*(1-EXP(X(2)*DTINT))) 291 : XARI(3)=ZI(INT) 292 : ICONT=INT 293 : ENDIF 294 : ENDIF 295 : ELSE 296 : C 297 : C ON A OBLIGATOIREMENT Q1+Q2 NON NUL DONC ON PEUT RAISONNER 298 : C COMME PRECEDEMMENT 299 : C ON TESTE TOUTES LES FACES MAIS L'EQUATION SUIVANT Y DIFFERE 300 : C (Q3+Q5=0) 301 : DO 23 I=1,2 302 : DO 22 I1=1,2 303 : C 304 : IF(UELEM(IFAC6(2*I+I1-2)).GT.0.D0)THEN 305 : VAR1=ZI(I1)-(Y(I)/X(I)) 306 : VAR2=XDEP(I)-(Y(I)/X(I)) 307 : C 308 : IF((VAR1/VAR2).LE.0.D0)GOTO 11 309 : DTI=(1.D0/X(I))*LOG(VAR1/VAR2) 310 : C 311 : IF(ABS(LOG(VAR1/VAR2)).LT.1.D-15)DTI=0.D0 312 : C 313 : IF(TEST(DTINT,DTI))THEN 314 : DTINT=DTI 315 : I2=IFAQ4(I) 316 : I3=I1 317 : I4=IFAT3(I) 318 : ENDIF 319 : ENDIF 320 : 11 CONTINUE 321 : 22 CONTINUE 322 : 23 CONTINUE 323 : XARI(I2)=ZI(I3) 324 : XARI(I4)=(XDEP(I4)*EXP(X(I4)*DTINT))+ 325 : * ((Y(I4)/X(I4))*(1-EXP(X(I4)*DTINT))) 326 : XARI(2)=XDEP(2)-(Y(2)*DTINT) 327 : ICONT=IFAC6((2*I2)+I3-2) 328 : C 329 : C ON TESTE MAINTENANT LES FACES Y=1 ET Y=-1 330 : C 331 : IF(UELEM(5).GT.0.D0)THEN 332 : DTI=(XDEP(2)-1.D0)/Y(2) 333 : INT=5 334 : ELSE 335 : DTI=(XDEP(2)+1.D0)/Y(2) 336 : INT=3 337 : ENDIF 338 : IF(TEST(DTINT,DTI))THEN 339 : DTINT=DTI 340 : XARI(1)=(XDEP(1)*EXP(X(1)*DTINT))+ 341 : * ((Y(1)/X(1))*(1-EXP(X(1)*DTINT))) 342 : XARI(3)=(XDEP(3)*EXP(X(3)*DTINT))+ 343 : * ((Y(3)/X(3))*(1-EXP(X(3)*DTINT))) 344 : XARI(2)=INT-4 345 : ICONT=INT 346 : ENDIF 347 : ENDIF 348 : C 349 : ELSE 350 : C ON TESTE TOUTE LES FACES MAIS L'EQUATION DE X DIFFERE 351 : C (Q4+Q6 =0) 352 : IF(ABS(X(2))/Q.GT.1.D-11)THEN 353 : DO 24 I=2,3 354 : DO 25 I1=1,2 355 : C 356 : IF(UELEM(IFAC8(2*I+I1-2)).GT.0.D0)THEN 357 : C 358 : VAR1=ZI(I1)-(Y(I)/X(I)) 359 : VAR2=XDEP(I)-(Y(I)/X(I)) 360 : C 361 : IF((VAR1/VAR2).LE.0.D0)GOTO 13 362 : C 363 : DTI=(1.D0/X(I))*LOG(VAR1/VAR2) 364 : IF(ABS(LOG(VAR1/VAR2)).LT.1.D-15)DTI=0.D0 365 : C 366 : IF(TEST(DTINT,DTI))THEN 367 : DTINT=DTI 368 : I2=I 369 : I3=I1 370 : I4=IJT4(2,I) 371 : ENDIF 372 : ENDIF 373 : 13 CONTINUE 374 : 25 CONTINUE 375 : 24 CONTINUE 376 : XARI(I2)=ZI(I3) 377 : XARI(I4)=(XDEP(I4)*EXP(X(I4)*DTINT))+ 378 : * ((Y(I4)/X(I4))*(1-EXP(X(I4)*DTINT))) 379 : XARI(1)=XDEP(1)-(Y(1)*DTINT) 380 : ICONT=IFAC8((2*I2)+I3-2) 381 : C ON TESTE MAINTENANT LES FACES X=1 ET X=-1 382 : IF(UELEM(6).GT.0.D0)THEN 383 : DTI=(XDEP(1)+1.D0)/Y(1) 384 : INT=6 385 : ELSE 386 : DTI=(XDEP(1)-1.D0)/Y(1) 387 : INT=4 388 : ENDIF 389 : IF(TEST(DTINT,DTI))THEN 390 : DTINT=DTI 391 : XARI(3)=(XDEP(3)*EXP(X(3)*DTINT))+ 392 : * ((Y(3)/X(3))*(1-EXP(X(3)*DTINT))) 393 : XARI(2)=(XDEP(2)*EXP(X(2)*DTINT))+ 394 : * ((Y(2)/X(2))*(1-EXP(X(2)*DTINT))) 395 : XARI(1)=ZI(INT-5) 396 : ICONT=INT 397 : ENDIF 398 : ELSE 399 : C 400 : C Q3+Q5 NEGLIGEABLE (DONC Q1+Q2 NEG) 401 : C ON CALCULE POUR TOUTES LES FACES 402 : DO 31 I=1,3 403 : DO 32 I1=1,2 404 : IF(UELEM(IFAC8(2*I+I1-2)).GT.0.D0)THEN 405 : C 406 : DTI=(XDEP(I)-ZI(I1))/Y(I) 407 : C 408 : IF(TEST(DTINT,DTI))THEN 409 : DTINT=DTI 410 : I2=I 411 : I3=I1 412 : I4=IJT4(1,I) 413 : I5=IJT4(2,I) 414 : ENDIF 415 : ENDIF 416 : 32 CONTINUE 417 : 31 CONTINUE 418 : XARI(I2)=ZI(I3) 419 : XARI(I4)=XDEP(I4)-(Y(I4)*DTINT) 420 : XARI(I5)=XDEP(I5)-(Y(I5)*DTINT) 421 : ICONT=IFAC8((2*I2)+I3-2) 422 : ENDIF 423 : ENDIF 424 : C 425 : C ON CHERCHE SI ON A AFFAIRE A UN NOEUD DU MAILLAGE 426 : C 427 : DO 45 I=1,IDIM 428 : IF(ABS(1.D0+XARI(I)).LT.EPS)ITEST=ITEST+1 429 : IF(ABS(1.D0-XARI(I)).LT.EPS)ITEST=ITEST+1 430 : 45 CONTINUE 431 : CC 432 : CC 433 : CC*** PRISME PRI6 434 : CC 435 : ELSEIF(ITYEL.EQ.16)THEN 436 : IDIM=3 437 : Q=MAX(ABS(UELEM(1)),ABS(UELEM(2)),ABS(UELEM(3)), 438 : * ABS(UELEM(4)),ABS(UELEM(5))) 439 : VAR=(UELEM(3)+UELEM(4)+UELEM(5))/2.D0 440 : X(1)=UELEM(5)/2.D0 441 : X(2)=UELEM(3)/2.D0 442 : X(3)=(UELEM(2)-UELEM(1))/2.D0 443 : C 444 : IF(ABS(VAR)/Q.GT.1.D-10)THEN 445 : C ON TESTE TOUTES LES FACES QUI ONT UN FLUX POSITIF 446 : C ET ON GARDE LE PLUS PETIT TEMPS POSITIF 447 : C ON TESTE LES FACES TRIANGULAIRES: FACE 4 448 : IF(UELEM(4).GT.0.D0)THEN 449 : A1=((1.D0-((X(1)+X(2))/VAR))) 450 : A2=(XDEP(1)+XDEP(2)-((X(1)+X(2))/VAR)) 451 : IF (A1/A2.LE.0.D0)GOTO 15 452 : DTINT=(1.D0/VAR)*LOG(A1/A2) 453 : XARI(2)=(XDEP(2)*EXP(VAR*DTINT))+((X(1)/VAR)* 454 : * (1.D0-EXP(VAR*DTINT))) 455 : XARI(3)=(XDEP(3)*EXP(-1.D0*VAR*DTINT))+((X(3)/VAR)* 456 : * (1.D0-EXP(-1.D0*VAR*DTINT))) 457 : XARI(1)=1.D0-XARI(2) 458 : ICONT=4 459 : 15 CONTINUE 460 : ENDIF 461 : C 462 : DO 8 I=1,IDIM-1 463 : C FACE 3 ET FACE 5 464 : IF(X(I).GT.0.D0)THEN 465 : A3=(-1.D0*X(I)/VAR) 466 : A4=(XDEP(I)-(X(I)/VAR)) 467 : IF(A3/A4.LE.0.D0)GOTO 16 468 : DTI=(1.D0/VAR)*LOG(A3/A4) 469 : C 470 : C ON NE GARDE QUE LE PLUS PETIT TEMPS POSITIF 471 : C 472 : IF(TEST(DTINT,DTI))THEN 473 : DTINT=DTI 474 : XARI(IJ(I))=(XDEP(IJ(I))*EXP(VAR*DTINT))+ 475 : * ((X(IJ(I))/VAR)*(1.D0-EXP(VAR*DTINT))) 476 : XARI(3)=(XDEP(3)*EXP(VAR*DTINT))+ 477 : * ((X(3)/VAR)*(1.D0-EXP(-1.D0*VAR*DTINT))) 478 : XARI(I)=0.D0 479 : ICONT=IFAT3(I)+2 480 : ENDIF 481 : 16 CONTINUE 482 : ENDIF 483 : C ON TESTE LES FACES OPPOSEES 484 : IF(UELEM(I).GT.0.D0)THEN 485 : A5=(ZI(I)-(X(3)/VAR)) 486 : A6=(XDEP(3)-(X(3)/VAR)) 487 : IF(A5/A6.LE.0.D0)GOTO 26 488 : DTI=(-5.D-1/VAR)*LOG(A5/A6) 489 : C 490 : C ON NE GARDE QUE LE PLUS PETIT TEMPS POSITIF 491 : C 492 : IF(TEST(DTINT,DTI))THEN 493 : DTINT=DTI 494 : XARI(I)=(XDEP(I)*EXP(VAR*DTINT))+ 495 : * ((X(I)/VAR)*(1.D0-EXP(VAR*DTINT))) 496 : XARI(IJ(I))=(XDEP(IJ(I))*EXP(VAR*DTINT))+ 497 : * ((X(IJ(I))/VAR)*(1.D0-EXP(VAR*DTINT))) 498 : XARI(3)=ZI(I) 499 : ICONT=I 500 : ENDIF 501 : 26 CONTINUE 502 : ENDIF 503 : 8 CONTINUE 504 : C 505 : ELSE 506 : C 507 : C Q3 NON NEGLIGEABLE ON VA TESTER LES FACES TRIANGULAIRES 508 : C FACE 4 509 : IF(UELEM(4).GT.0.D0)THEN 510 : DTINT=(2.D0/UELEM(4))*(1.D0-XDEP(1)-XDEP(2)) 511 : XARI(2)=XDEP(2)-(X(2)*DTINT) 512 : XARI(1)=1.D0-XARI(2) 513 : ICONT=4 514 : ENDIF 515 : C 516 : DO 12 I=1,IDIM-1 517 : IF(X(I).GT.0.D0)THEN 518 : DTI=XDEP(I)/X(I) 519 : C ON NE GARDE QUE LE PLUS PETIT TEMPS POSITIF 520 : IF(TEST(DTINT,DTI))THEN 521 : DTINT=DTI 522 : XARI(IJ(I))=XDEP(IJ(I))-(X(IJ(I))*DTINT) 523 : XARI(I)=0.D0 524 : ICONT=IFAT3(I)+2 525 : ENDIF 526 : ENDIF 527 : 12 CONTINUE 528 : XARI(3)=XDEP(3)+(2.D0*X(3)*DTINT) 529 : C 530 : C SI Q1 NON NEGLIGEABLE ON TESTE LES FACES OPPOSEES A FLUX POSITIF 531 : IF(UELEM(1).GT.0.D0)THEN 532 : DTI=(-1.D0-XDEP(3))/(2.D0*X(3)) 533 : IND=1 534 : ELSE 535 : DTI=(1.D0-XDEP(3))/(2.D0*X(3)) 536 : IND=2 537 : ENDIF 538 : IF(TEST(DTINT,DTI))THEN 539 : DTINT=DTI 540 : DO 14 I1=1,IDIM-1 541 : XARI(I1)=XDEP(I1)-(X(I1)*DTINT) 542 : 14 CONTINUE 543 : XARI(3)=ZI(IND) 544 : ICONT=IND 545 : ENDIF 546 : ENDIF 547 : C 548 : C ON CHERCHE SI ON A AFFAIRE A UN NOEUD DU MAILLAGE 549 : C 550 : IF((ABS(XARI(3)-1.D0).LT.EPS).OR. 551 : * (ABS(XARI(3)+1.D0).LT.EPS))ITEST=ITEST+1 552 : DO 41 I=1,2 553 : IF(ABS(XARI(I)).LT.EPS)ITEST=ITEST+1 554 : IF(ABS(1.D0-XARI(I)).LT.EPS)ITEST=ITEST+1 555 : 41 CONTINUE 556 : C 557 : C 558 : C ****************** TETRAEDRE 559 : C * FACE 1 Y=0 560 : C * FACE 2 Z=0 561 : C * FACE 4 X=0 562 : C * FACE 3 1-X-Y-Z=0 563 : ELSEIF(ITYEL.EQ.23)THEN 564 : IDIM=3 565 : C*** FACES 1 2 ET 4 566 : C 567 : C** ON TESTE L'INTERSECTION AVEC LA FACE 3 ET TOUTES LES 568 : C FACES AYANT UN FLUX POSITIF 569 : C 570 : IF(UELEM(3).GT.0.D0)THEN 571 : DTINT=(1-XDEP(3)-XDEP(2)-XDEP(1))/(2.D0*UELEM(3)) 572 : DO 56 I=1,IDIM 573 : XARI(I)=XDEP(I)-(2.D0*UELEM(IFAT4(I))*DTINT) 574 : 56 CONTINUE 575 : ICONT=3 576 : ENDIF 577 : C 578 : DO 55 I=1,IDIM 579 : C 580 : IF(UELEM(IFAT4(I)).GT.0.D0)THEN 581 : DTI=XDEP(I)/(2.D0*UELEM(IFAT4(I))) 582 : C 583 : C ON NE GARDE QUE LE PLUS PETIT TEMPS POSITIF 584 : C 585 : IF(TEST(DTINT,DTI))THEN 586 : DTINT=DTI 587 : I2=IJT4(1,I) 588 : I3=IJT4(2,I) 589 : C ON GARDE LES COORDONNES CORRESPONDANT AUPLUS PETIT TEMPS 590 : XARI(I2)=XDEP(I2)-(2.D0*UELEM(IFAT4(I2))*DTINT) 591 : XARI(I3)=XDEP(I3)-(2.D0*UELEM(IFAT4(I3))*DTINT) 592 : XARI(I)=0.D0 593 : ICONT=IFAT4(I) 594 : ENDIF 595 : ENDIF 596 : 55 CONTINUE 597 : C 598 : 599 : C ON CHERCHE LE NUMERO DU NOEUD OU DE L'ARETE 600 : C 601 : DO 51 I=1,IDIM 602 : IF(ABS(XARI(I)).LT.EPS)ITEST=ITEST+1 603 : IF(ABS(1.D0-XARI(I)).LT.EPS)ITEST=ITEST+1 604 : 605 : 51 CONTINUE 606 : ENDIF 607 : C 608 : C TRAITEMENT DES CAS PARTICULIERS (NOEUDS/ARETES) 609 : C 610 : 611 : C SI LE POINT DE SORTIE EST UN NOEUD 612 : IF(ITEST.EQ.IDIM)THEN 612 : CALL NUMNOE(XARI,ITYEL,INOEU) 613 : ELSE 614 : C SINON ON VERIFIE EN 3D QU'ON TOMBE PAS SUR UNE ARETE 615 : IF(IDIM.EQ.3)THEN 615 : CALL ARTREF(XARI,ITYEL,IART) 616 : ENDIF 617 : ENDIF 618 : END 619 : 620 : 621 : 622 : 623 :
© Cast3M 2003 - All rights reserved.
Disclaimer