Télécharger la plaquette Cast3M


fr en

Télécharger trjint.eso

Retour à la liste

Numérotation des lignes : 

trjint

  1. C TRJINT    SOURCE    CHAT      05/01/13    03:49:59     5004
  2.       SUBROUTINE TRJINT(XARI,XDEP,XINT,UELEM,DTCA,DTINT,
  3.      *                                              IDIM,ICONT,ITYEL,IO)
  4. C
  5. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  6. C
  7. C     INTERSECTION TRAJECTOIRE-ELEMENT  ( issu de TRIOEF)
  8. C
  9. C   ENTREES
  10. C     XDEP POSITION INITIALE
  11. C     XARI POSITION FINALE
  12. C     UELEM VITESSE
  13. C     DTCA PAS DE TEMPS
  14. C     IDIM DIMENSION DE L ESPACE
  15. C     ITYEL TYPE DE L ELEMENT
  16. C
  17. C   SORTIES
  18. C     XINT POSITION A L INTERSECTION
  19. C     DTINT TEMPS ECOULE JUSQU A LA SORTIE
  20. C     ICONT NUMERO DE LA PORTE DE SORTIE
  21. C     IO   DIFFERENT DE 0 ON EST SORTI DE L ELEMENT MAIS ON NE SAIT
  22. C          PAS PAR QUEL COTE
  23. C
  24. C   TYPES D ELEMENTS CONSIDERES
  25. C         TRI3 TRI6 TRI7 QUA4 QUA9 CUB8 PRI6
  26. C         4    6    7    8    11   14   16
  27. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  28. C
  29.       IMPLICIT INTEGER(I-N)
  30.       IMPLICIT REAL*8 (A-H,O-Z)
  31. C
  32.       DIMENSION IJ(2),ZI(2),IFAT3(3),IFAQ4(4),IFAC8(6)
  33.       DIMENSION IJKL(3,3),IJT4(2,3),IFAT4(4)
  34.       DIMENSION XARI(3),XDEP(3),XINT(3),UELEM(3),X(4)
  35. C
  36.       DATA IJKL/1,2,3,2,1,3,3,2,1/
  37.       DATA IJ /2,1/
  38.       DATA IJT4 /2,3, 1,3, 1,2/
  39.       DATA IFAT3 /3,1,2/
  40.       DATA IFAQ4 /4,2,1,3/
  41.       DATA IFAC8 /6,4,3,5,1,2/
  42.       DATA IFAT4 /4,1,2,3/
  43.       DATA ZI/-1.D0,1.D0/
  44. C
  45.       ICONT=0
  46.       IO=0
  47.       DTT=DTCA
  48.       DO 50 I=1,IDIM
  49.       XINT(I)=XDEP(I)
  50.    50 CONTINUE
  51. C     write(6,*)' inttrj ',XDEP,XARI,UELEM,DTCA,IDIM,ITYEL
  52. C
  53. C***  TRIANGLES TRI3 TRI6 TRI7             *
  54. C                                          * *
  55. C                               FACES    3 *  * 2
  56. C                                          *   *
  57. C                                          ******
  58. C                                            1
  59.       IF(ITYEL.EQ.4.OR.ITYEL.EQ.6.OR.ITYEL.EQ.7)THEN
  60.         IDIM=2
  61.         XARI(3)=1.D0-XARI(1)-XARI(2)
  62.         ITEST=0
  63.         DO 1 I=1,3
  64.           IF(XARI(I)*(1.D0-XARI(I)).GE.0.D0)ITEST=ITEST+1
  65.    1    CONTINUE
  66.         IF(ITEST.NE.3)THEN
  67. C            ON SORT DU TRIANGLE
  68. C
  69. C***  COTE 1 ET 3
  70. C
  71.           DO 2 I1=1,IDIM
  72.             IF(UELEM(I1).NE.0.AND.XDEP(I1).GT.0.D0)THEN
  73.               DTI=-XDEP(I1)/UELEM(I1)
  74.                   DTINT=DTI
  75.               IF(DTI*(DTT-DTI).GE.0.D0)THEN
  76.                 ITEST=0
  77.                 I3=IJ(I1)
  78.                 XINT(I3)=XDEP(I3)+UELEM(I3)*DTI
  79.                 X3=1.D0-XINT(I3)
  80.                 IF(XINT(I3)*(1.D0-XINT(I3)).GE.0.D0)ITEST=ITEST+1
  81.                 IF(X3*(1.D0-X3).GE.0.D0)ITEST=ITEST+1
  82.                 IF(ITEST.EQ.2)THEN
  83.                   ICONT=IFAT3(I1)
  84.                   XINT(I1)=0.D0
  85.                   DTINT=DTI
  86.                   RETURN
  87.                 ENDIF
  88.               ENDIF
  89.             ENDIF
  90.    2      CONTINUE
  91. C
  92. C***  COTE 2
  93. C
  94.           IF((UELEM(1)+UELEM(2)).NE.0.D0.AND.
  95.      *                      (1.D0-XDEP(1)-XDEP(2)).GT.0.D0)
  96.      *     THEN
  97.             DTI=(1.D0-XDEP(1)-XDEP(2))/(UELEM(1)+UELEM(2))
  98.                   DTINT=DTI
  99.             IF(DTI*(DTT-DTI).GE.0.D0)THEN
  100.               XINT(1)=XDEP(1)+UELEM(1)*DTI
  101.               XINT(2)=XDEP(2)+UELEM(2)*DTI
  102.               ITEST=0
  103.               IF(XINT(1)*(1.D0-XINT(1)).GE.0.D0)ITEST=ITEST+1
  104.               IF(XINT(2)*(1.D0-XINT(2)).GE.0.D0)ITEST=ITEST+1
  105.               IF(ITEST.EQ.2)THEN
  106.                 ICONT=2
  107.                 DTINT=DTI
  108.                 RETURN
  109.               ENDIF
  110.             ENDIF
  111.           ENDIF
  112. C
  113. C***  ??????????
  114. C
  115.           IO=3
  116.         ENDIF
  117. C
  118. C***  QUADRANGLE       QUA4 QUA9             3
  119. C                                           *****
  120. C                                           *   *
  121. C                                 FACES  4  *   * 2
  122. C                                           *****
  123. C                                             1
  124.       ELSEIF(ITYEL.EQ.8.OR.ITYEL.EQ.11)THEN
  125.         IDIM=2
  126.         ITEST=0
  127.         DO 3 I=1,IDIM
  128.           XY=(1.D0+XARI(I))*(1.D0-XARI(I))
  129.           IF(XY.GE.0.D0)ITEST=ITEST+1
  130.     3   CONTINUE
  131.         IF(ITEST.NE.IDIM)THEN
  132.           ID=IDIM-1
  133.           DO 4 I1=1,IDIM
  134.             IF(UELEM(I1).NE.0.D0)THEN
  135.               DO 5 I4=1,2
  136. C               IF(ZI(I4)*(ZI(I4)-XDEP(I1)).le.0.D0)
  137. C    * write(6,*)' xdep i1',i1,i4,xdep(i1)
  138.                 IF(ZI(I4)*(ZI(I4)-XDEP(I1)).GT.0.D0)THEN
  139.                   DTI=(ZI(I4)-XDEP(I1))/UELEM(I1)
  140.                   DTINT=DTI
  141.                   JI=I1
  142. C                 IF(DTI*(DTT-DTI).lt.0.D0)
  143. C    *   write(6,*)' dti dtt',dti,dtt
  144.                   IF(DTI*(DTT-DTI).GE.0.D0)THEN
  145.                     ITEST=0
  146.                     DO 6 I5=2,IDIM
  147.                       I7=IJKL(I5,I1)
  148.                       XINT(I7)=XDEP(I7)+UELEM(I7)*DTI
  149.                       XXX=(XINT(I7)+1.D0)*(1.D0-XINT(I7))
  150.                       IF(XXX.GE.0.D0)ITEST=ITEST+1
  151. C                write(6,*)' itest id ',itest,id,xxx
  152.    6                CONTINUE
  153.                     IF(ITEST.EQ.ID)THEN
  154.                       DTINT=DTI
  155.                       XINT(I1)=ZI(I4)
  156.                       ICONT=IFAQ4(2*I1+I4-2)
  157. C         WRITE(6,*)' SORTIE TRJINT ',XINT,DTINT,ICONT,IO
  158.                       RETURN
  159.                     ENDIF
  160.                   ENDIF
  161.                 ENDIF
  162.    5          CONTINUE
  163.             ENDIF
  164.    4      CONTINUE
  165. C
  166. C***  ??????????
  167. C
  168.           IO=4
  169.         ENDIF
  170. C
  171. C***   CUBE      CUB8
  172. C
  173.       ELSEIF(ITYEL.EQ.14)THEN
  174.         IDIM=3
  175.         ITEST=0
  176.         DO 13 I=1,IDIM
  177.           XY=(1.D0+XARI(I))*(1.D0-XARI(I))
  178.           IF(XY.GE.0.D0)ITEST=ITEST+1
  179.    13   CONTINUE
  180.         IF(ITEST.NE.IDIM)THEN
  181.           ID=IDIM-1
  182.           DO 14 I1=1,IDIM
  183.             IF(UELEM(I1).NE.0.D0)THEN
  184.               DO 15 I4=1,2
  185.                 IF(ZI(I4)*(ZI(I4)-XDEP(I1)).GT.0.D0)THEN
  186.                   DTI=(ZI(I4)-XDEP(I1))/UELEM(I1)
  187.                   DTINT=DTI
  188.                   JI=I1
  189.                   IF(DTI*(DTT-DTI).GE.0.D0)THEN
  190.                     ITEST=0
  191.                     DO 16 I5=2,IDIM
  192.                       I7=IJKL(I5,I1)
  193.                       XINT(I7)=XDEP(I7)+UELEM(I7)*DTI
  194.                       XXX=(XINT(I7)+1.D0)*(1.D0-XINT(I7))
  195.                       IF(XXX.GE.0)ITEST=ITEST+1
  196.   16                CONTINUE
  197.                     IF(ITEST.EQ.ID)THEN
  198.                       DTINT=DTI
  199.                       XINT(I1)=ZI(I4)
  200.                       ICONT=IFAC8(2*I1+I4-2)
  201.                       RETURN
  202.                     ENDIF
  203.                   ENDIF
  204.                 ENDIF
  205.   15          CONTINUE
  206.             ENDIF
  207.   14      CONTINUE
  208. C
  209. C***  ??????????
  210. C
  211.           IO=4
  212.         ENDIF
  213. C
  214. C***  PRISME    PRI6
  215. C
  216.       ELSEIF(ITYEL.EQ.16)THEN
  217.         IDIM=3
  218.         ITEST=0
  219.         X(1)=XARI(1)
  220.         X(2)=XARI(2)
  221.         X(3)=1.D0-X(1)-X(2)
  222.         X(4)=XARI(3)
  223.         DO 7 I=1,3
  224.           IF(X(I)*(1.D0-X(I)).GE.0.D0)ITEST=ITEST+1
  225.    7    CONTINUE
  226.         IF((X(4)+1.D0)*(1.D0-X(4)).GE.0.D0)ITEST=ITEST+1
  227.         IF(ITEST.NE.4)THEN
  228. C
  229. C***  FACE 5 ET 3
  230. C
  231.           DO 8 I1=1,2
  232.             IF(UELEM(I1).NE.0.D0.AND.XDEP(I1).GT.0.D0)THEN
  233.               DTI=-XDEP(I1)/UELEM(I1)
  234.                   DTINT=DTI
  235.               IF(DTI*(DTT-DTI).GE.0.D0)THEN
  236.                 ITEST=0
  237.                 I3=IJ(I1)
  238.                 XINT(I3)=XDEP(I3)+UELEM(I3)*DTI
  239.                 X3=1.D0-XINT(I3)
  240.                 XINT(3)=XDEP(3)+UELEM(3)*DTI
  241.                 IF(XINT(I3)*(1.D0-XINT(I3)).GE.0.D0)ITEST=ITEST+1
  242.                 IF((XINT(3)+1.D0)*(1.D0-XINT(3)).GE.0.D0)ITEST=ITEST+1
  243.                 IF(X3*(1.D0-X3).GE.0.D0)ITEST=ITEST+1
  244.                 IF(ITEST.EQ.3)THEN
  245.                   XINT(I1)=0.D0
  246.                   ICONT=IFAT3(I1)+2
  247.                   DTINT=DTI
  248.                   RETURN
  249.                 ENDIF
  250.               ENDIF
  251.             ENDIF
  252.    8      CONTINUE
  253. C
  254. C***  FACE 4
  255. C
  256.           IF((UELEM(1)+UELEM(2)).NE.0.D0.AND.
  257.      *                (1.D0-XDEP(1)-XDEP(2)).GT.0.D0)
  258.      *     THEN
  259.             DTI=(1.D0-XDEP(1)-XDEP(2))/(UELEM(1)+UELEM(2))
  260.                   DTINT=DTI
  261.             IF(DTI*(DTT-DTI).GE.0.D0)THEN
  262.               XINT(1)=XDEP(1)+UELEM(1)*DTI
  263.               XINT(2)=XDEP(2)+UELEM(2)*DTI
  264.               XINT(3)=XDEP(3)+UELEM(3)*DTI
  265.               ITEST=0
  266.               IF(XINT(1)*(1.D0-XINT(1)).GE.0.D0)ITEST=ITEST+1
  267.               IF(XINT(2)*(1.D0-XINT(2)).GE.0.D0)ITEST=ITEST+1
  268.               IF((1.D0-XINT(3))*(XINT(3)+1.D0).GE.0.D0)ITEST=ITEST+1
  269.               IF(ITEST.EQ.3)THEN
  270.                 ICONT=4
  271.                 DTINT=DTI
  272.                 RETURN
  273.               ENDIF
  274.             ENDIF
  275.           ENDIF
  276. C
  277. C***  FACE 1 ET 2
  278. C
  279.           DO 9 I2=1,2
  280.             IF(UELEM(3).NE.0.D0.AND.ZI(I2)*(ZI(I2)-XDEP(3))
  281.      *                                             .GT.0.D0)THEN
  282.               DTI=(ZI(I2)-XDEP(3))/UELEM(3)
  283.                   DTINT=DTI
  284.               IF(DTI*(DTT-DTI).GE.0.D0)THEN
  285.                 ITEST=0
  286.                 DO 10 I5=1,2
  287.                   XINT(I5)=XDEP(I5)+UELEM(I5)*DTI
  288.                   IF((1.D0-XINT(I5))*(XINT(I5)+1.D0).GE.0.D0)
  289.      *                                              ITEST=ITEST+1
  290.   10            CONTINUE
  291.                 X3=1.D0-XINT(1)-XINT(2)
  292.                 IF(X3*(1.D0-X3).GE.0.D0)ITEST=ITEST+1
  293.                 IF(ITEST.EQ.3)THEN
  294.                   XINT(3)=ZI(I2)
  295.                   ICONT=I2
  296.                   DTINT=DTI
  297.                   RETURN
  298.                 ENDIF
  299.               ENDIF
  300.             ENDIF
  301.    9      CONTINUE
  302. C
  303. C*** ?????????????
  304. C
  305.           IO=5
  306.         ENDIF
  307. C
  308. C ******************            TETRAEDRE
  309. C *               FACE 1 Y=0
  310. C *               FACE 2 Z=0
  311. C *               FACE 4 X=0
  312. C *               FACE 3 1-X-Y-Z=0
  313.       ELSEIF(ITYEL.EQ.23)THEN
  314.            ITEST=0
  315.            DO 20 I=1,3
  316.             IF(XARI(I)*(1.D0-XARI(I)).GE.0.D0)ITEST=ITEST+1
  317.    20      CONTINUE
  318.            XXX=XARI(1)+XARI(2)+XARI(3)
  319.             IF(XXX*(1.D0-XXX).GE.0.D0)ITEST=ITEST+1
  320. C          write(6,*)' trjint ', XXX,XARI(1),XARI(2),XARI(3),itest
  321.             IF(ITEST.NE.4)THEN
  322. C          write(6,*)' trjint ', XXX,XARI(1),XARI(2),XARI(3),itest
  323. C   ON SORT DU TETRAEDRE
  324. C
  325. C***  FACES 1  2 ET 4
  326. C
  327.           DO 25 I1=1,IDIM
  328. C          write(6,*)'axes',i1,uelem,xdep
  329.             IF(UELEM(I1).NE.0.AND.XDEP(I1).GT.0.D0)THEN
  330.               DTI=-XDEP(I1)/UELEM(I1)
  331.                   DTINT=DTI
  332.               IF(DTI*(DTT-DTI).GE.0.D0)THEN
  333.                 ITEST=0
  334.                 I2=IJT4(1,I1)
  335.                 I3=IJT4(2,I1)
  336.                 XINT(I2)=XDEP(I2)+UELEM(I2)*DTI
  337.                 XINT(I3)=XDEP(I3)+UELEM(I3)*DTI
  338.                 X3=1.D0-XINT(I2)-XINT(I3)
  339.                 IF(XINT(I2)*(1.D0-XINT(I2)).GE.0.D0)ITEST=ITEST+1
  340.                 IF(XINT(I3)*(1.D0-XINT(I3)).GE.0.D0)ITEST=ITEST+1
  341.                 IF(X3*(1.D0-X3).GE.0.D0)ITEST=ITEST+1
  342. C          write(6,*)'xint',xint,itest
  343.                 IF(ITEST.EQ.3)THEN
  344.                   ICONT=IFAT4(I1)
  345.                   XINT(I1)=0.D0
  346.                   DTINT=DTI
  347. C                 write(6,*)' sorti face ',icont,xint,io
  348.                   RETURN
  349.                 ENDIF
  350.               ENDIF
  351.             ENDIF
  352.    25     CONTINUE
  353. C
  354. C***  FACE 3
  355. C
  356.           IF((UELEM(1)+UELEM(2)+UELEM(3)).NE.0.D0.AND.
  357.      *                      (1.D0-XDEP(1)-XDEP(2)-XDEP(3)).GT.0.D0)
  358.      *     THEN
  359.             DTI=(1.D0-XDEP(1)-XDEP(2)-XDEP(3))/
  360.      *                                (UELEM(1)+UELEM(2)+UELEM(3))
  361.                   DTINT=DTI
  362.             IF(DTI*(DTT-DTI).GE.0.D0)THEN
  363.               XINT(1)=XDEP(1)+UELEM(1)*DTI
  364.               XINT(2)=XDEP(2)+UELEM(2)*DTI
  365.               XINT(3)=XDEP(3)+UELEM(3)*DTI
  366.               ITEST=0
  367.               IF(XINT(1)*(1.D0-XINT(1)).GE.0.D0)ITEST=ITEST+1
  368.               IF(XINT(2)*(1.D0-XINT(2)).GE.0.D0)ITEST=ITEST+1
  369.               IF(XINT(3)*(1.D0-XINT(3)).GE.0.D0)ITEST=ITEST+1
  370.               IF(ITEST.EQ.3)THEN
  371.                 ICONT=3
  372.                 DTINT=DTI
  373. C               XINT(3)=1.D0-XINT(1)-XINT(2)
  374. C                 write(6,*)' sorti face ',icont,xint,io
  375.                 RETURN
  376.               ENDIF
  377.             ENDIF
  378.           ENDIF
  379. C
  380. C***  ??????????
  381. C
  382.           IO=3
  383.         ENDIF
  384.       ENDIF
  385. C
  386.       END
  387.  
  388.  
  389.  

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