C TRJSAU SOURCE CHAT 05/01/13 03:51:09 5004 SUBROUTINE TRJSAU(XARI,XDEP,UELEM,DTINT,IDIM,ITYEL) C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C FORMULAIRE DES TRAJECTOIRES POUR LES SAUTS CONVECTIFS C DANS UNE MAILLE ( calcul analytique formulation EFMH) C C ENTREES C XDEP POSITION INITIALE C UELEM FLUX AUX FACES C IDIM DIMENSION DE L ESPACE C ITYEL TYPE DE L ELEMENT C DTINT DUREE DU SAUT C SORTIES C XARI POSITION APRES LE SAUT C LTEST INDIQUE SI ON SORT DU TRIANGLE C C TYPES D ELEMENTS CONSIDERES C TRI3 QUA4 CUB8 PRI6 TET C 4 8 14 16 23 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C IMPLICIT INTEGER(I-N) IMPLICIT REAL*8 (A-H,O-Z) C DIMENSION XARI(3),XDEP(3),X(3),Y(3),UELEM(6),XY(4) C C C C*** QUADRANGLE QUA4 QUA9 3 C ***** C * * C FACES 4 * * 2 C ***** C 1 IF(ITYEL.EQ.8)THEN C Q=MAX(ABS(UELEM(1)),ABS(UELEM(2)),ABS(UELEM(3)),ABS(UELEM(4))) C LES FLUX OPPOSES NE SE COMPENSENT PAS C X(1)=UELEM(2)+UELEM(4) X(2)=UELEM(3)+UELEM(1) Y(1)=UELEM(2)-UELEM(4) Y(2)=UELEM(3)-UELEM(1) C IF(ABS(X(1))/Q.GT.1D-10)THEN c write(6,*) 'coucou1' DO 8 I2=1,IDIM XARI(I2)=(XDEP(I2)*EXP(25.D-2*X(I2)*DTINT))- * ((Y(I2)/X(I2))*(1-EXP(25.D-2*X(I2)*DTINT))) 8 CONTINUE ELSE C C LES FLUX OPPOSES SE COMPENSENT C DO 2 I3=1,IDIM XARI(I3)=XDEP(I3)+(25.D-2*Y(I3)*DTINT) 2 CONTINUE ENDIF C*** CUBE CUB8 C C LES FLUX OPPOSES NE SE COMPENSENT PAS C ELSEIF(ITYEL.EQ.14)THEN Q=MAX(ABS(UELEM(1)),ABS(UELEM(2)),ABS(UELEM(3)), * ABS(UELEM(4)),ABS(UELEM(5)),ABS(UELEM(6))) X(1)=125.D-3*(UELEM(4)+UELEM(6)) X(2)=125.D-3*(UELEM(3)+UELEM(5)) X(3)=125.D-3*(UELEM(2)+UELEM(1)) C Y(1)=125.D-3*(UELEM(6)-UELEM(4)) Y(2)=125.D-3*(UELEM(3)-UELEM(5)) Y(3)=125.D-3*(UELEM(1)-UELEM(2)) C IF(ABS(X(1))/Q.GT.1.D-11)THEN IF(ABS(X(2))/Q.GT.1.D-11)THEN IF(ABS(X(3))/Q.GT.1.D-11)THEN DO 3 I4=1,IDIM XARI(I4)=(XDEP(I4)*EXP(125.D-3*X(I4)*DTINT))+ * ((Y(I4)/X(I4))*(1-EXP(125.D-3*X(I4)*DTINT))) 3 CONTINUE C C SI L'EQUATION SUIVANT Z DIFFERE C ELSE XARI(1)=(XDEP(1)*EXP(125.D-3*X(1)*DTINT))+ * ((Y(1)/X(1))*(1-EXP(125.D-3*X(1)*DTINT))) XARI(2)=(XDEP(2)*EXP(125.D-3*X(2)*DTINT))+ * ((Y(2)/X(2))*(1-EXP(125.D-3*X(2)*DTINT))) XARI(3)=XDEP(3)-(Y(3)*DTINT) C C SI L'EQUATION SUIVANT Y DIFFERE C ENDIF ELSE XARI(1)=(XDEP(1)*EXP(125.D-3*X(1)*DTINT))+ * ((Y(1)/X(1))*(1-EXP(125.D-3*X(1)*DTINT))) XARI(3)=(XDEP(3)*EXP(125.D-3*X(3)*DTINT))+ * ((Y(3)/X(3))*(1-EXP(125.D-3*X(3)*DTINT))) XARI(2)=XDEP(2)-(Y(2)*DTINT) C C SI L'EQUATION SUIVANT X DIFFERE C ENDIF ELSE IF(ABS(X(2))/Q.GT.1.D-11)THEN XARI(3)=(XDEP(3)*EXP(125.D-3*X(3)*DTINT))+ * ((Y(3)/X(3))*(1-EXP(125.D-3*X(3)*DTINT))) XARI(2)=(XDEP(2)*EXP(125.D-3*X(2)*DTINT))+ * ((Y(2)/X(2))*(1-EXP(125.D-3*X(2)*DTINT))) XARI(1)=XDEP(1)-(Y(1)*DTINT) C C TOUS LES FLUX SE COMPENSENT ELSE DO 4 I5=1,IDIM XARI(I5)=XDEP(I5)-(Y(I5)*DTINT) 4 CONTINUE ENDIF ENDIF CC CC CC*** PRISME PRI6 CC ELSEIF(ITYEL.EQ.16)THEN C Q=MAX(ABS(UELEM(1)),ABS(UELEM(2)),ABS(UELEM(3)), * ABS(UELEM(4)),ABS(UELEM(5))) VAR=(UELEM(3)+UELEM(4)+UELEM(5))/2.D0 C SI LES FLUX NE SE COMPENSENT PAS C IDIM=3 VAR=(UELEM(3)+UELEM(4)+UELEM(5)) X(1)=UELEM(5) X(2)=UELEM(3) X(3)=(UELEM(2)-UELEM(1)) C IF(ABS(VAR)/Q.GT.1.D-10)THEN XARI(1)=(XDEP(1)*EXP(5.D-1*VAR*DTINT))+((X(1)/VAR)* * (1.D0-EXP(5.D-1*VAR*DTINT))) XARI(2)=(XDEP(2)*EXP(5.D-1*VAR*DTINT))+((X(2)/VAR)* * (1.D0-EXP(5.D-1*VAR*DTINT))) XARI(3)=(XDEP(3)*EXP(-1.D0*VAR*DTINT))+((X(3)/VAR)* * (1.D0-EXP(-1.D0*VAR*DTINT))) C ELSE C C SI LES FLUX SE COMPENSENT C XARI(1)=XDEP(1)-(5.D-1*X(1)*DTINT) XARI(2)=XDEP(2)-(5.D-1*X(2)*DTINT) XARI(3)=XDEP(3)+(X(3)*DTINT) C ENDIF ENDIF END