Télécharger la plaquette Cast3M


fr en

Télécharger trjma2.eso

Retour à la liste

Numérotation des lignes : 

trjma2

  1. C TRJMA2    SOURCE    PV        08/09/11    21:17:12     6150
  2.       SUBROUTINE TRJMA2(XARI,XDEP,UELEM,DTINT,
  3.      *                  IDIM,ICONT,ITYEL,IART,INOEU)
  4. C
  5. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  6. C     CALCUL DES TRAJECTOIRES ET
  7. C     INTERSECTION TRAJECTOIRE-ELEMENT ( calcul analytique pour une
  8. C     formulation EFMH)
  9. C
  10. C   ENTREES
  11. C     XDEP POSITION INITIALE
  12. C     UELEM FLUX AUX FACES
  13. C     IDIM DIMENSION DE L ESPACE
  14. C     ITYEL TYPE DE L ELEMENT
  15. C
  16. C   SORTIES
  17. C     XARI POSITION A L INTERSECTION
  18. C     DTINT TEMPS ECOULE JUSQU A LA SORTIE
  19. C     ICONT NUMERO DE LA PORTE DE SORTIE
  20. C     IART   DIFFERENT DE 0 ON EST SORTI DE L ELEMENT 3D PAR UNE ARETE
  21. C     INOEU   DIFFERENT DE 0 ON EST SORTI DE L ELEMENT NOEUD
  22. C
  23. C   TYPES D ELEMENTS CONSIDERES
  24. C         TRI3  QUA4  CUB8 PRI6 TET
  25. C         4     8     14   16    23
  26. C
  27. C         Auteur Patrick Meyniel
  28. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  29. C
  30.       IMPLICIT INTEGER(I-N)
  31.       IMPLICIT REAL*8 (A-H,O-Z)
  32.       LOGICAL TEST
  33. C
  34.       DIMENSION IJ(2),ZI(2),IFAT3(3),YI(3)
  35.       DIMENSION IFAQ4(4),IJKL(3,3),IJT4(2,3),IFAT4(4),IFAC6(4),IFAC8(6)
  36.       DIMENSION XARI(3),XDEP(3),UELEM(6),X(3),Y(3)
  37. C
  38.       TEST(DTINT,DTI)=(DTINT*DTI*(DTINT-DTI).GE.0.D0)
  39.       DATA IJKL/1,2,3,2,1,3,3,2,1/
  40.       DATA IJ /2,1/
  41.       DATA IJT4 /2,3, 1,3, 1,2/
  42.       DATA IFAT3 /3,1,2/
  43.       DATA IFAQ4 /1,3,4,2/
  44.       DATA IFAC6 /6,4,1,2/
  45.       DATA IFAC8 /6,4,3,5,1,2/
  46.       DATA IFAT4 /4,1,2,3/
  47.       DATA YI/5.D-1,5.D-1,-1.D0/
  48.       DATA ZI/-1.D0,1.D0/
  49.  
  50. C
  51.       ICONT=0
  52.       IART=0
  53.       INOEU=0
  54.       ITEST=0
  55.       EPS=1.D-14
  56.       DTINT=1.D30
  57.       DO 50 I=1,IDIM
  58.       XARI(I)=XDEP(I)
  59.    50 CONTINUE
  60.  
  61. C
  62. C***  TRIANGLES TRI3                       *
  63. C                                          * *
  64. C                               FACES    3 *  * 2
  65. C                                          *   *
  66. C                                          ******
  67. C                                            1
  68.        IF(ITYEL.EQ.4)THEN
  69.          IDIM=2
  70. C
  71. C             A PARTIR DES FLUX ON RECHERCHE LE POINT DE SORTIE DE
  72. C             LA TRAJECTOIRE DANS L'ELEMENT
  73. C
  74. C             ON TESTE L'INTERSECTION AVEC LA FACE 2 ET TOUTES LES
  75. C                        FACES AYANT UN FLUX POSITIF
  76. C  FACE 2
  77. C
  78.                 IF(UELEM(2).GT.0.D0)THEN
  79.                    DTINT=(XDEP(2)+XDEP(1)-1)/(UELEM(1)+UELEM(3))
  80.                    XARI(1)=XDEP(1)-(UELEM(3)*DTINT)
  81.                    XARI(2)=XDEP(2)-(UELEM(1)*DTINT)
  82.                    ICONT=2
  83.                 ENDIF
  84. C
  85. C  FACE 1 ET 3
  86. C
  87.                 DO 1 I=1,IDIM
  88. C
  89.                    IF(UELEM(IFAT3(I)).GT.0.D0)THEN
  90.                       DTI=XDEP(I)/UELEM(IFAT3(I))
  91. C
  92. C      ON NE GARDE QUE LE PLUS PETIT TEMPS POSITIF
  93. C
  94.                       IF(TEST(DTINT,DTI))THEN
  95.                          DTINT=DTI
  96. C       ON GARDE LES COORD DU POINT ET LA FACE DE SORTIE
  97. C       CORRESPONDANTS AU TEMPS LE PLUS PETIT
  98.                          XARI(IJ(I))=XDEP(IJ(I))-(UELEM(IFAT3(IJ(I)))
  99.      *                        *DTINT)
  100.                          XARI(I)=0.D0
  101.                          ICONT=IFAT3(I)
  102.                       ENDIF
  103.                    ENDIF
  104.  1              CONTINUE
  105. C
  106. C        ON CHERCHE SI ON A AFFAIRE A UN NOEUD DU MAILLAGE
  107. C
  108.  
  109.              DO 30 I=1,IDIM
  110.              IF(ABS(XARI(I)).LT.EPS)ITEST=ITEST+1
  111.              IF(ABS(1.D0-XARI(I)).LT.EPS)ITEST=ITEST+1
  112.    30        CONTINUE
  113. C
  114. C
  115. C
  116. C***  QUADRANGLE       QUA4 QUA9             3
  117. C                                           *****
  118. C                                           *   *
  119. C                                 FACES  4  *   * 2
  120. C                                           *****
  121. C                                             1
  122.        ELSEIF(ITYEL.EQ.8)THEN
  123.        IDIM=2
  124.         Q=MAX(ABS(UELEM(1)),ABS(UELEM(2)),ABS(UELEM(3)),ABS(UELEM(4)))
  125. C        SI LES FLUX OPPOSES NE SE COMPENSENT PAS
  126. C
  127.         X(1)=UELEM(3)+UELEM(1)
  128.         X(2)=UELEM(2)+UELEM(4)
  129.         Y(1)=UELEM(3)-UELEM(1)
  130.         Y(2)=UELEM(2)-UELEM(4)
  131. C
  132.         IF(ABS(X(1))/Q.GT.1D-10)THEN
  133. C             ON TESTE TOUTES LES FACES DE SORTIE QUI ONT UN FLUX POSITIF
  134.  
  135.           DO 2 I=1,IDIM
  136.              DO 3 I1=1,IDIM
  137.                 IF(UELEM(IFAQ4(2*I+I1-2)).GT.0.D0)THEN
  138.                    VAR1=ZI(I1)+(Y(I)/X(I))
  139.                    VAR2=XDEP(IJ(I))+(Y(I)/X(I))
  140. C       ON ELIMINE LES CAS OÙ LE LOG EST NEG OU NON DEFINI
  141.                    IF((VAR1/VAR2).LE.0.D0)GOTO 7
  142.                    DTI=(4.D0/X(I))*LOG(VAR1/VAR2)
  143.                    IF(ABS(LOG(VAR1/VAR2)).LT.1.D-15)DTI=0.D0
  144. C             ON GARDE LE PLUS PETIT PAS DE TEMPS
  145.                    IF(TEST(DTINT,DTI))THEN
  146.                       DTINT=DTI
  147.                       I2=I
  148.                       I3=I1
  149.                    ENDIF
  150.                 ENDIF
  151.  7              CONTINUE
  152.  3           CONTINUE
  153.  2        CONTINUE
  154. C
  155. C        ON EXHIBE LES COORDONNEES CORRESPONDANT AU PLUS PETIT TEMPS
  156.           XARI(I2)=(XDEP(I2)*EXP(25.D-2*X(IJ(I2))*DTINT))-
  157.      *         ((Y(IJ(I2))/X(IJ(I2)))*(1.-EXP(25.D-2*X(IJ(I2))*DTINT)))
  158. C
  159. C         SINON (LES FLUX OPPOSES SE COMPENSENT)
  160. C
  161.        ELSE
  162. C
  163. C
  164. C           ON TESTE TOUTES LES FACES AYANT UN FLUX POSITIF
  165. C
  166.           DO 4 I=1,IDIM
  167.              DO 5 I1=1,IDIM
  168.                 IF(UELEM(IFAQ4(2*I+I1-2)).GT.0.D0)THEN
  169. C
  170.                    DTI=(ZI(I1)-XDEP(IJ(I)))/(25.D-2*Y(I))
  171. C           ON GARDE LE PLUS PETIT TEMPS POSITIF
  172.                    IF(TEST(DTINT,DTI))THEN
  173.                       DTINT=DTI
  174.                       I2=I
  175.                       I3=I1
  176.                    ENDIF
  177.                 ENDIF
  178.  5           CONTINUE
  179.  4        CONTINUE
  180. C
  181. C     ON GARDE LES COORDONNEES CORRESPONDANT AU PLUS PETIT TEMPS
  182.           XARI(I2)=XDEP(I2)+25.D-2*Y(IJ(I2))*DTINT
  183.        ENDIF
  184.        XARI(IJ(I2))=ZI(I3)
  185. C     ON EXHIBE LE NUMERO DE LA FACE DE SORTIE
  186.        ICONT=IFAQ4((2*I2)+I3-2)
  187. C
  188. C        ON CHERCHE SI ON A AFFAIRE    A UN NOEUD DU MAILLAGE
  189.  
  190.        DO 35 I=1,IDIM
  191.           IF(ABS(1.D0-XARI(I)).LT.EPS)ITEST=ITEST+1
  192.           IF(ABS(1.D0+XARI(I)).LT.EPS)ITEST=ITEST+1
  193.  35    CONTINUE
  194. C
  195. C
  196. C***   CUBE      CUB8
  197.       ELSEIF(ITYEL.EQ.14)THEN
  198.          IDIM=3
  199.          Q=MAX(ABS(UELEM(1)),ABS(UELEM(2)),ABS(UELEM(3)),
  200.      *        ABS(UELEM(4)),ABS(UELEM(5)),ABS(UELEM(6)))
  201.         X(1)=125.D-3*(UELEM(4)+UELEM(6))
  202.         X(2)=125.D-3*(UELEM(3)+UELEM(5))
  203.         X(3)=125.D-3*(UELEM(2)+UELEM(1))
  204. C
  205.         Y(1)=125.D-3*(UELEM(6)-UELEM(4))
  206.         Y(2)=125.D-3*(UELEM(3)-UELEM(5))
  207.         Y(3)=125.D-3*(UELEM(1)-UELEM(2))
  208. C
  209. C    SI LES FLUX NE SECOMPENSENT PAS DEUX A DEUX
  210. C
  211.         IF(ABS(X(1))/Q.GT.1.D-11)THEN
  212.            IF(ABS(X(2))/Q.GT.1.D-11)THEN
  213.                IF(ABS(X(3))/Q.GT.1.D-11)THEN
  214. C
  215. C             ON TESTE TOUTES LES FACES DE SORTIE
  216. C
  217.                   DO 17 I=1,IDIM
  218.                      DO 18 I1=1,2
  219. C             ON TESTE TOUTES LES FACES QUI ONT UN FLUX POSITIF
  220.                         IF(UELEM(IFAC8(2*I+I1-2)).GT.0.D0)THEN
  221.                            VAR1=ZI(I1)-(Y(I)/X(I))
  222.                            VAR2=XDEP(I)-(Y(I)/X(I))
  223. C             ON ELIMINE LES CAS OÙ LE LOG EST NEG OU NON DEFINI
  224.                            IF((VAR1/VAR2).LE.0.D0)GOTO 6
  225.                            DTI=(1.D0/X(I))*LOG(VAR1/VAR2)
  226.                            IF(ABS(LOG(VAR1/VAR2)).LT.1.D-15)DTI=0.D0
  227. C             ON GARDE LE PLUS PETIT TEMPS POSITIF
  228.                            IF(TEST(DTINT,DTI))THEN
  229.                               DTINT=DTI
  230.                               I2=I
  231.                               I3=I1
  232.                               I4=IJT4(1,I)
  233.                               I5=IJT4(2,I)
  234.                            ENDIF
  235.                         ENDIF
  236.  6                      CONTINUE
  237.  18                 CONTINUE
  238.  17               CONTINUE
  239. C    ON GARDE LES COORDONNEES CORRESPONDANT AU TEMPS LE PLUS PETIT
  240.                   XARI(I2)=ZI(I3)
  241.                   XARI(I4)=(XDEP(I4)*EXP(X(I4)*DTINT))+
  242.      *                 ((Y(I4)/X(I4))*(1-EXP(X(I4)*DTINT)))
  243.                   XARI(I5)=(XDEP(I5)*EXP(X(I5)*DTINT))+
  244.      *                 ((Y(I5)/X(I5))*(1-EXP(X(I5)*DTINT)))
  245.                   ICONT=IFAC8((2*I2)+I3-2)
  246.  
  247.                ELSE
  248. C      ON TESTE TOUTES LES FACES MAIS L'EQUATION SUIVANT Z DIFFERE
  249. C                                (Q1+Q2=O)
  250.                   DO 21 I=1,2
  251.                      DO 19 I1=1,2
  252.                         IF(UELEM(IFAC8(2*I+I1-2)).GT.0.D0)THEN
  253.                            VAR1=ZI(I1)-(Y(I)/X(I))
  254.                            VAR2=XDEP(I)-(Y(I)/X(I))
  255. C              ON ELIMINE LES CAS OÙ LE LOG EST NEG OU NON DEFINI
  256.                       IF((VAR1/VAR2).LE.0.D0)GOTO 9
  257.                            DTI=(1.D0/X(I))*LOG(VAR1/VAR2)
  258.                            IF(ABS(LOG(VAR1/VAR2)).LT.1.D-15)DTI=0.D0
  259. c              ON GARDE LE PLUS PETIT TEMPS POSITIF
  260.                            IF(TEST(DTINT,DTI))THEN
  261.                               DTINT=DTI
  262.                               I2=I
  263.                               I3=I1
  264.                               I4=IJT4(1,I)
  265.                            ENDIF
  266.                         ENDIF
  267.  9                      CONTINUE
  268.  19                  CONTINUE
  269.  21               CONTINUE
  270. C           ON GARDE LES COORDONNEES ET LA FACE CORRESPONDANT AU TEMPS
  271.                   XARI(I2)=ZI(I3)
  272.                   XARI(I4)=(XDEP(I4)*EXP(X(I4)*DTINT))+
  273.      *                 ((Y(I4)/X(I4))*(1-EXP(X(I4)*DTINT)))
  274.                   XARI(3)=XDEP(3)-(Y(3)*DTINT)
  275.                   ICONT=IFAC8((2*I2)+I3-2)
  276. C          ON TESTE MAINTENANT LES FACES  Z=1 ET Z=-1
  277.                   IF(UELEM(2).GT.0.D0)THEN
  278.                      DTI=(XDEP(3)-1.D0)/Y(3)
  279.                      INT=2
  280.                   ELSE
  281.                      DTI=(XDEP(3)+1.D0)/Y(3)
  282.                      INT=1
  283.                   ENDIF
  284. C                    ON TESTE LE TEMPS
  285.                   IF(TEST(DTINT,DTI))THEN
  286.                      DTINT=DTI
  287.                      XARI(1)=(XDEP(1)*EXP(X(1)*DTINT))+
  288.      *                    ((Y(1)/X(1))*(1-EXP(X(1)*DTINT)))
  289.                      XARI(2)=(XDEP(2)*EXP(X(2)*DTINT))+
  290.      *                    ((Y(2)/X(2))*(1-EXP(X(2)*DTINT)))
  291.                      XARI(3)=ZI(INT)
  292.                      ICONT=INT
  293.                   ENDIF
  294.                ENDIF
  295.            ELSE
  296. C
  297. C   ON A OBLIGATOIREMENT Q1+Q2 NON NUL DONC ON PEUT RAISONNER
  298. C                          COMME PRECEDEMMENT
  299. C     ON TESTE TOUTES LES FACES MAIS L'EQUATION SUIVANT Y DIFFERE
  300. C                         (Q3+Q5=0)
  301.               DO 23 I=1,2
  302.                  DO 22 I1=1,2
  303. C
  304.                     IF(UELEM(IFAC6(2*I+I1-2)).GT.0.D0)THEN
  305.                        VAR1=ZI(I1)-(Y(I)/X(I))
  306.                        VAR2=XDEP(I)-(Y(I)/X(I))
  307. C
  308.                       IF((VAR1/VAR2).LE.0.D0)GOTO 11
  309.                       DTI=(1.D0/X(I))*LOG(VAR1/VAR2)
  310. C
  311.                       IF(ABS(LOG(VAR1/VAR2)).LT.1.D-15)DTI=0.D0
  312. C
  313.                       IF(TEST(DTINT,DTI))THEN
  314.                          DTINT=DTI
  315.                          I2=IFAQ4(I)
  316.                          I3=I1
  317.                          I4=IFAT3(I)
  318.                       ENDIF
  319.                    ENDIF
  320.  11                CONTINUE
  321.  22             CONTINUE
  322.  23          CONTINUE
  323.               XARI(I2)=ZI(I3)
  324.               XARI(I4)=(XDEP(I4)*EXP(X(I4)*DTINT))+
  325.      *             ((Y(I4)/X(I4))*(1-EXP(X(I4)*DTINT)))
  326.               XARI(2)=XDEP(2)-(Y(2)*DTINT)
  327.               ICONT=IFAC6((2*I2)+I3-2)
  328. C
  329. C           ON TESTE MAINTENANT LES FACES Y=1 ET Y=-1
  330. C
  331.               IF(UELEM(5).GT.0.D0)THEN
  332.                  DTI=(XDEP(2)-1.D0)/Y(2)
  333.                  INT=5
  334.               ELSE
  335.                  DTI=(XDEP(2)+1.D0)/Y(2)
  336.                  INT=3
  337.               ENDIF
  338.               IF(TEST(DTINT,DTI))THEN
  339.                  DTINT=DTI
  340.                  XARI(1)=(XDEP(1)*EXP(X(1)*DTINT))+
  341.      *                ((Y(1)/X(1))*(1-EXP(X(1)*DTINT)))
  342.                  XARI(3)=(XDEP(3)*EXP(X(3)*DTINT))+
  343.      *                ((Y(3)/X(3))*(1-EXP(X(3)*DTINT)))
  344.                  XARI(2)=INT-4
  345.                  ICONT=INT
  346.               ENDIF
  347.            ENDIF
  348. C
  349.         ELSE
  350. C        ON TESTE TOUTE LES FACES MAIS L'EQUATION DE X DIFFERE
  351. C                   (Q4+Q6 =0)
  352.            IF(ABS(X(2))/Q.GT.1.D-11)THEN
  353.               DO 24 I=2,3
  354.                  DO 25 I1=1,2
  355. C
  356.                     IF(UELEM(IFAC8(2*I+I1-2)).GT.0.D0)THEN
  357. C
  358.                        VAR1=ZI(I1)-(Y(I)/X(I))
  359.                        VAR2=XDEP(I)-(Y(I)/X(I))
  360. C
  361.                       IF((VAR1/VAR2).LE.0.D0)GOTO 13
  362. C
  363.                        DTI=(1.D0/X(I))*LOG(VAR1/VAR2)
  364.                        IF(ABS(LOG(VAR1/VAR2)).LT.1.D-15)DTI=0.D0
  365. C
  366.                        IF(TEST(DTINT,DTI))THEN
  367.                           DTINT=DTI
  368.                           I2=I
  369.                           I3=I1
  370.                           I4=IJT4(2,I)
  371.                        ENDIF
  372.                     ENDIF
  373.  13                 CONTINUE
  374.  25              CONTINUE
  375.  24           CONTINUE
  376.               XARI(I2)=ZI(I3)
  377.               XARI(I4)=(XDEP(I4)*EXP(X(I4)*DTINT))+
  378.      *             ((Y(I4)/X(I4))*(1-EXP(X(I4)*DTINT)))
  379.               XARI(1)=XDEP(1)-(Y(1)*DTINT)
  380.               ICONT=IFAC8((2*I2)+I3-2)
  381. C          ON TESTE MAINTENANT LES FACES X=1 ET X=-1
  382.               IF(UELEM(6).GT.0.D0)THEN
  383.                  DTI=(XDEP(1)+1.D0)/Y(1)
  384.                  INT=6
  385.               ELSE
  386.                  DTI=(XDEP(1)-1.D0)/Y(1)
  387.                  INT=4
  388.               ENDIF
  389.               IF(TEST(DTINT,DTI))THEN
  390.                  DTINT=DTI
  391.                  XARI(3)=(XDEP(3)*EXP(X(3)*DTINT))+
  392.      *                ((Y(3)/X(3))*(1-EXP(X(3)*DTINT)))
  393.                  XARI(2)=(XDEP(2)*EXP(X(2)*DTINT))+
  394.      *                ((Y(2)/X(2))*(1-EXP(X(2)*DTINT)))
  395.                  XARI(1)=ZI(INT-5)
  396.                  ICONT=INT
  397.               ENDIF
  398.            ELSE
  399. C
  400. C            Q3+Q5  NEGLIGEABLE (DONC Q1+Q2  NEG)
  401. C                     ON CALCULE POUR TOUTES LES FACES
  402.               DO 31 I=1,3
  403.                  DO 32 I1=1,2
  404.                     IF(UELEM(IFAC8(2*I+I1-2)).GT.0.D0)THEN
  405. C
  406.                        DTI=(XDEP(I)-ZI(I1))/Y(I)
  407. C
  408.                        IF(TEST(DTINT,DTI))THEN
  409.                           DTINT=DTI
  410.                           I2=I
  411.                           I3=I1
  412.                           I4=IJT4(1,I)
  413.                           I5=IJT4(2,I)
  414.                        ENDIF
  415.                     ENDIF
  416.  32              CONTINUE
  417.  31           CONTINUE
  418.               XARI(I2)=ZI(I3)
  419.               XARI(I4)=XDEP(I4)-(Y(I4)*DTINT)
  420.               XARI(I5)=XDEP(I5)-(Y(I5)*DTINT)
  421.               ICONT=IFAC8((2*I2)+I3-2)
  422.            ENDIF
  423.         ENDIF
  424. C
  425. C        ON CHERCHE SI ON A AFFAIRE A UN NOEUD DU MAILLAGE
  426. C
  427.         DO 45 I=1,IDIM
  428.            IF(ABS(1.D0+XARI(I)).LT.EPS)ITEST=ITEST+1
  429.            IF(ABS(1.D0-XARI(I)).LT.EPS)ITEST=ITEST+1
  430.  45     CONTINUE
  431. CC
  432. CC
  433. CC***  PRISME    PRI6
  434. CC
  435.       ELSEIF(ITYEL.EQ.16)THEN
  436.          IDIM=3
  437.          Q=MAX(ABS(UELEM(1)),ABS(UELEM(2)),ABS(UELEM(3)),
  438.      *        ABS(UELEM(4)),ABS(UELEM(5)))
  439.          VAR=(UELEM(3)+UELEM(4)+UELEM(5))/2.D0
  440.          X(1)=UELEM(5)/2.D0
  441.          X(2)=UELEM(3)/2.D0
  442.          X(3)=(UELEM(2)-UELEM(1))/2.D0
  443. C
  444.          IF(ABS(VAR)/Q.GT.1.D-10)THEN
  445. C     ON TESTE TOUTES LES FACES QUI ONT UN FLUX POSITIF
  446. C          ET ON GARDE LE PLUS PETIT TEMPS POSITIF
  447. C          ON TESTE LES FACES TRIANGULAIRES: FACE 4
  448.             IF(UELEM(4).GT.0.D0)THEN
  449.                A1=((1.D0-((X(1)+X(2))/VAR)))
  450.                A2=(XDEP(1)+XDEP(2)-((X(1)+X(2))/VAR))
  451.                IF (A1/A2.LE.0.D0)GOTO 15
  452.                DTINT=(1.D0/VAR)*LOG(A1/A2)
  453.                   XARI(2)=(XDEP(2)*EXP(VAR*DTINT))+((X(1)/VAR)*
  454.      *                 (1.D0-EXP(VAR*DTINT)))
  455.                   XARI(3)=(XDEP(3)*EXP(-1.D0*VAR*DTINT))+((X(3)/VAR)*
  456.      *                 (1.D0-EXP(-1.D0*VAR*DTINT)))
  457.                XARI(1)=1.D0-XARI(2)
  458.                ICONT=4
  459.  15            CONTINUE
  460.             ENDIF
  461. C
  462.             DO 8 I=1,IDIM-1
  463. C                        FACE 3 ET FACE 5
  464.                IF(X(I).GT.0.D0)THEN
  465.                  A3=(-1.D0*X(I)/VAR)
  466.                   A4=(XDEP(I)-(X(I)/VAR))
  467.                   IF(A3/A4.LE.0.D0)GOTO 16
  468.                   DTI=(1.D0/VAR)*LOG(A3/A4)
  469. C
  470. C              ON NE GARDE QUE LE PLUS PETIT TEMPS POSITIF
  471. C
  472.                   IF(TEST(DTINT,DTI))THEN
  473.                      DTINT=DTI
  474.                      XARI(IJ(I))=(XDEP(IJ(I))*EXP(VAR*DTINT))+
  475.      *                    ((X(IJ(I))/VAR)*(1.D0-EXP(VAR*DTINT)))
  476.                      XARI(3)=(XDEP(3)*EXP(VAR*DTINT))+
  477.      *                    ((X(3)/VAR)*(1.D0-EXP(-1.D0*VAR*DTINT)))
  478.                      XARI(I)=0.D0
  479.                      ICONT=IFAT3(I)+2
  480.                   ENDIF
  481.  16               CONTINUE
  482.                ENDIF
  483. C              ON TESTE LES FACES OPPOSEES
  484.                IF(UELEM(I).GT.0.D0)THEN
  485.                   A5=(ZI(I)-(X(3)/VAR))
  486.                   A6=(XDEP(3)-(X(3)/VAR))
  487.                    IF(A5/A6.LE.0.D0)GOTO 26
  488.                    DTI=(-5.D-1/VAR)*LOG(A5/A6)
  489. C
  490. C              ON NE GARDE QUE LE PLUS PETIT TEMPS POSITIF
  491. C
  492.                   IF(TEST(DTINT,DTI))THEN
  493.                      DTINT=DTI
  494.                      XARI(I)=(XDEP(I)*EXP(VAR*DTINT))+
  495.      *                    ((X(I)/VAR)*(1.D0-EXP(VAR*DTINT)))
  496.                      XARI(IJ(I))=(XDEP(IJ(I))*EXP(VAR*DTINT))+
  497.      *                  ((X(IJ(I))/VAR)*(1.D0-EXP(VAR*DTINT)))
  498.                      XARI(3)=ZI(I)
  499.                      ICONT=I
  500.                   ENDIF
  501.  26                CONTINUE
  502.                 ENDIF
  503.  8          CONTINUE
  504. C
  505.          ELSE
  506. C
  507. C      Q3 NON NEGLIGEABLE ON VA TESTER LES FACES TRIANGULAIRES
  508. C                 FACE 4
  509.             IF(UELEM(4).GT.0.D0)THEN
  510.                DTINT=(2.D0/UELEM(4))*(1.D0-XDEP(1)-XDEP(2))
  511.                XARI(2)=XDEP(2)-(X(2)*DTINT)
  512.                XARI(1)=1.D0-XARI(2)
  513.                ICONT=4
  514.             ENDIF
  515. C
  516.             DO 12 I=1,IDIM-1
  517.                IF(X(I).GT.0.D0)THEN
  518.                   DTI=XDEP(I)/X(I)
  519. C      ON NE GARDE QUE LE PLUS PETIT TEMPS POSITIF
  520.                   IF(TEST(DTINT,DTI))THEN
  521.                      DTINT=DTI
  522.                      XARI(IJ(I))=XDEP(IJ(I))-(X(IJ(I))*DTINT)
  523.                      XARI(I)=0.D0
  524.                      ICONT=IFAT3(I)+2
  525.                   ENDIF
  526.                ENDIF
  527.  12         CONTINUE
  528.             XARI(3)=XDEP(3)+(2.D0*X(3)*DTINT)
  529. C
  530. C   SI Q1 NON NEGLIGEABLE  ON TESTE  LES FACES OPPOSEES A FLUX POSITIF
  531.             IF(UELEM(1).GT.0.D0)THEN
  532.                DTI=(-1.D0-XDEP(3))/(2.D0*X(3))
  533.                IND=1
  534.             ELSE
  535.                DTI=(1.D0-XDEP(3))/(2.D0*X(3))
  536.                IND=2
  537.             ENDIF
  538.             IF(TEST(DTINT,DTI))THEN
  539.                DTINT=DTI
  540.                DO 14 I1=1,IDIM-1
  541.                   XARI(I1)=XDEP(I1)-(X(I1)*DTINT)
  542.  14            CONTINUE
  543.                XARI(3)=ZI(IND)
  544.                ICONT=IND
  545.             ENDIF
  546.          ENDIF
  547. C
  548. C        ON CHERCHE SI ON A AFFAIRE    A UN NOEUD DU MAILLAGE
  549. C
  550.          IF((ABS(XARI(3)-1.D0).LT.EPS).OR.
  551.      *                 (ABS(XARI(3)+1.D0).LT.EPS))ITEST=ITEST+1
  552.          DO 41 I=1,2
  553.             IF(ABS(XARI(I)).LT.EPS)ITEST=ITEST+1
  554.             IF(ABS(1.D0-XARI(I)).LT.EPS)ITEST=ITEST+1
  555.  41       CONTINUE
  556. C
  557. C
  558. C ******************            TETRAEDRE
  559. C *               FACE 1 Y=0
  560. C *               FACE 2 Z=0
  561. C *               FACE 4 X=0
  562. C *               FACE 3 1-X-Y-Z=0
  563.       ELSEIF(ITYEL.EQ.23)THEN
  564.          IDIM=3
  565. C***  FACES 1  2 ET 4
  566. C
  567. C**                 ON TESTE L'INTERSECTION AVEC LA FACE 3 ET TOUTES LES
  568. C                        FACES AYANT UN FLUX POSITIF
  569. C
  570.          IF(UELEM(3).GT.0.D0)THEN
  571.             DTINT=(1-XDEP(3)-XDEP(2)-XDEP(1))/(2.D0*UELEM(3))
  572.             DO 56  I=1,IDIM
  573.                XARI(I)=XDEP(I)-(2.D0*UELEM(IFAT4(I))*DTINT)
  574.  56         CONTINUE
  575.             ICONT=3
  576.          ENDIF
  577. C
  578.          DO 55 I=1,IDIM
  579. C
  580.             IF(UELEM(IFAT4(I)).GT.0.D0)THEN
  581.                DTI=XDEP(I)/(2.D0*UELEM(IFAT4(I)))
  582. C
  583. C              ON NE GARDE QUE LE PLUS PETIT TEMPS POSITIF
  584. C
  585.                IF(TEST(DTINT,DTI))THEN
  586.                   DTINT=DTI
  587.                   I2=IJT4(1,I)
  588.                   I3=IJT4(2,I)
  589. C       ON GARDE LES COORDONNES CORRESPONDANT AUPLUS PETIT TEMPS
  590.                   XARI(I2)=XDEP(I2)-(2.D0*UELEM(IFAT4(I2))*DTINT)
  591.                   XARI(I3)=XDEP(I3)-(2.D0*UELEM(IFAT4(I3))*DTINT)
  592.                   XARI(I)=0.D0
  593.                   ICONT=IFAT4(I)
  594.                ENDIF
  595.             ENDIF
  596.  55      CONTINUE
  597. C
  598.  
  599. C         ON CHERCHE LE NUMERO DU NOEUD OU DE L'ARETE
  600. C
  601.          DO 51 I=1,IDIM
  602.             IF(ABS(XARI(I)).LT.EPS)ITEST=ITEST+1
  603.             IF(ABS(1.D0-XARI(I)).LT.EPS)ITEST=ITEST+1
  604.  
  605.  51      CONTINUE
  606.       ENDIF
  607. C
  608. C              TRAITEMENT DES CAS PARTICULIERS (NOEUDS/ARETES)
  609. C
  610.  
  611. C      SI LE POINT DE SORTIE EST UN NOEUD
  612.       IF(ITEST.EQ.IDIM)THEN
  613.          CALL NUMNOE(XARI,ITYEL,INOEU)
  614.       ELSE
  615. C       SINON ON VERIFIE EN 3D QU'ON TOMBE PAS SUR UNE ARETE
  616.          IF(IDIM.EQ.3)THEN
  617.             CALL ARTREF(XARI,ITYEL,IART)
  618.          ENDIF
  619.       ENDIF
  620.       END
  621.  
  622.  
  623.  
  624.  
  625.  
  626.  

© Cast3M 2003 - Tous droits réservés.
Mentions légales