1 : C TRJRAT SOURCE CHAT 05/01/13 03:50:55 5004 2 : SUBROUTINE TRJRAT(IEL1,IEL2,XDEP,XARI,IZCOU,MELEME,IZVIT,DT, 3 : * IZCENT,IELTFA,XINT,DTINT,ICONT,IZSH,TTEMP) 4 : C 5 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 6 : C 7 : C APPELE PAR TRJAVA 8 : C RATTRAPE LES COORDONNEES DE REFERENCES ET LES APPARTENANCES 9 : C QUAND ON S'EST PERDU 10 : C ( CE CAS SE PRESENTE LORSQUE LA TRAJECTOIRE PASSE PRES DES COINS) 11 : C 12 : C IEL1 NUMERO DE L ELEMENT AVANT DE SE PERDRE 13 : C APRES 14 : C XDEP COORDONNEES DE REFERENCES AU POINT DE PERDITION 15 : C APRES 16 : C IZVIT POINTEUR DU SEGMENT DECRIVANT LES VITESSES 17 : C DT PAS DE TEMPS DU CALCUL 18 : C IZCENT POINTEUR DU MAILLAGE CENTRE DES ELEMENTS 19 : C IELTFA POINTEUR DES CONNECTIONS FACES ELEMENTS 20 : C XINT COORDONNEES DE L'INTERSECTION FACE TRAJECTOIRE 21 : C DTINT TEMPS ELEMENTAIRE AUQUEL LA TRAJECTOIRE TRAVERSE LA FACE 22 : C ICONT NUMERO DE LA FACE PAR LAQUELLE ON SORT DE L'ELEMENT 23 : C TTEMP TEMP AUQUEL EST FAIT LE CALCUL(utilisé en transitoire) 24 : C 25 : C issu de RATTAR de TRIO-EF adapté par F Auriol 26 : C 27 : CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 28 : C 29 : IMPLICIT INTEGER(I-N) 30 : IMPLICIT REAL*8 (A-H,O-Z) 31 : C 32 : C 33 : -INC CCOPTIO 34 : -INC SMCOORD 35 : -INC SMELEME 36 : -INC SMCHPOI 37 : C 38 : POINTEUR IZCENT.MELEME,IELTFA.MELEME 39 : C 40 : SEGMENT IZPART 41 : INTEGER NLEPA(NPART),NUMPA(NPART) 42 : REAL*8 COORPA(NDIM,NPART) 43 : ENDSEGMENT 44 : POINTEUR IZL.IZPART,IZN.IZPART 45 : SEGMENT IZCOU 46 : REAL*8 DTCO(NEL),COU 47 : ENDSEGMENT 48 : C 49 : SEGMENT IZSH 50 : REAL*8 SHP(6,MNO9),SHY(12,MNO9),XYZL(3,MNO9) 51 : ENDSEGMENT 52 : C 53 : SEGMENT IZNOEU 54 : REAL*8 XELE(IDIM,NOEL) 55 : INTEGER NOEGLO(NOEL) 56 : ENDSEGMENT 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 : SEGMENT IZUN 67 : REAL*8 UN(I1,I2,I3) 68 : ENDSEGMENT 69 : C 70 : C 71 : DIMENSION XDEP(3),XARI(3),UE(3),XYL(3),XINT(3) 72 : DIMENSION XYZ(3) 73 : C 74 : C 75 : NDIM=IDIM 76 : NPART=1 77 : ICONT=0 78 : NBNUL=0 79 : SEGINI IZN,IZL 80 : IZN.NLEPA(1)=0 81 : IZN.NUMPA(1)=1 82 : IZL.NLEPA(1)=0 83 : IZL.NUMPA(1)=1 84 : NEL0=0 85 : IPT1=MELEME 86 : NBSOUS=LISOUS(/1) 87 : IF(NBSOUS.NE.0)THEN 88 : DO 10 ISOU=1,NBSOUS 89 : IPT1=LISOUS(ISOU) 90 : SEGACT IPT1 91 : NEL=IPT1.NUM(/2) 92 : IF(NEL+NEL0.GE.IEL1)GO TO 15 93 : NEL0=NEL+NEL0 94 : SEGDES IPT1 95 : 10 CONTINUE 96 : ENDIF 97 : 15 CONTINUE 98 : ITYP=IPT1.ITYPEL 99 : NOEL=IPT1.NUM(/1) 100 : IELL=IEL1-NEL0 101 : NVPT=IPVPT(/1) 102 : IVPT=1 103 : IF(NVPT.NE.1)THEN 104 : IVPT=2 105 : SEGACT IZVIT 105 : CALL TRJTPT(IZVIT,TTEMP,IVPT) 106 : ENDIF 106 : CALL TRJVEL(IZVIT,IZUN,IEL1,IVPT,TTEMP) 107 : 19 CONTINUE 108 : DO 18 I=1,3 109 : XYL(I)=0.D0 110 : UE(I)=0.D0 111 : XYZ(I)=0.D0 112 : 18 CONTINUE 113 : DTI=DT 114 : DO 12 I=1,NDIM 115 : XYZ(I)=XDEP(I) 116 : 12 CONTINUE 117 : C ON CHERCHE SI ON A AFFAIRE A UN NOEUD DU MAILLAGE 118 : ITEST=0 118 : CALL NOEREF(XYZ,ITYP,ITEST) 119 : IF(ITEST.NE.0)THEN 120 : C CAS D'UN NOEUD DU MAILLAGE 121 : C WRITE(6,*)' CAS D UN NOEUD DU MAILLAGE ',XYZ,IEL1 122 : IF(IFORML.EQ.1)THEN 122 : CALL TRJVPO(UN(1,1,IELL),UE,XYZ,ITYP,NDIM,NOEL,SHP) 123 : GO TO 50 124 : ELSEIF(IFORML.EQ.2)THEN 125 : SEGINI IZNOEU 125 : CALL PREXN(IPT1,IELL,IZNOEU) 125 : CALL VHYNOE(UE,IZNOEU,IZUN,NDIM,XYZ,ITYP,IELL) 126 : SEGSUP IZNOEU 127 : GO TO 50 128 : ENDIF 129 : ENDIF 130 : IF(IFORML.EQ.1)THEN 130 : CALL TRJVPO(UN(1,1,IELL),UE,XYZ,ITYP,NDIM,NOEL,SHP) 131 : GO TO 50 132 : ELSEIF(IFORML.EQ.2)THEN 133 : C ON CHERCHE SI XINT EST UN NOEUD DU MAILLAGE 134 : DO 22 I=1,NDIM 135 : XYZ(I)=XINT(I) 136 : 22 CONTINUE 136 : CALL NOEREF(XYZ,ITYP,ITEST) 137 : IF(ITEST.NE.0)THEN 138 : SEGINI IZNOEU 138 : CALL PREXN(IPT1,IELL,IZNOEU) 138 : CALL VHYNOE(UE,IZNOEU,IZUN,NDIM,XYZ,ITYP,IELL) 139 : SEGSUP IZNOEU 140 : DTI=DT-DTINT 141 : GO TO 50 142 : ENDIF 143 : C ON TRAITE LE CAS OU XINT EST SUR UNE ARETE 3D 144 : C ON IRA VERS LA FACE OU LE FLUX EST POSITIF 145 : IF(IDIM.EQ.3)THEN 145 : CALL ARTREF(XDEP,ITYP,IART) 146 : IF(IART.NE.0)THEN 146 : CALL TRJARF(IART,NF1,NF2,ITYP) 147 : IF(UN(1,NF1,IELL).GT.0.D0)THEN 148 : ICONT=NF1 149 : ELSEIF(UN(1,NF2,IELL).GT.0.D0)THEN 150 : ICONT=NF2 151 : ELSE 152 : INTERR(2)= IEL1 153 : INTERR(3)= IART 153 : CALL ERREUR(-316) 154 : IEL2=0 155 : ICONT=0 156 : IF(IIMPI.GT.0)THEN 157 : WRITE(6,*) ' ON RESTE DANS L ELEMENT ',IEL1 158 : WRITE(6,*)' ON EST SUR L ARETE ',IART ,'FLUX',UN(1,NF1,IELL), 159 : * UN(1,NF2,IELL) 160 : ENDIF 161 : SEGSUP IZN,IZL 162 : RETURN 163 : ENDIF 164 : DO 48 ID1=1,NDIM 165 : XINT(ID1)=XDEP(ID1) 166 : 48 CONTINUE 167 : DTINT=0.D0 168 : SEGSUP IZN,IZL 169 : RETURN 170 : ENDIF 171 : ENDIF 172 : C FIN DU TRAITEMENT DES ARETES 173 : C 173 : CALL FACREF(XDEP,ITYP,IFAC) 174 : C WRITE(6,*)' ON EST SUR LA FACE ' ,IFAC ,IEL1,ITYP,XDEP 175 : IF(IFAC.NE.0)THEN 176 : C ON EST SUR LA FACE IFAC 177 : C WRITE(6,*)' ON EST SUR LA FACE ' ,IFAC ,UN(1,IFAC,IELL) 178 : C LA PARTICULE SORTIRA DE L' ELEMENT SI LE FLUX EST POSITIF 179 : IF(UN(1,IFAC,IELL).GT.0.D0)THEN 180 : ICONT=IFAC 181 : DO 45 ID1=1,NDIM 182 : XINT(ID1)=XDEP(ID1) 183 : 45 CONTINUE 184 : DTINT=0.D0 185 : SEGSUP IZN,IZL 186 : RETURN 187 : ELSEIF(UN(1,IFAC,IELL).LT.0.D0)THEN 188 : C ON RESTE DANS L ELEMENT 189 : C revoir ce cas s il se presente(a priori il est difficile a imaginer) 190 : C Il se produit lorsque la vitesse s'annule tres pres de la face sur 191 : C laquelle on se trouve 192 : C ( cas de forte difference de flux et orientation opposé ) 193 : IF(IIMPI.GT.0)THEN 194 : WRITE(6,*) ' ON RESTE DANS L ELEMENT ',IEL1 195 : WRITE(6,*)' ON EST SUR LA FACE ',IFAC ,'FLUX',UN(1,IFAC,IELL) 196 : ENDIF 197 : INTERR(2)= IEL1 198 : INTERR(3)= IFAC 198 : CALL ERREUR(-315) 199 : IEL2=0 200 : ICONT=0 201 : SEGSUP IZN,IZL 202 : RETURN 203 : ELSEIF(UN(1,IFAC,IELL).EQ.0.D0)THEN 204 : C CAS OU LA TRAJECTOIRE SUIT LA FACE 205 : IF(IIMPI.GT.0) 206 : * WRITE(6,*) 'CAS OU LA TRAJECTOIRE SUIT LA FACE ' 207 : ENDIF 208 : IEL2=IEL1 209 : SEGSUP IZN,IZL 210 : RETURN 211 : ELSE 212 : IF(IIMPI.GT.0)THEN 213 : WRITE(6,*)' ON S EST PERDU AU POINT ',XDEP,' DANS L ELEMENT', 214 : * IEL1 215 : ENDIF 216 : INTERR(2)= IEL1 216 : CALL ERREUR(-317) 217 : IEL2=0 218 : ICONT=0 219 : SEGSUP IZN,IZL 220 : RETURN 221 : ENDIF 222 : ENDIF 223 : 50 CONTINUE 223 : CALL SHAPE(XYZ(1),XYZ(2),XYZ(3),ITYP,SHP,IRET) 224 : DO 20 J=1,NOEL 225 : NN=(IPT1.NUM(J,IELL)-1)*(IDIM+1) 226 : DO 25 K=1,NDIM 227 : XYL(K)=XYL(K)+XCOOR(NN+K)*SHP(1,J) 228 : 25 CONTINUE 229 : 20 CONTINUE 230 : DO 2 I=1,NDIM 231 : IZL.COORPA(I,1)=XYL(I)+DTI*UE(I) 232 : C write(6,*)' trjrat xyl ',IZL.COORPA(I,1),XYL(I),DTI,UE(I) 233 : 2 CONTINUE 233 : CALL TRJPEL(IZL,IZN,MELEME,IZVIT,IZCOU,IZCENT,IELTFA,IZSH,TTEMP) 234 : DO 3 I=1,NDIM 235 : XARI(I)=IZN.COORPA(I,1) 236 : 3 CONTINUE 237 : IEL2=IZN.NLEPA(1) 238 : IF(IEL2.NE.0)THEN 239 : C ON VERIFIE QUE IEL2 A AU MOINS UN POINT COMMUN 240 : C AVEC IEL1 SINON ON DIVISE LE PAS DE TEMPS PAR 2 (4 FOIS ) 240 : CALL MELNEL(IEL1,MELEME,IPT2,NEL2,1) 240 : CALL MELNEL(IEL2,MELEME,IPT3,NEL3,1) 241 : NBNN2=IPT2.NUM(/1) 242 : NBNN3=IPT3.NUM(/1) 243 : DO 80 I=1,NBNN2 244 : DO 85 J=1,NBNN3 245 : IF(IPT2.NUM(I,IEL1).EQ.IPT3.NUM(J,IEL2))GO TO 90 246 : 85 CONTINUE 247 : 80 CONTINUE 248 : DT=DT*0.5 249 : GO TO 19 250 : 90 CONTINUE 251 : ELSE 252 : NBNUL= NBNUL+1 253 : IF(NBNUL.LE.4)THEN 254 : DT=DT*0.5 255 : GO TO 19 256 : ELSE 257 : INTERR(2)=IEL1 257 : CALL ERREUR(-318) 258 : ENDIF 259 : ENDIF 260 : SEGSUP IZN,IZL 261 : C 262 : END 263 : 264 : 265 : 266 :
© Cast3M 2003 - All rights reserved.
Disclaimer