Télécharger la plaquette Cast3M


fr en

Télécharger trjrnd.eso

Retour à la liste

Numérotation des lignes : 

trjrnd

  1. C TRJRND    SOURCE    CB215821  20/11/25    13:41:39     10792          
  2.       SUBROUTINE TRJRND(IZN3,IZVIT,MELEME,TTEMP,NOE,IEL,NDIM,ITER,
  3.      *                  IELL2,XARI,DTI,ICONT,INOEU)
  4. C
  5. *******************************************************************
  6. *
  7. *   CE SOUS PROGRAMME CHERCHE UN VOISIN ET CALCULE LA TRAJECTOIRE
  8. *   LORSQU'ON TOMBE SUR UN NOEUD.
  9. *   ENTREES:
  10. *     IZN3  SEGMENT QUI STOCKE LES RESULTATS
  11. *     IZVIT SEGMENT QUI CONTIENT LES INFOS SUR LES FLUX
  12. *     MELEME MAILLAGE
  13. *     TTEMP  TEMPS COURANT PLUS LE PAS DE LA MAILLE IEL1
  14. *     INOEU  NOEUD PAR LEQUEL ON SORT DE IEL1
  15. *     IEL1   ELEMENT DUQUEL ON SORT
  16. *     NDIM
  17. *   SORTIES:
  18. *     IEL2 ELEMENT DANS LEQUEL PASSE LA TRAJECTOIRE
  19. *     XARI POSITION DE SORTIE DE L'ELEMENT
  20. *     DTINT TEMPS MIS POUR PARCOURIR IEL2
  21. *     ICONT FACE DE SORTIE DE IEL2
  22. *
  23. *******************************************************************
  24. C
  25.       IMPLICIT INTEGER(I-N)
  26.             IMPLICIT REAL*8 (A-H,O-Z)
  27. C
  28.  
  29. -INC PPARAM
  30. -INC CCOPTIO
  31. -INC SMCOORD
  32. -INC SMELEME
  33. -INC SMCHPOI
  34. C
  35.       SEGMENT  IZN3
  36.            INTEGER NAPAR3(NPOS),NUM3(NPOS)
  37.            REAL*8 CREF3(NDIM,NPOS),TPAR(NPOS)
  38.       ENDSEGMENT
  39. C
  40.       SEGMENT IZVIT
  41.            REAL*8 TEMTRA(NVIPT)
  42.            INTEGER IPUN(NBS),IDUN(NBS),IPVPT(NVIPT),IFORML
  43.       ENDSEGMENT
  44. C             IDUN(I) nombre d elements avant le sous maillage I
  45. C             IPVPT  pointeurs de izvpt pour chaque pas de temps
  46. C
  47.       SEGMENT IZUN
  48.            REAL*8 UN(I1,I2,I3)
  49.       ENDSEGMENT
  50.       POINTEUR IZUN1.IZUN,IZUN2.IZUN
  51. C
  52.       DIMENSION XARI(3),XDEP(3),X(4)
  53. C
  54.       DIMENSION XTRI3(2,3),XQUA4(2,4),XTET4(3,4)
  55.       DIMENSION XPRI6(3,6),XCUB8(3,8)
  56. C
  57.       DATA XTRI3 /0.D0,0.D0, 1.D0,0.D0, 0.D0,1.D0/
  58. C
  59.       DATA XQUA4 /-1.D0,-1.D0, 1.D0,-1.D0, 1.D0,1.D0, -1.D0,1.D0/
  60. C
  61.       DATA XTET4 /0.D0,0.D0,0.D0, 1.D0,0.D0,0.D0, 0.D0,1.D0,0.D0,
  62.      *            0.D0,0.D0,1.D0 /
  63. C
  64.       DATA XCUB8 /-1.D0,-1.D0,-1.D0, 1.D0,-1.D0,-1.D0, 1.D0,1.D0,-1.D0,
  65.      *            -1.D0,1.D0,-1.D0, -1.D0,-1.D0,1.D0, 1.D0,-1.D0,1.D0,
  66.      *             1.D0,1.D0,1.D0, -1.D0,1.D0,1.D0/
  67. C
  68.       DATA XPRI6 /0.D0,0.D0,-1.D0, 1.D0,0.D0,-1.D0, 0.D0,1.D0,-1.D0,
  69.      *            0.D0,0.D0,1.D0, 1.D0,0.D0,1.D0, 0.D0,1.D0,1.D0/
  70.  
  71. C
  72. C***                     RATTRAPAGE
  73. C       ON VA CHERCHER DANS LE MAILLAGE LES ELEMENTS QUI ONT UN
  74. C                 NOEUD EN COMMUN AVEC L'ELEMENT DE CALCUL
  75. C
  76. C***  ON CHERCHE DANS LE MAILLAGE UN ELEMENT IEL2 AYANT LE NOEUD
  77. C          D'ENTREE EN COMMUN AVEC LE NOEUD DE SORTIE DE IEL1
  78. C
  79. C                     INITIALISATION
  80.       NPAPAR=0
  81.       IEL2=0
  82.       IND=0
  83.       IVPT=1
  84.       ITEST=4
  85.       NPOS=CREF3(/2)
  86. C
  87. C***   ON ACTIVE LE MAILLAGE ET ON FAIT UNE RECHERCHE PAR SOUS MAILLAGE
  88. C      IPT5 EST LE SOUS DOMAINE CONTENANT LA MAILLE PAR LAQUELLE ON
  89. C      SORT PAR UN NOEUD
  90. C      IPT1 REPRESENTE LE SOUS MAILLAGE DANS LEQUEL ON RECHERCHE UN VOISIN
  91. C
  92.       SEGACT MELEME
  93.       NBSOUS=LISOUS(/1)
  94.       NBS=NBSOUS
  95.       IF(NBSOUS.EQ.0) NBS=1
  96.       IPT1=MELEME
  97.       CALL MELNEL(IEL,MELEME,IPT5,NEL5,1)
  98.       SEGACT IPT5
  99.       IELE=IEL-NEL5
  100.       DO 100 ISOUS=1,NBS
  101.          IF(NBSOUS.GT.0)THEN
  102.             IPT1=LISOUS(ISOUS)
  103.          ENDIF
  104.          SEGACT IPT1
  105.          NEL=IPT1.NUM(/2)
  106.          NBN=IPT1.NUM(/1)
  107. C
  108. C   POUR CHAQUE NOEUD DE CHAQUE ELEMENT ON REGARDE SI CE NOEUD
  109. C   CORRESPOND AU NOEUD DE L'ELEMENT COURANT PAR LEQUEL PASSE
  110. C   LA TRAJECTOIRE.DES QU'ON TROUVE UN ELEMENT AYANT LE NOEUD
  111. C   EN COMMUN ON CALCULE LA TRAJECTOIRE DE LA PARTICULE DANS CET
  112. C   ELEMENT.
  113. C
  114. C
  115. C            BOUCLE SUR LES ELEMENTS
  116.          DO 1 IEL2=1,NEL
  117. C            ON VA TESTER LES ELEMENTS DIFFERENTS DE L'ELEMENT COURANT
  118.             IF((IEL2.NE.IELE).OR.(IPT1.NE.IPT5))THEN
  119. C            BOUCLE SUR LES NOEUDS
  120.                DO 2 INO=1,NBN
  121. C            TEST DE COMPARAISON DES NOEUDS
  122.                   IF((IPT1.NUM(INO,IEL2)).EQ.(IPT5.NUM(NOE,IELE)))THEN
  123.                      IELL2=IEL2+IND
  124.                      ITEST=0
  125.                      ITYG=NUMGEO(IPT1.ITYPEL)
  126. C
  127. C        ON DETERMINE LES COORDONNES DU POINT D'ENTREE DE LA PARTICULE
  128. C        A PARTIR DU NOEUD TROUVE.ON CALCULE ENSUITE LA TRAJECTOIRE ET
  129. C        ON TESTE SI LES COORDONNEES DE SORTIE SONT VALABLES
  130. C
  131.                   CALL TRJVEL(IZVIT,IZUN,IELL2,IVPT,TTEMP)
  132. C
  133. C                            ON TESTE LES TRIANGLES
  134. C
  135.                   IF(ITYG.EQ.4)THEN
  136.                      DO 7 I=1,2
  137.                         XDEP(I)=XTRI3(I,INO)
  138.  7                   CONTINUE
  139.                      CALL TRJMA2(XARI,XDEP,UN(1,1,IEL2),DTI,2,
  140.      *                ICONT,4,IART,INOEU)
  141.                      XARI(3)=1.D0-XARI(1)-XARI(2)
  142.                      DO 5 I=1,3
  143.                         IF(XARI(I)*(1.D0-XARI(I)).GE.0.D0)ITEST=ITEST+1
  144.  5                      CONTINUE
  145.                      ITEST=ITEST-1
  146. C
  147.                   ELSEIF(ITYG.EQ.8)THEN
  148. C
  149. C                           ON TESTE LES QUADRANGLES
  150. C
  151.                      DO 8 I=1,2
  152.                         XDEP(I)=XQUA4(I,INO)
  153.  8                   CONTINUE
  154.                      CALL TRJMA2(XARI,XDEP,UN(1,1,IEL2),DTI,2,
  155.      *                ICONT,8,IART,INOEU)
  156.                      DO 3 I=1,NDIM
  157.                         XY=(1.D0+XARI(I))*(1.D0-XARI(I))
  158.                         IF(XY.GE.0.D0)ITEST=ITEST+1
  159.  3                   CONTINUE
  160. C
  161. C                     ON EST EN DIMENSION 3 ON TESTE LES CUBES
  162. C
  163.                   ELSEIF(ITYG.EQ.14)THEN
  164.                      DO 9 I=1,3
  165.                         XDEP(I)=XCUB8(I,INO)
  166.  9                   CONTINUE
  167.                      CALL TRJMA2(XARI,XDEP,UN(1,1,IEL2),DTI,3,
  168.      *                ICONT,14,IART,INOEU)
  169.                      DO 13 I=1,NDIM
  170.                         XY=(1.D0+XARI(I))*(1.D0-XARI(I))
  171.                         IF(XY.GE.0.D0)ITEST=ITEST+1
  172.  13                  CONTINUE
  173. C
  174. C                        ON TESTE LES PRISMES
  175.                   ELSEIF(ITYG.EQ.16)THEN
  176. C
  177.                      DO 10 I=1,3
  178.                         XDEP(I)=XPRI6(I,INO)
  179.  10                  CONTINUE
  180.                      CALL TRJMA2(XARI,XDEP,UN(1,1,IEL2),DTI,3,
  181.      *                ICONT,16,IART,INOEU)
  182.                      X(1)=XARI(1)
  183.                      X(2)=XARI(2)
  184.                      X(3)=1.D0-X(1)-X(2)
  185.                      X(4)=XARI(3)
  186.                      DO 14 I=1,3
  187.                         IF(X(I)*(1.D0-X(I)).GE.0.D0)ITEST=ITEST+1
  188.  14                  CONTINUE
  189.                      IF((X(4)+1.D0)*(1.D0-X(4)).GE.0.D0)ITEST=ITEST+1
  190.                      ITEST=ITEST-1
  191. C
  192. C                        ON TESTE LES TETRAEDRES
  193.                   ELSEIF(ITYG.EQ.23)THEN
  194. C
  195.                      DO 11 I=1,3
  196.                         XDEP(I)=XTET4(I,INO)
  197.  11                  CONTINUE
  198.                      CALL TRJMA2(XARI,XDEP,UN(1,1,IEL2),DTI,3,
  199.      *                ICONT,23,IART,INOEU)
  200.                      DO 20 I=1,3
  201.                         IF(XARI(I)*(1.D0-XARI(I)).GE.0.D0)ITEST=ITEST+1
  202.  20                  CONTINUE
  203.                      XXX=XARI(1)+XARI(2)+XARI(3)
  204.                      IF(XXX*(1.D0-XXX).GE.0.D0)ITEST=ITEST+1
  205.                      ITEST=ITEST-1
  206. C
  207.                   ENDIF
  208. C
  209. C            UNE FOIS QUE L'ON A CALCULER LA TRAJECTOIRE POUR UN ELEMENT
  210. C               ON TESTE SI LA TRAJECTOIRE APPARTIENT A CET ELEMENT
  211. C                   DANS CE CAS ON STOCKE LES VALEURS A CONSERVER
  212. C
  213.                   IF((ITEST.EQ.NDIM).AND.(INO.NE.INOEU))THEN
  214.                      TPAR(ITER)=TTEMP
  215.                      DO 4 ID1=1,NDIM
  216.                         CREF3(ID1,ITER)=XDEP(ID1)
  217.  4                   CONTINUE
  218.                      NAPAR3(ITER)=IELL2
  219.                      IF ((NPOS-ITER).LE.1) THEN
  220.                         NPOS=NPOS+50
  221.                         SEGADJ IZN3
  222.                         ITER=ITER+1
  223.                      ENDIF
  224.                      ITER=ITER+1
  225. C       SI ON A TROUVE LE BON ELEMENT ON SORT
  226.                      GOTO 17
  227.                   ENDIF
  228. C
  229. C   SINON ON RECOMMENCE AVEC UN AUTRE NOEUD OU UN AUTRE ELEMENT
  230. C
  231.                ENDIF
  232.  2       CONTINUE
  233.             ENDIF
  234.  1    CONTINUE
  235.       IND=IND+NEL
  236.       IF(IPT1.NE.IPT5)SEGDES IPT1
  237.  100  CONTINUE
  238.       SEGDES IPT5
  239. C
  240. C               CAS OU ON ARRIVE EN BORD DE DOMAINE
  241. C         SI ON TROUVE PAS DE VOISIN ON RETOURNE DANS TRJAVA
  242.           DTI=0.D0
  243.           IF(ITEST.EQ.4)GOTO 17
  244. C         CAS OU LA TRAJECTOIRE NE PASSE PAS PAR LE VOISIN TROUVE
  245.           DO 18 I=1,NDIM
  246.              XARI(I)=XDEP(I)
  247.  18       CONTINUE
  248.  17   CONTINUE
  249.       RETURN
  250.       END
  251.  
  252.  
  253.  
  254.  
  255.  
  256.  
  257.  
  258.  
  259.  
  260.  

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