C TRJRAT SOURCE CHAT 05/01/13 03:50:55 5004 SUBROUTINE TRJRAT(IEL1,IEL2,XDEP,XARI,IZCOU,MELEME,IZVIT,DT, * IZCENT,IELTFA,XINT,DTINT,ICONT,IZSH,TTEMP) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C C APPELE PAR TRJAVA C RATTRAPE LES COORDONNEES DE REFERENCES ET LES APPARTENANCES C QUAND ON S'EST PERDU C ( CE CAS SE PRESENTE LORSQUE LA TRAJECTOIRE PASSE PRES DES COINS) C C IEL1 NUMERO DE L ELEMENT AVANT DE SE PERDRE C APRES C XDEP COORDONNEES DE REFERENCES AU POINT DE PERDITION C APRES C IZVIT POINTEUR DU SEGMENT DECRIVANT LES VITESSES C DT PAS DE TEMPS DU CALCUL C IZCENT POINTEUR DU MAILLAGE CENTRE DES ELEMENTS C IELTFA POINTEUR DES CONNECTIONS FACES ELEMENTS C XINT COORDONNEES DE L'INTERSECTION FACE TRAJECTOIRE C DTINT TEMPS ELEMENTAIRE AUQUEL LA TRAJECTOIRE TRAVERSE LA FACE C ICONT NUMERO DE LA FACE PAR LAQUELLE ON SORT DE L'ELEMENT C TTEMP TEMP AUQUEL EST FAIT LE CALCUL(utilisé en transitoire) C C issu de RATTAR de TRIO-EF adapté par F Auriol C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C C -INC CCOPTIO -INC SMCOORD -INC SMELEME -INC SMCHPOI C POINTEUR IZCENT.MELEME,IELTFA.MELEME C SEGMENT IZPART INTEGER NLEPA(NPART),NUMPA(NPART) REAL*8 COORPA(NDIM,NPART) ENDSEGMENT POINTEUR IZL.IZPART,IZN.IZPART SEGMENT IZCOU REAL*8 DTCO(NEL),COU ENDSEGMENT C SEGMENT IZSH REAL*8 SHP(6,MNO9),SHY(12,MNO9),XYZL(3,MNO9) ENDSEGMENT C SEGMENT IZNOEU REAL*8 XELE(IDIM,NOEL) INTEGER NOEGLO(NOEL) ENDSEGMENT SEGMENT IZVIT REAL*8 TEMTRA(NVIPT) INTEGER IPUN(NBS),IDUN(NBS),IPVPT(NVIPT),IFORML ENDSEGMENT C IDUN(I) nombre d elements avant le sous maillage I C IPVPT pointeurs de izvpt pour chaque pas de temps SEGMENT IZVPT INTEGER IPUN1(NBS),IPUMAX ENDSEGMENT SEGMENT IZUN REAL*8 UN(I1,I2,I3) ENDSEGMENT C C DIMENSION XDEP(3),XARI(3),UE(3),XYL(3),XINT(3) DIMENSION XYZ(3) C C NDIM=IDIM NPART=1 ICONT=0 NBNUL=0 SEGINI IZN,IZL IZN.NLEPA(1)=0 IZN.NUMPA(1)=1 IZL.NLEPA(1)=0 IZL.NUMPA(1)=1 NEL0=0 IPT1=MELEME NBSOUS=LISOUS(/1) IF(NBSOUS.NE.0)THEN DO 10 ISOU=1,NBSOUS IPT1=LISOUS(ISOU) SEGACT IPT1 NEL=IPT1.NUM(/2) IF(NEL+NEL0.GE.IEL1)GO TO 15 NEL0=NEL+NEL0 SEGDES IPT1 10 CONTINUE ENDIF 15 CONTINUE ITYP=IPT1.ITYPEL NOEL=IPT1.NUM(/1) IELL=IEL1-NEL0 NVPT=IPVPT(/1) IVPT=1 IF(NVPT.NE.1)THEN IVPT=2 SEGACT IZVIT CALL TRJTPT(IZVIT,TTEMP,IVPT) ENDIF CALL TRJVEL(IZVIT,IZUN,IEL1,IVPT,TTEMP) 19 CONTINUE DO 18 I=1,3 XYL(I)=0.D0 UE(I)=0.D0 XYZ(I)=0.D0 18 CONTINUE DTI=DT DO 12 I=1,NDIM XYZ(I)=XDEP(I) 12 CONTINUE C ON CHERCHE SI ON A AFFAIRE A UN NOEUD DU MAILLAGE ITEST=0 CALL NOEREF(XYZ,ITYP,ITEST) IF(ITEST.NE.0)THEN C CAS D'UN NOEUD DU MAILLAGE C WRITE(6,*)' CAS D UN NOEUD DU MAILLAGE ',XYZ,IEL1 IF(IFORML.EQ.1)THEN CALL TRJVPO(UN(1,1,IELL),UE,XYZ,ITYP,NDIM,NOEL,SHP) GO TO 50 ELSEIF(IFORML.EQ.2)THEN SEGINI IZNOEU CALL PREXN(IPT1,IELL,IZNOEU) CALL VHYNOE(UE,IZNOEU,IZUN,NDIM,XYZ,ITYP,IELL) SEGSUP IZNOEU GO TO 50 ENDIF ENDIF IF(IFORML.EQ.1)THEN CALL TRJVPO(UN(1,1,IELL),UE,XYZ,ITYP,NDIM,NOEL,SHP) GO TO 50 ELSEIF(IFORML.EQ.2)THEN C ON CHERCHE SI XINT EST UN NOEUD DU MAILLAGE DO 22 I=1,NDIM XYZ(I)=XINT(I) 22 CONTINUE CALL NOEREF(XYZ,ITYP,ITEST) IF(ITEST.NE.0)THEN SEGINI IZNOEU CALL PREXN(IPT1,IELL,IZNOEU) CALL VHYNOE(UE,IZNOEU,IZUN,NDIM,XYZ,ITYP,IELL) SEGSUP IZNOEU DTI=DT-DTINT GO TO 50 ENDIF C ON TRAITE LE CAS OU XINT EST SUR UNE ARETE 3D C ON IRA VERS LA FACE OU LE FLUX EST POSITIF IF(IDIM.EQ.3)THEN CALL ARTREF(XDEP,ITYP,IART) IF(IART.NE.0)THEN CALL TRJARF(IART,NF1,NF2,ITYP) IF(UN(1,NF1,IELL).GT.0.D0)THEN ICONT=NF1 ELSEIF(UN(1,NF2,IELL).GT.0.D0)THEN ICONT=NF2 ELSE INTERR(2)= IEL1 INTERR(3)= IART CALL ERREUR(-316) IEL2=0 ICONT=0 IF(IIMPI.GT.0)THEN WRITE(6,*) ' ON RESTE DANS L ELEMENT ',IEL1 WRITE(6,*)' ON EST SUR L ARETE ',IART ,'FLUX',UN(1,NF1,IELL), * UN(1,NF2,IELL) ENDIF SEGSUP IZN,IZL RETURN ENDIF DO 48 ID1=1,NDIM XINT(ID1)=XDEP(ID1) 48 CONTINUE DTINT=0.D0 SEGSUP IZN,IZL RETURN ENDIF ENDIF C FIN DU TRAITEMENT DES ARETES C CALL FACREF(XDEP,ITYP,IFAC) C WRITE(6,*)' ON EST SUR LA FACE ' ,IFAC ,IEL1,ITYP,XDEP IF(IFAC.NE.0)THEN C ON EST SUR LA FACE IFAC C WRITE(6,*)' ON EST SUR LA FACE ' ,IFAC ,UN(1,IFAC,IELL) C LA PARTICULE SORTIRA DE L' ELEMENT SI LE FLUX EST POSITIF IF(UN(1,IFAC,IELL).GT.0.D0)THEN ICONT=IFAC DO 45 ID1=1,NDIM XINT(ID1)=XDEP(ID1) 45 CONTINUE DTINT=0.D0 SEGSUP IZN,IZL RETURN ELSEIF(UN(1,IFAC,IELL).LT.0.D0)THEN C ON RESTE DANS L ELEMENT C revoir ce cas s il se presente(a priori il est difficile a imaginer) C Il se produit lorsque la vitesse s'annule tres pres de la face sur C laquelle on se trouve C ( cas de forte difference de flux et orientation opposé ) IF(IIMPI.GT.0)THEN WRITE(6,*) ' ON RESTE DANS L ELEMENT ',IEL1 WRITE(6,*)' ON EST SUR LA FACE ',IFAC ,'FLUX',UN(1,IFAC,IELL) ENDIF INTERR(2)= IEL1 INTERR(3)= IFAC CALL ERREUR(-315) IEL2=0 ICONT=0 SEGSUP IZN,IZL RETURN ELSEIF(UN(1,IFAC,IELL).EQ.0.D0)THEN C CAS OU LA TRAJECTOIRE SUIT LA FACE IF(IIMPI.GT.0) * WRITE(6,*) 'CAS OU LA TRAJECTOIRE SUIT LA FACE ' ENDIF IEL2=IEL1 SEGSUP IZN,IZL RETURN ELSE IF(IIMPI.GT.0)THEN WRITE(6,*)' ON S EST PERDU AU POINT ',XDEP,' DANS L ELEMENT', * IEL1 ENDIF INTERR(2)= IEL1 CALL ERREUR(-317) IEL2=0 ICONT=0 SEGSUP IZN,IZL RETURN ENDIF ENDIF 50 CONTINUE CALL SHAPE(XYZ(1),XYZ(2),XYZ(3),ITYP,SHP,IRET) DO 20 J=1,NOEL NN=(IPT1.NUM(J,IELL)-1)*(IDIM+1) DO 25 K=1,NDIM XYL(K)=XYL(K)+XCOOR(NN+K)*SHP(1,J) 25 CONTINUE 20 CONTINUE DO 2 I=1,NDIM IZL.COORPA(I,1)=XYL(I)+DTI*UE(I) C write(6,*)' trjrat xyl ',IZL.COORPA(I,1),XYL(I),DTI,UE(I) 2 CONTINUE CALL TRJPEL(IZL,IZN,MELEME,IZVIT,IZCOU,IZCENT,IELTFA,IZSH,TTEMP) DO 3 I=1,NDIM XARI(I)=IZN.COORPA(I,1) 3 CONTINUE IEL2=IZN.NLEPA(1) IF(IEL2.NE.0)THEN C ON VERIFIE QUE IEL2 A AU MOINS UN POINT COMMUN C AVEC IEL1 SINON ON DIVISE LE PAS DE TEMPS PAR 2 (4 FOIS ) CALL MELNEL(IEL1,MELEME,IPT2,NEL2,1) CALL MELNEL(IEL2,MELEME,IPT3,NEL3,1) NBNN2=IPT2.NUM(/1) NBNN3=IPT3.NUM(/1) DO 80 I=1,NBNN2 DO 85 J=1,NBNN3 IF(IPT2.NUM(I,IEL1).EQ.IPT3.NUM(J,IEL2))GO TO 90 85 CONTINUE 80 CONTINUE DT=DT*0.5 GO TO 19 90 CONTINUE ELSE NBNUL= NBNUL+1 IF(NBNUL.LE.4)THEN DT=DT*0.5 GO TO 19 ELSE INTERR(2)=IEL1 CALL ERREUR(-318) ENDIF ENDIF SEGSUP IZN,IZL C END