1 : C TRJINT SOURCE CHAT 05/01/13 03:49:59 5004 2 : SUBROUTINE TRJINT(XARI,XDEP,XINT,UELEM,DTCA,DTINT, 3 : * IDIM,ICONT,ITYEL,IO) 4 : C 5 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 6 : C 7 : C INTERSECTION TRAJECTOIRE-ELEMENT ( issu de TRIOEF) 8 : C 9 : C ENTREES 10 : C XDEP POSITION INITIALE 11 : C XARI POSITION FINALE 12 : C UELEM VITESSE 13 : C DTCA PAS DE TEMPS 14 : C IDIM DIMENSION DE L ESPACE 15 : C ITYEL TYPE DE L ELEMENT 16 : C 17 : C SORTIES 18 : C XINT POSITION A L INTERSECTION 19 : C DTINT TEMPS ECOULE JUSQU A LA SORTIE 20 : C ICONT NUMERO DE LA PORTE DE SORTIE 21 : C IO DIFFERENT DE 0 ON EST SORTI DE L ELEMENT MAIS ON NE SAIT 22 : C PAS PAR QUEL COTE 23 : C 24 : C TYPES D ELEMENTS CONSIDERES 25 : C TRI3 TRI6 TRI7 QUA4 QUA9 CUB8 PRI6 26 : C 4 6 7 8 11 14 16 27 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 28 : C 29 : IMPLICIT INTEGER(I-N) 30 : IMPLICIT REAL*8 (A-H,O-Z) 31 : C 32 : DIMENSION IJ(2),ZI(2),IFAT3(3),IFAQ4(4),IFAC8(6) 33 : DIMENSION IJKL(3,3),IJT4(2,3),IFAT4(4) 34 : DIMENSION XARI(3),XDEP(3),XINT(3),UELEM(3),X(4) 35 : C 36 : DATA IJKL/1,2,3,2,1,3,3,2,1/ 37 : DATA IJ /2,1/ 38 : DATA IJT4 /2,3, 1,3, 1,2/ 39 : DATA IFAT3 /3,1,2/ 40 : DATA IFAQ4 /4,2,1,3/ 41 : DATA IFAC8 /6,4,3,5,1,2/ 42 : DATA IFAT4 /4,1,2,3/ 43 : DATA ZI/-1.D0,1.D0/ 44 : C 45 : ICONT=0 46 : IO=0 47 : DTT=DTCA 48 : DO 50 I=1,IDIM 49 : XINT(I)=XDEP(I) 50 : 50 CONTINUE 51 : C write(6,*)' inttrj ',XDEP,XARI,UELEM,DTCA,IDIM,ITYEL 52 : C 53 : C*** TRIANGLES TRI3 TRI6 TRI7 * 54 : C * * 55 : C FACES 3 * * 2 56 : C * * 57 : C ****** 58 : C 1 59 : IF(ITYEL.EQ.4.OR.ITYEL.EQ.6.OR.ITYEL.EQ.7)THEN 60 : IDIM=2 61 : XARI(3)=1.D0-XARI(1)-XARI(2) 62 : ITEST=0 63 : DO 1 I=1,3 64 : IF(XARI(I)*(1.D0-XARI(I)).GE.0.D0)ITEST=ITEST+1 65 : 1 CONTINUE 66 : IF(ITEST.NE.3)THEN 67 : C ON SORT DU TRIANGLE 68 : C 69 : C*** COTE 1 ET 3 70 : C 71 : DO 2 I1=1,IDIM 72 : IF(UELEM(I1).NE.0.AND.XDEP(I1).GT.0.D0)THEN 73 : DTI=-XDEP(I1)/UELEM(I1) 74 : DTINT=DTI 75 : IF(DTI*(DTT-DTI).GE.0.D0)THEN 76 : ITEST=0 77 : I3=IJ(I1) 78 : XINT(I3)=XDEP(I3)+UELEM(I3)*DTI 79 : X3=1.D0-XINT(I3) 80 : IF(XINT(I3)*(1.D0-XINT(I3)).GE.0.D0)ITEST=ITEST+1 81 : IF(X3*(1.D0-X3).GE.0.D0)ITEST=ITEST+1 82 : IF(ITEST.EQ.2)THEN 83 : ICONT=IFAT3(I1) 84 : XINT(I1)=0.D0 85 : DTINT=DTI 86 : RETURN 87 : ENDIF 88 : ENDIF 89 : ENDIF 90 : 2 CONTINUE 91 : C 92 : C*** COTE 2 93 : C 94 : IF((UELEM(1)+UELEM(2)).NE.0.D0.AND. 95 : * (1.D0-XDEP(1)-XDEP(2)).GT.0.D0) 96 : * THEN 97 : DTI=(1.D0-XDEP(1)-XDEP(2))/(UELEM(1)+UELEM(2)) 98 : DTINT=DTI 99 : IF(DTI*(DTT-DTI).GE.0.D0)THEN 100 : XINT(1)=XDEP(1)+UELEM(1)*DTI 101 : XINT(2)=XDEP(2)+UELEM(2)*DTI 102 : ITEST=0 103 : IF(XINT(1)*(1.D0-XINT(1)).GE.0.D0)ITEST=ITEST+1 104 : IF(XINT(2)*(1.D0-XINT(2)).GE.0.D0)ITEST=ITEST+1 105 : IF(ITEST.EQ.2)THEN 106 : ICONT=2 107 : DTINT=DTI 108 : RETURN 109 : ENDIF 110 : ENDIF 111 : ENDIF 112 : C 113 : C*** ?????????? 114 : C 115 : IO=3 116 : ENDIF 117 : C 118 : C*** QUADRANGLE QUA4 QUA9 3 119 : C ***** 120 : C * * 121 : C FACES 4 * * 2 122 : C ***** 123 : C 1 124 : ELSEIF(ITYEL.EQ.8.OR.ITYEL.EQ.11)THEN 125 : IDIM=2 126 : ITEST=0 127 : DO 3 I=1,IDIM 128 : XY=(1.D0+XARI(I))*(1.D0-XARI(I)) 129 : IF(XY.GE.0.D0)ITEST=ITEST+1 130 : 3 CONTINUE 131 : IF(ITEST.NE.IDIM)THEN 132 : ID=IDIM-1 133 : DO 4 I1=1,IDIM 134 : IF(UELEM(I1).NE.0.D0)THEN 135 : DO 5 I4=1,2 136 : C IF(ZI(I4)*(ZI(I4)-XDEP(I1)).le.0.D0) 137 : C * write(6,*)' xdep i1',i1,i4,xdep(i1) 138 : IF(ZI(I4)*(ZI(I4)-XDEP(I1)).GT.0.D0)THEN 139 : DTI=(ZI(I4)-XDEP(I1))/UELEM(I1) 140 : DTINT=DTI 141 : JI=I1 142 : C IF(DTI*(DTT-DTI).lt.0.D0) 143 : C * write(6,*)' dti dtt',dti,dtt 144 : IF(DTI*(DTT-DTI).GE.0.D0)THEN 145 : ITEST=0 146 : DO 6 I5=2,IDIM 147 : I7=IJKL(I5,I1) 148 : XINT(I7)=XDEP(I7)+UELEM(I7)*DTI 149 : XXX=(XINT(I7)+1.D0)*(1.D0-XINT(I7)) 150 : IF(XXX.GE.0.D0)ITEST=ITEST+1 151 : C write(6,*)' itest id ',itest,id,xxx 152 : 6 CONTINUE 153 : IF(ITEST.EQ.ID)THEN 154 : DTINT=DTI 155 : XINT(I1)=ZI(I4) 156 : ICONT=IFAQ4(2*I1+I4-2) 157 : C WRITE(6,*)' SORTIE TRJINT ',XINT,DTINT,ICONT,IO 158 : RETURN 159 : ENDIF 160 : ENDIF 161 : ENDIF 162 : 5 CONTINUE 163 : ENDIF 164 : 4 CONTINUE 165 : C 166 : C*** ?????????? 167 : C 168 : IO=4 169 : ENDIF 170 : C 171 : C*** CUBE CUB8 172 : C 173 : ELSEIF(ITYEL.EQ.14)THEN 174 : IDIM=3 175 : ITEST=0 176 : DO 13 I=1,IDIM 177 : XY=(1.D0+XARI(I))*(1.D0-XARI(I)) 178 : IF(XY.GE.0.D0)ITEST=ITEST+1 179 : 13 CONTINUE 180 : IF(ITEST.NE.IDIM)THEN 181 : ID=IDIM-1 182 : DO 14 I1=1,IDIM 183 : IF(UELEM(I1).NE.0.D0)THEN 184 : DO 15 I4=1,2 185 : IF(ZI(I4)*(ZI(I4)-XDEP(I1)).GT.0.D0)THEN 186 : DTI=(ZI(I4)-XDEP(I1))/UELEM(I1) 187 : DTINT=DTI 188 : JI=I1 189 : IF(DTI*(DTT-DTI).GE.0.D0)THEN 190 : ITEST=0 191 : DO 16 I5=2,IDIM 192 : I7=IJKL(I5,I1) 193 : XINT(I7)=XDEP(I7)+UELEM(I7)*DTI 194 : XXX=(XINT(I7)+1.D0)*(1.D0-XINT(I7)) 195 : IF(XXX.GE.0)ITEST=ITEST+1 196 : 16 CONTINUE 197 : IF(ITEST.EQ.ID)THEN 198 : DTINT=DTI 199 : XINT(I1)=ZI(I4) 200 : ICONT=IFAC8(2*I1+I4-2) 201 : RETURN 202 : ENDIF 203 : ENDIF 204 : ENDIF 205 : 15 CONTINUE 206 : ENDIF 207 : 14 CONTINUE 208 : C 209 : C*** ?????????? 210 : C 211 : IO=4 212 : ENDIF 213 : C 214 : C*** PRISME PRI6 215 : C 216 : ELSEIF(ITYEL.EQ.16)THEN 217 : IDIM=3 218 : ITEST=0 219 : X(1)=XARI(1) 220 : X(2)=XARI(2) 221 : X(3)=1.D0-X(1)-X(2) 222 : X(4)=XARI(3) 223 : DO 7 I=1,3 224 : IF(X(I)*(1.D0-X(I)).GE.0.D0)ITEST=ITEST+1 225 : 7 CONTINUE 226 : IF((X(4)+1.D0)*(1.D0-X(4)).GE.0.D0)ITEST=ITEST+1 227 : IF(ITEST.NE.4)THEN 228 : C 229 : C*** FACE 5 ET 3 230 : C 231 : DO 8 I1=1,2 232 : IF(UELEM(I1).NE.0.D0.AND.XDEP(I1).GT.0.D0)THEN 233 : DTI=-XDEP(I1)/UELEM(I1) 234 : DTINT=DTI 235 : IF(DTI*(DTT-DTI).GE.0.D0)THEN 236 : ITEST=0 237 : I3=IJ(I1) 238 : XINT(I3)=XDEP(I3)+UELEM(I3)*DTI 239 : X3=1.D0-XINT(I3) 240 : XINT(3)=XDEP(3)+UELEM(3)*DTI 241 : IF(XINT(I3)*(1.D0-XINT(I3)).GE.0.D0)ITEST=ITEST+1 242 : IF((XINT(3)+1.D0)*(1.D0-XINT(3)).GE.0.D0)ITEST=ITEST+1 243 : IF(X3*(1.D0-X3).GE.0.D0)ITEST=ITEST+1 244 : IF(ITEST.EQ.3)THEN 245 : XINT(I1)=0.D0 246 : ICONT=IFAT3(I1)+2 247 : DTINT=DTI 248 : RETURN 249 : ENDIF 250 : ENDIF 251 : ENDIF 252 : 8 CONTINUE 253 : C 254 : C*** FACE 4 255 : C 256 : IF((UELEM(1)+UELEM(2)).NE.0.D0.AND. 257 : * (1.D0-XDEP(1)-XDEP(2)).GT.0.D0) 258 : * THEN 259 : DTI=(1.D0-XDEP(1)-XDEP(2))/(UELEM(1)+UELEM(2)) 260 : DTINT=DTI 261 : IF(DTI*(DTT-DTI).GE.0.D0)THEN 262 : XINT(1)=XDEP(1)+UELEM(1)*DTI 263 : XINT(2)=XDEP(2)+UELEM(2)*DTI 264 : XINT(3)=XDEP(3)+UELEM(3)*DTI 265 : ITEST=0 266 : IF(XINT(1)*(1.D0-XINT(1)).GE.0.D0)ITEST=ITEST+1 267 : IF(XINT(2)*(1.D0-XINT(2)).GE.0.D0)ITEST=ITEST+1 268 : IF((1.D0-XINT(3))*(XINT(3)+1.D0).GE.0.D0)ITEST=ITEST+1 269 : IF(ITEST.EQ.3)THEN 270 : ICONT=4 271 : DTINT=DTI 272 : RETURN 273 : ENDIF 274 : ENDIF 275 : ENDIF 276 : C 277 : C*** FACE 1 ET 2 278 : C 279 : DO 9 I2=1,2 280 : IF(UELEM(3).NE.0.D0.AND.ZI(I2)*(ZI(I2)-XDEP(3)) 281 : * .GT.0.D0)THEN 282 : DTI=(ZI(I2)-XDEP(3))/UELEM(3) 283 : DTINT=DTI 284 : IF(DTI*(DTT-DTI).GE.0.D0)THEN 285 : ITEST=0 286 : DO 10 I5=1,2 287 : XINT(I5)=XDEP(I5)+UELEM(I5)*DTI 288 : IF((1.D0-XINT(I5))*(XINT(I5)+1.D0).GE.0.D0) 289 : * ITEST=ITEST+1 290 : 10 CONTINUE 291 : X3=1.D0-XINT(1)-XINT(2) 292 : IF(X3*(1.D0-X3).GE.0.D0)ITEST=ITEST+1 293 : IF(ITEST.EQ.3)THEN 294 : XINT(3)=ZI(I2) 295 : ICONT=I2 296 : DTINT=DTI 297 : RETURN 298 : ENDIF 299 : ENDIF 300 : ENDIF 301 : 9 CONTINUE 302 : C 303 : C*** ????????????? 304 : C 305 : IO=5 306 : ENDIF 307 : C 308 : C ****************** TETRAEDRE 309 : C * FACE 1 Y=0 310 : C * FACE 2 Z=0 311 : C * FACE 4 X=0 312 : C * FACE 3 1-X-Y-Z=0 313 : ELSEIF(ITYEL.EQ.23)THEN 314 : ITEST=0 315 : DO 20 I=1,3 316 : IF(XARI(I)*(1.D0-XARI(I)).GE.0.D0)ITEST=ITEST+1 317 : 20 CONTINUE 318 : XXX=XARI(1)+XARI(2)+XARI(3) 319 : IF(XXX*(1.D0-XXX).GE.0.D0)ITEST=ITEST+1 320 : C write(6,*)' trjint ', XXX,XARI(1),XARI(2),XARI(3),itest 321 : IF(ITEST.NE.4)THEN 322 : C write(6,*)' trjint ', XXX,XARI(1),XARI(2),XARI(3),itest 323 : C ON SORT DU TETRAEDRE 324 : C 325 : C*** FACES 1 2 ET 4 326 : C 327 : DO 25 I1=1,IDIM 328 : C write(6,*)'axes',i1,uelem,xdep 329 : IF(UELEM(I1).NE.0.AND.XDEP(I1).GT.0.D0)THEN 330 : DTI=-XDEP(I1)/UELEM(I1) 331 : DTINT=DTI 332 : IF(DTI*(DTT-DTI).GE.0.D0)THEN 333 : ITEST=0 334 : I2=IJT4(1,I1) 335 : I3=IJT4(2,I1) 336 : XINT(I2)=XDEP(I2)+UELEM(I2)*DTI 337 : XINT(I3)=XDEP(I3)+UELEM(I3)*DTI 338 : X3=1.D0-XINT(I2)-XINT(I3) 339 : IF(XINT(I2)*(1.D0-XINT(I2)).GE.0.D0)ITEST=ITEST+1 340 : IF(XINT(I3)*(1.D0-XINT(I3)).GE.0.D0)ITEST=ITEST+1 341 : IF(X3*(1.D0-X3).GE.0.D0)ITEST=ITEST+1 342 : C write(6,*)'xint',xint,itest 343 : IF(ITEST.EQ.3)THEN 344 : ICONT=IFAT4(I1) 345 : XINT(I1)=0.D0 346 : DTINT=DTI 347 : C write(6,*)' sorti face ',icont,xint,io 348 : RETURN 349 : ENDIF 350 : ENDIF 351 : ENDIF 352 : 25 CONTINUE 353 : C 354 : C*** FACE 3 355 : C 356 : IF((UELEM(1)+UELEM(2)+UELEM(3)).NE.0.D0.AND. 357 : * (1.D0-XDEP(1)-XDEP(2)-XDEP(3)).GT.0.D0) 358 : * THEN 359 : DTI=(1.D0-XDEP(1)-XDEP(2)-XDEP(3))/ 360 : * (UELEM(1)+UELEM(2)+UELEM(3)) 361 : DTINT=DTI 362 : IF(DTI*(DTT-DTI).GE.0.D0)THEN 363 : XINT(1)=XDEP(1)+UELEM(1)*DTI 364 : XINT(2)=XDEP(2)+UELEM(2)*DTI 365 : XINT(3)=XDEP(3)+UELEM(3)*DTI 366 : ITEST=0 367 : IF(XINT(1)*(1.D0-XINT(1)).GE.0.D0)ITEST=ITEST+1 368 : IF(XINT(2)*(1.D0-XINT(2)).GE.0.D0)ITEST=ITEST+1 369 : IF(XINT(3)*(1.D0-XINT(3)).GE.0.D0)ITEST=ITEST+1 370 : IF(ITEST.EQ.3)THEN 371 : ICONT=3 372 : DTINT=DTI 373 : C XINT(3)=1.D0-XINT(1)-XINT(2) 374 : C write(6,*)' sorti face ',icont,xint,io 375 : RETURN 376 : ENDIF 377 : ENDIF 378 : ENDIF 379 : C 380 : C*** ?????????? 381 : C 382 : IO=3 383 : ENDIF 384 : ENDIF 385 : C 386 : END 387 : 388 :
© Cast3M 2003 - All rights reserved.
Disclaimer