1 : C TRJAVA SOURCE CHAT 05/01/13 03:48:49 5004 2 : SUBROUTINE TRJAVA(IZVIT,IZPART,IZN3,IPART,TMIN,TMAX,IZCOU, 3 : * MELEME,IELTFA,IZCENT,IFACEL,DELTAT,IZSH,IEROR,DT1) 4 : C 5 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 6 : C Issu du sp AVATAR de TRIO-EF ( adapté par F Auriol) 7 : C 8 : C FAIT AVANCER UNE PARTICULE ( COORDONNEES DE REFERENCES ) 9 : C JUSQU'A L'INSTANT TMAX OU A L'INSTANT DE SORTIE 10 : C IZVIT SEGMENT DECRIVANT LES VITESSES ( <--- TRJVIT TRJFLU) 11 : C IZPART POSITIONS INITIALES DES PARTICULES 12 : C IZN3 SEGMENT RESULTANT (AJUSTE ICI) 13 : C IPART NUMERO DE LA PARTICULE 14 : C TMIN INSTANT DE DEPART 15 : C TMAX INSTANT A NE PAS DEPASSER 16 : C IZCOU SEGMENT DES DT DONNANT UN NOMBRE DE COURANT DE 1 17 : C COU VALEUR DU NOMBRE DE COURANT VOULU 18 : C 2 PERM 19 : C 3 TRANS 20 : C MELEME POINTEUR DU MAILLAGE 21 : C IELTFA POINTEUR DE DOMAINE.ELTFA 22 : C IZCENT POINTEUR DE DOMAINE.CENTRE 23 : C IFACEL POINTEUR DE DOMAINE.FACEL 24 : C DELTAT PAS DE TEMPS AVEC LEQUEL ON CONSERVE LES RESULTATS 25 : C 26 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 27 : C 28 : IMPLICIT INTEGER(I-N) 29 : IMPLICIT REAL*8 (A-H,O-Z) 30 : C 31 : -INC CCOPTIO 32 : -INC SMCOORD 33 : -INC SMELEME 34 : -INC SMCHPOI 35 : C 36 : POINTEUR IZCENT.MELEME,IELTFA.MELEME,IZFAC1.MELEME 37 : POINTEUR IFACEL.MELEME 38 : C 39 : SEGMENT IZPART 40 : INTEGER NLEPA(NPART),NUMPA(NPART) 41 : REAL*8 COORPA(NDIM,NPART) 42 : ENDSEGMENT 43 : C 44 : SEGMENT IZN3 45 : INTEGER NAPAR3(NPOS),NUM3(NPOS) 46 : REAL*8 CREF3(NDIM,NPOS),TPAR(NPOS) 47 : ENDSEGMENT 48 : C 49 : SEGMENT IZCOU 50 : REAL*8 DTCO(NEL),COU 51 : ENDSEGMENT 52 : C 53 : SEGMENT IZSH 54 : REAL*8 SHP(6,MNO9),SHY(12,MNO9),XYZL(3,MNO9) 55 : ENDSEGMENT 56 : C 57 : SEGMENT IZVIT 58 : REAL*8 TEMTRA(NVIPT) 59 : INTEGER IPUN(NBS),IDUN(NBS),IPVPT(NVIPT),IFORML 60 : ENDSEGMENT 61 : C IDUN(I) nombre d elements avant le sous maillage I 62 : C IPVPT pointeurs de izvpt pour chaque pas de temps 63 : SEGMENT IZVPT 64 : INTEGER IPUN1(NBS),IPUMAX 65 : ENDSEGMENT 66 : C 67 : SEGMENT IZUN 68 : REAL*8 UN(I1,I2,I3) 69 : ENDSEGMENT 70 : POINTEUR IZUN1.IZUN,IZUN2.IZUN 71 : SEGMENT IZUMAX 72 : REAL*8 UMAX(NBREL) 73 : ENDSEGMENT 74 : C 75 : DIMENSION UELEM(3),XARI(3),XDEP(3),XINT(3) 76 : C 77 : C*** 78 : C 79 : ISENS=0 80 : IEROR=0 81 : KIO=0 82 : NDIM=IDIM 83 : NPOS=50 84 : SEGINI IZN3 85 : SEGACT IZVIT 86 : IFORMU=IFORML 87 : IEL1=NLEPA(IPART) 88 : DO 2 ID1=1,NDIM 89 : XDEP(ID1)=COORPA(ID1,IPART) 90 : 2 CONTINUE 91 : ITER=2 92 : NAPAR3(1)=IEL1 93 : DO 3 ID1=1,NDIM 94 : CREF3(ID1,1)=XDEP(ID1) 95 : 3 CONTINUE 96 : TCOUR=TMIN 97 : TPAR(1)=TMIN 98 : NUAPAR=IEL1 99 : C DT1=0.D0 100 : SENS=SIGN(1.D0,(TMAX-TMIN)) 101 : C 102 : SEGACT IELTFA,IZCENT,IFACEL 103 : SEGACT IZSH 104 : NVIPT=TEMTRA(/1) 105 : 1 CONTINUE 106 : C IVPT 1 EN PERMANENT 107 : C EN TRANSITOIRE NUMERO DU PAS DE TEMPS DIRECTEMENT 108 : C SUPERIEUR AU TEMPS CONSIDERE 109 : 110 : IF(NVIPT.EQ.1)THEN 111 : IVPT=1 112 : ELSE 113 : IVPT=2 113 : CALL TRJTPT(IZVIT,TCOUR,IVPT) 114 : ENDIF 114 : CALL MELNEL(IEL1,MELEME,IPT1,NEL0,0) 115 : NOEL1=IPT1.NUM(/1) 116 : IELL=IEL1-NEL0 116 : CALL DOXE(XCOOR,IDIM,NOEL1,IPT1.NUM,IELL,XYZL) 117 : ITY1=IPT1.ITYPEL 118 : C 119 : C*** NOMBRE DE COURANT 120 : C 121 : C IF(DTCO(IEL1).EQ.0.D0) THEN 122 : C WRITE(6,*) 'APPEL A TRJCN5 DANS TRJAVA ' 122 : CALL TRJCN5(ITY1,IZSH,IJK) 123 : C IF(IERR.GT.0) RETURN 124 : IF(ISENS.NE.IJK)THEN 125 : IF(ISENS.EQ.0)THEN 126 : ISENS=IJK 127 : ELSE 128 : IEROR=1 129 : GO TO 10 130 : ENDIF 131 : ENDIF 132 : IZVPT=IPVPT(IVPT) 133 : SEGACT IZVPT 134 : IZUMAX=IPUMAX 135 : SEGACT IZUMAX 136 : UEM=UMAX(IEL1) 137 : SEGDES IZUMAX 138 : SEGDES IZVPT 139 : IF(IVPT.NE.1)THEN 140 : IZVPT=IPVPT(IVPT-1) 141 : SEGACT IZVPT 142 : IZUMAX=IPUMAX 143 : SEGACT IZUMAX 144 : UEM=MAX(UEM,UMAX(IEL1)) 145 : SEGDES IZUMAX 146 : SEGDES IZVPT 147 : ENDIF 148 : NUCENT=IZCENT.NUM(1,IEL1) 148 : CALL MELNEL(IEL1,IELTFA,IZFAC1,NEL0,1) 148 : CALL TRJCOU(UEM,IZCOU,IZFAC1,IELL,IEL1,NUCENT) 149 : C ENDIF 150 : C 151 : C*** PAS DE TEMPS 152 : C 153 : TPREC=TCOUR 154 : IF(ABS(DELTAT).GT.1.D-15)THEN 155 : DTT=DELTAT 156 : ELSE 157 : DTT=10.D0*COU*DTCO(IEL1)*SENS 158 : IF (ABS(TMAX-TPREC).LE.ABS(DTT))DTT=TMAX-TPREC 159 : DT1=0.D0 160 : ENDIF 161 : DTI=DTCO(IEL1)*COU*SENS 162 : C write(6,*) ' dtco ',DTCO(IEL1),IEL1,COU,DTT 163 : IF(ABS(DTT-DT1).LT.ABS(DTI*1.1D0))THEN 164 : DTI=DTT-DT1 165 : NI=1 166 : ELSE 167 : NI=INT((DTT-DT1)/DTI)+1 168 : DTI=(DTT-DT1)/NI 169 : ENDIF 170 : C write(6,*)iel1,ni,dti,dt1,dtt,li 171 : LI=1 172 : C DT1=0.D0 173 : C 174 : 6 CONTINUE 175 : DTC=DTI*0.5D0 176 : DTC1=DTC 177 : C write(6,*)' trjava li ni ',li,ni,iel1,dti,dtc,dtt 178 : TTEMP=TPREC+(LI-1)*DTI 178 : IF(IVPT.NE.1)CALL TRJTPT(IZVIT,TTEMP,IVPT) 178 : CALL TRJVEL(IZVIT,IZUN,IEL1,IVPT,TTEMP) 179 : NOUN=UN(/2) 179 : CALL TRJULO(XYZL,UN(1,1,IELL),UELEM,XDEP,ITY1,NDIM,NOEL1,NOUN, 180 : * IFORMU,SHP,SHY) 181 : DO 9 ID1=1,NDIM 182 : XARI(ID1)=XDEP(ID1)+UELEM(ID1)*DTC 183 : 9 CONTINUE 184 : ITYG=NUMGEO(ITY1) 184 : CALL TRJINT(XARI,XDEP,XINT,UELEM,DTT,DTINT,NDIM,ICONT,ITYG,IO) 185 : IF(IO.NE.0)THEN 186 : C 187 : C*** RATTRAPAGE 188 : C 189 : KIO=KIO+1 190 : C IF(MOD(KIO,2).EQ.0)WRITE(6,999)IPART,IEL1,IEL2 191 : IF(MOD(KIO,2).EQ.0)THEN 192 : INTERR(1)=IPART 193 : INTERR(2)=IEL1 193 : CALL ERREUR(-298) 194 : ENDIF 195 : C WRITE(6,*)'TRJRAT 1 ',IPART,IEL1,IEL2, 196 : C * XARI,XDEP,XINT,UELEM,DTT,DTINT,NDIM,ICONT,ITYG,IO 197 : C IF(MOD(KIO,10).EQ.0) STOP 198 : 199 : INTERR(1)=IPART 199 : CALL TRJRAT(IEL1,IEL2,XDEP,XARI,IZCOU,MELEME,IZVIT,DTC1, 200 : * IZCENT,IELTFA,XINT,DTINT,ICONT,IZSH,TTEMP) 201 : C write(6,*) ' apres trjrat ', icont ,iel1,iel2 202 : IF(ICONT.NE.0)GO TO 16 203 : IF(IEL2.EQ.0)THEN 204 : DO 99 ID1=1,NDIM 205 : XINT(ID1)=XDEP(ID1) 206 : 99 CONTINUE 207 : DTINT=0.D0 208 : GO TO 10 209 : ENDIF 210 : DO 12 ID1=1,NDIM 211 : XDEP(ID1)=XARI(ID1) 212 : 12 CONTINUE 213 : TCOUR=TPREC+(LI-1)*DTI+DTC1 214 : NUAPAR=IEL2 215 : DT1=0.D0 216 : GO TO 13 217 : ENDIF 218 : 16 CONTINUE 219 : IF(ICONT.NE.0) THEN 220 : C 221 : C*** LA PARTICULE I EST SORTIE DE L'ELEMENT COURANT 222 : C 222 : CALL TRJIEL(IEL1,IEL2,ICONT,NF,IFACEL,IZCENT,IELTFA) 223 : C write(6,*)' apres TRJIEL 1 ',IEL1,IEL2,' nuapar' ,NUAPAR ,ITER 224 : IF(IEL2 .EQ.0)GO TO 10 225 : GO TO 14 226 : ENDIF 227 : DTC=DTI 228 : DTC1=DTI 229 : TTEMP=TPREC+(LI-0.5D0)*DTI 229 : IF(IVPT.NE.1)CALL TRJTPT(IZVIT,TTEMP,IVPT) 229 : CALL TRJVEL(IZVIT,IZUN,IEL1,IVPT,TTEMP) 230 : NOUN=UN(/2) 230 : CALL TRJULO(XYZL,UN(1,1,IELL),UELEM,XARI,ITY1,NDIM,NOEL1,NOUN, 231 : * IFORMU,SHP,SHY) 232 : DO 17 ID1=1,NDIM 233 : XARI(ID1)=XDEP(ID1)+DTC*UELEM(ID1) 234 : 17 CONTINUE 235 : ITYG=NUMGEO(ITY1) 235 : CALL TRJINT(XARI,XDEP,XINT,UELEM,DTC,DTINT,NDIM,ICONT,ITYG,IO) 236 : C write(6,*)' apres trjint ',io 237 : IF(IO.NE.0)THEN 238 : C 239 : C*** RATTRAPAGE 240 : C 241 : KIO=KIO+1 242 : C IF(MOD(KIO,2).EQ.0)WRITE(6,999)IPART,IEL1,IEL2 243 : IF(MOD(KIO,2).EQ.0)THEN 244 : INTERR(1)=IPART 245 : INTERR(2)=IEL1 245 : CALL ERREUR(-298) 246 : ENDIF 247 : C WRITE(6,*)'TRJRAT 2 ',IPART,IEL1,IEL2, 248 : C * XARI,XDEP,XINT,UELEM,DTT,DTINT,NDIM,ICONT,ITYG,IO 249 : C IF(MOD(KIO,10).EQ.0) STOP 250 : INTERR(1)=IPART 250 : CALL TRJRAT(IEL1,IEL2,XDEP,XARI,IZCOU,MELEME,IZVIT,DTC1, 251 : * IZCENT,IELTFA,XINT,DTINT,ICONT,IZSH,TTEMP) 252 : IF(ICONT.NE.0)GO TO 28 253 : IF(IEL2.EQ.0)THEN 254 : DO 98 ID1=1,NDIM 255 : XINT(ID1)=XDEP(ID1) 256 : 98 CONTINUE 257 : DTINT=0.D0 258 : GO TO 10 259 : ENDIF 260 : IF(IEL1.EQ.IEL2)GO TO 11 261 : DO 18 ID1=1,NDIM 262 : XDEP(ID1)=XARI(ID1) 263 : 18 CONTINUE 264 : TCOUR=TPREC+(LI-1)*DTI+DTC1 265 : NUAPAR=IEL2 266 : DT1=0.D0 267 : GO TO 13 268 : ENDIF 269 : 28 CONTINUE 270 : IF(ICONT.NE.0) THEN 271 : C 272 : C*** LA PARTICULE I EST SORTIE DE L'ELEMENT COURANT 273 : C 273 : CALL TRJIEL(IEL1,IEL2,ICONT,NF,IFACEL,IZCENT,IELTFA) 274 : IF(IEL2 .EQ.0)GO TO 10 275 : GO TO 14 276 : ENDIF 277 : C 278 : 11 CONTINUE 279 : LI=LI+1 280 : DO 24 ID1=1,NDIM 281 : IF(ABS(XARI(ID1)-XDEP(ID1)).GT.1.D-15)GO TO 25 282 : 24 CONTINUE 283 : IF(IIMPI.GT.0)WRITE(6,*)' LA PARTICULE ',IPART , 284 : * ' N AVANCE PLUS ON ARETE LE CALCUL ' 285 : NUAPAR=IEL1 286 : TCOUR=TPREC+NI*DTI 287 : GO TO 21 288 : 25 CONTINUE 289 : DO 19 ID1=1,NDIM 290 : XDEP(ID1)=XARI(ID1) 291 : 19 CONTINUE 292 : IF(LI.LE.NI)GO TO 6 293 : C 294 : C*** LA PARTICULE I EST RESTE DANS L'ELEMENT COURANT 295 : C 296 : NUAPAR=IEL1 297 : TCOUR=TPREC+NI*DTI 298 : IF (ABS((TMAX-TCOUR)/TMAX).LE.1.D-04)GO TO 21 299 : DT1=0.D0 300 : GO TO 13 301 : C 302 : C*** LA PARTICULE I RESTE DANS LE DOMAINE DE CALCUL 303 : C 304 : 14 CONTINUE 304 : CALL MELNEL(IEL2,MELEME,IPT2,NEL2,1) 305 : IELL2=IEL2-NEL2 306 : JTY=IPT2.ITYPEL 307 : CCCCCCCCCCC 308 : C write(6,*) 'iel2 ',iel2,'iel1' ,iel1,'ITER',ITER 309 : 310 : NOEL2=IPT2.NUM(/1) 310 : CALL DOXE(XCOOR,IDIM,NOEL2,IPT2.NUM,IELL2,XYZL) 310 : CALL TRJCN5(JTY,IZSH,IJK) 311 : C IF(IERR.GT.0) RETURN 312 : IF(ISENS.NE.IJK) THEN 313 : IEROR= 1 314 : IEL2= 0 315 : C WRITE (6,*) ' JE VAIS EN 10 ', iel1,ITER,ieror 316 : GO TO 10 317 : ENDIF 318 : CCCCCCCCCCC 318 : CALL TRJFAC(IEL2,NF,JFA,IELTFA) 319 : IOR=1 319 : IF(NDIM.EQ.3)CALL TRJIOR(IELL,IELL2,IPT1,IPT2,ICONT,JFA,IOR) 320 : C WRITE(6,*)' IOR ',IOR,TCOUR 321 : ITYG=NUMGEO(ITY1) 322 : JTYG=NUMGEO(JTY) 322 : CALL TRJTRJ(XDEP,XINT,ITYG,JTYG,JFA,ICONT,IOR) 323 : TCOUR=TPREC+(LI-1)*DTI+DTINT 324 : NUAPAR=IEL2 325 : C DT1=DTINT+(LI-1)*DTI 326 : IF(ABS(DELTAT).GT.1.D-15)THEN 327 : DT1=DTINT+(LI-1)*DTI+DT1 328 : IEL1=NUAPAR 329 : GO TO 1 330 : ENDIF 331 : 13 CONTINUE 332 : DO 22 ID1=1,NDIM 333 : CREF3(ID1,ITER)=XDEP(ID1) 334 : 22 CONTINUE 335 : TPAR(ITER)=TCOUR 336 : NAPAR3(ITER)=NUAPAR 337 : IEL1=NUAPAR 338 : IF ((NPOS-ITER).LE.1) THEN 339 : NPOS=NPOS+50 340 : SEGADJ IZN3 341 : ENDIF 342 : ITER=ITER+1 343 : GO TO 1 344 : C 345 : C*** LA PARTICULE I EST SORTIE DU DOMAINE DE CALCUL 346 : C 347 : 10 CONTINUE 348 : DO 23 ID1=1,NDIM 349 : XDEP(ID1)=XINT(ID1) 350 : 23 CONTINUE 351 : TCOUR=TPREC+(LI-1)*DTI+DTINT 352 : NUAPAR=IEL1 353 : 21 CONTINUE 354 : DO 20 ID1=1,NDIM 355 : CREF3(ID1,ITER)=XDEP(ID1) 356 : 20 CONTINUE 357 : C write(6,*)' IEL1 NUAPAR ITER ', IEL1,NUAPAR,ITER 358 : TPAR(ITER)=TCOUR 359 : NAPAR3(ITER)=NUAPAR 360 : NPOS=ITER 361 : SEGADJ IZN3 362 : C 363 : 999 FORMAT('PARTICULE',I4,': PROBLEMES DANS LE COIN D''UN ELEMENT !!', 364 : * 2I7 ) 365 : C 366 : END 367 : 368 : 369 : 370 : 371 :
© Cast3M 2003 - All rights reserved.
Disclaimer