Sources Esope | Cast3M

trjsau
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