Télécharger la plaquette Cast3M


fr en

Télécharger trjrat.eso

Retour à la liste

Numérotation des lignes : 

trjrat

  1. C TRJRAT    SOURCE    CB215821  20/11/25    13:41:38     10792          
  2.       SUBROUTINE TRJRAT(IEL1,IEL2,XDEP,XARI,IZCOU,MELEME,IZVIT,DT,
  3.      *       IZCENT,IELTFA,XINT,DTINT,ICONT,IZSH,TTEMP)
  4. C
  5. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  6. C
  7. C       APPELE PAR TRJAVA
  8. C     RATTRAPE LES COORDONNEES DE REFERENCES ET LES APPARTENANCES
  9. C       QUAND ON S'EST PERDU
  10. C     ( CE CAS SE PRESENTE LORSQUE LA TRAJECTOIRE PASSE PRES DES COINS)
  11. C
  12. C     IEL1 NUMERO DE L ELEMENT AVANT DE SE PERDRE
  13. C                              APRES
  14. C     XDEP COORDONNEES DE REFERENCES AU POINT DE PERDITION
  15. C                                    APRES
  16. C     IZVIT POINTEUR DU SEGMENT DECRIVANT LES VITESSES
  17. C     DT PAS DE TEMPS DU CALCUL
  18. C     IZCENT POINTEUR DU MAILLAGE CENTRE DES ELEMENTS
  19. C     IELTFA POINTEUR DES CONNECTIONS FACES ELEMENTS
  20. C     XINT COORDONNEES DE L'INTERSECTION FACE TRAJECTOIRE
  21. C     DTINT TEMPS ELEMENTAIRE AUQUEL LA TRAJECTOIRE TRAVERSE LA FACE
  22. C     ICONT NUMERO DE LA FACE PAR LAQUELLE ON SORT DE L'ELEMENT
  23. C     TTEMP TEMP AUQUEL EST FAIT LE CALCUL(utilisé en transitoire)
  24. C
  25. C          issu de RATTAR de TRIO-EF  adapté par F Auriol
  26. C
  27. CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
  28. C
  29.       IMPLICIT INTEGER(I-N)
  30.       IMPLICIT REAL*8 (A-H,O-Z)
  31. C
  32. C
  33.  
  34. -INC PPARAM
  35. -INC CCOPTIO
  36. -INC SMCOORD
  37. -INC SMELEME
  38. -INC SMCHPOI
  39. C
  40.       POINTEUR IZCENT.MELEME,IELTFA.MELEME
  41. C
  42.       SEGMENT IZPART
  43.            INTEGER NLEPA(NPART),NUMPA(NPART)
  44.            REAL*8  COORPA(NDIM,NPART)
  45.       ENDSEGMENT
  46.       POINTEUR IZL.IZPART,IZN.IZPART
  47.       SEGMENT  IZCOU
  48.            REAL*8 DTCO(NEL),COU
  49.       ENDSEGMENT
  50. C
  51.       SEGMENT IZSH
  52.             REAL*8 SHP(6,MNO9),SHY(12,MNO9),XYZL(3,MNO9)
  53.       ENDSEGMENT
  54. C
  55.       SEGMENT  IZNOEU
  56.             REAL*8 XELE(IDIM,NOEL)
  57.             INTEGER NOEGLO(NOEL)
  58.       ENDSEGMENT
  59.       SEGMENT IZVIT
  60.            REAL*8 TEMTRA(NVIPT)
  61.            INTEGER IPUN(NBS),IDUN(NBS),IPVPT(NVIPT),IFORML
  62.       ENDSEGMENT
  63. C             IDUN(I) nombre d elements avant le sous maillage I
  64. C             IPVPT  pointeurs de izvpt pour chaque pas de temps
  65.       SEGMENT IZVPT
  66.            INTEGER IPUN1(NBS),IPUMAX
  67.       ENDSEGMENT
  68.       SEGMENT IZUN
  69.            REAL*8 UN(I1,I2,I3)
  70.       ENDSEGMENT
  71. C
  72. C
  73.       DIMENSION XDEP(3),XARI(3),UE(3),XYL(3),XINT(3)
  74.       DIMENSION XYZ(3)
  75. C
  76. C
  77.       NDIM=IDIM
  78.       NPART=1
  79.       ICONT=0
  80.       NBNUL=0
  81.       SEGINI IZN,IZL
  82.       IZN.NLEPA(1)=0
  83.       IZN.NUMPA(1)=1
  84.       IZL.NLEPA(1)=0
  85.       IZL.NUMPA(1)=1
  86.       NEL0=0
  87.       IPT1=MELEME
  88.       NBSOUS=LISOUS(/1)
  89.       IF(NBSOUS.NE.0)THEN
  90.            DO 10 ISOU=1,NBSOUS
  91.                 IPT1=LISOUS(ISOU)
  92.                 SEGACT IPT1
  93.                 NEL=IPT1.NUM(/2)
  94.                 IF(NEL+NEL0.GE.IEL1)GO TO 15
  95.                 NEL0=NEL+NEL0
  96.                 SEGDES IPT1
  97.    10      CONTINUE
  98.       ENDIF
  99.    15 CONTINUE
  100.       ITYP=IPT1.ITYPEL
  101.       NOEL=IPT1.NUM(/1)
  102.       IELL=IEL1-NEL0
  103.       NVPT=IPVPT(/1)
  104.       IVPT=1
  105.       IF(NVPT.NE.1)THEN
  106.           IVPT=2
  107.           SEGACT IZVIT
  108.           CALL TRJTPT(IZVIT,TTEMP,IVPT)
  109.       ENDIF
  110.       CALL TRJVEL(IZVIT,IZUN,IEL1,IVPT,TTEMP)
  111.    19 CONTINUE
  112.       DO 18 I=1,3
  113.           XYL(I)=0.D0
  114.           UE(I)=0.D0
  115.           XYZ(I)=0.D0
  116.    18 CONTINUE
  117.       DTI=DT
  118.       DO 12 I=1,NDIM
  119.           XYZ(I)=XDEP(I)
  120.    12 CONTINUE
  121. C        ON CHERCHE SI ON A AFFAIRE    A UN NOEUD DU MAILLAGE
  122.        ITEST=0
  123.       CALL NOEREF(XYZ,ITYP,ITEST)
  124.       IF(ITEST.NE.0)THEN
  125. C                          CAS D'UN NOEUD DU MAILLAGE
  126. C          WRITE(6,*)'     CAS D UN NOEUD DU MAILLAGE ',XYZ,IEL1
  127.           IF(IFORML.EQ.1)THEN
  128.               CALL TRJVPO(UN(1,1,IELL),UE,XYZ,ITYP,NDIM,NOEL,SHP)
  129.               GO TO 50
  130.           ELSEIF(IFORML.EQ.2)THEN
  131.               SEGINI IZNOEU
  132.               CALL PREXN(IPT1,IELL,IZNOEU)
  133.               CALL VHYNOE(UE,IZNOEU,IZUN,NDIM,XYZ,ITYP,IELL)
  134.               SEGSUP IZNOEU
  135.               GO TO 50
  136.           ENDIF
  137.       ENDIF
  138.           IF(IFORML.EQ.1)THEN
  139.               CALL TRJVPO(UN(1,1,IELL),UE,XYZ,ITYP,NDIM,NOEL,SHP)
  140.               GO TO 50
  141.           ELSEIF(IFORML.EQ.2)THEN
  142. C               ON CHERCHE SI XINT EST UN NOEUD DU MAILLAGE
  143.           DO 22 I=1,NDIM
  144.           XYZ(I)=XINT(I)
  145.    22     CONTINUE
  146.           CALL NOEREF(XYZ,ITYP,ITEST)
  147.           IF(ITEST.NE.0)THEN
  148.               SEGINI IZNOEU
  149.               CALL PREXN(IPT1,IELL,IZNOEU)
  150.               CALL VHYNOE(UE,IZNOEU,IZUN,NDIM,XYZ,ITYP,IELL)
  151.               SEGSUP IZNOEU
  152.               DTI=DT-DTINT
  153.               GO TO 50
  154.           ENDIF
  155. C                ON TRAITE LE CAS OU XINT EST SUR UNE ARETE 3D
  156. C        ON IRA VERS LA FACE OU LE FLUX EST POSITIF
  157.           IF(IDIM.EQ.3)THEN
  158.              CALL ARTREF(XDEP,ITYP,IART)
  159.              IF(IART.NE.0)THEN
  160.                 CALL TRJARF(IART,NF1,NF2,ITYP)
  161.                 IF(UN(1,NF1,IELL).GT.0.D0)THEN
  162.                    ICONT=NF1
  163.                 ELSEIF(UN(1,NF2,IELL).GT.0.D0)THEN
  164.                    ICONT=NF2
  165.                 ELSE
  166.                   INTERR(2)= IEL1
  167.                   INTERR(3)= IART
  168.                   CALL ERREUR(-316)
  169.                    IEL2=0
  170.                    ICONT=0
  171.             IF(IIMPI.GT.0)THEN
  172.             WRITE(6,*) ' ON RESTE DANS L ELEMENT ',IEL1
  173.        WRITE(6,*)' ON EST SUR L ARETE ',IART ,'FLUX',UN(1,NF1,IELL),
  174.      *        UN(1,NF2,IELL)
  175.             ENDIF
  176.                    SEGSUP IZN,IZL
  177.                    RETURN
  178.                 ENDIF
  179.                 DO 48 ID1=1,NDIM
  180.                    XINT(ID1)=XDEP(ID1)
  181.    48           CONTINUE
  182.                 DTINT=0.D0
  183.                 SEGSUP IZN,IZL
  184.                 RETURN
  185.              ENDIF
  186.           ENDIF
  187. C                       FIN DU TRAITEMENT DES ARETES
  188. C
  189.           CALL FACREF(XDEP,ITYP,IFAC)
  190. C         WRITE(6,*)' ON EST SUR LA FACE ' ,IFAC ,IEL1,ITYP,XDEP
  191.               IF(IFAC.NE.0)THEN
  192. C                     ON EST SUR LA FACE IFAC
  193. C         WRITE(6,*)' ON EST SUR LA FACE ' ,IFAC ,UN(1,IFAC,IELL)
  194. C        LA PARTICULE SORTIRA DE L' ELEMENT SI LE FLUX EST POSITIF
  195.           IF(UN(1,IFAC,IELL).GT.0.D0)THEN
  196.               ICONT=IFAC
  197.               DO 45 ID1=1,NDIM
  198.                    XINT(ID1)=XDEP(ID1)
  199.    45         CONTINUE
  200.               DTINT=0.D0
  201.               SEGSUP IZN,IZL
  202.               RETURN
  203.           ELSEIF(UN(1,IFAC,IELL).LT.0.D0)THEN
  204. C               ON RESTE DANS L ELEMENT
  205. C  revoir ce cas s il se presente(a priori il est difficile a imaginer)
  206. C  Il se produit lorsque la vitesse s'annule tres pres de la face sur
  207. C  laquelle on  se trouve
  208. C      ( cas de forte difference de flux et orientation opposé  )
  209.             IF(IIMPI.GT.0)THEN
  210.             WRITE(6,*) ' ON RESTE DANS L ELEMENT ',IEL1
  211.        WRITE(6,*)' ON EST SUR LA FACE ',IFAC ,'FLUX',UN(1,IFAC,IELL)
  212.             ENDIF
  213.              INTERR(2)= IEL1
  214.              INTERR(3)= IFAC
  215.              CALL ERREUR(-315)
  216.              IEL2=0
  217.              ICONT=0
  218.              SEGSUP IZN,IZL
  219.              RETURN
  220.           ELSEIF(UN(1,IFAC,IELL).EQ.0.D0)THEN
  221. C                     CAS OU LA TRAJECTOIRE SUIT LA FACE
  222.             IF(IIMPI.GT.0)
  223.      *    WRITE(6,*) 'CAS OU LA TRAJECTOIRE SUIT LA FACE '
  224.           ENDIF
  225.           IEL2=IEL1
  226.           SEGSUP IZN,IZL
  227.           RETURN
  228.               ELSE
  229.        IF(IIMPI.GT.0)THEN
  230.        WRITE(6,*)' ON S EST PERDU AU POINT ',XDEP,' DANS L ELEMENT',
  231.      * IEL1
  232.        ENDIF
  233.        INTERR(2)= IEL1
  234.        CALL ERREUR(-317)
  235.              IEL2=0
  236.              ICONT=0
  237.        SEGSUP IZN,IZL
  238.        RETURN
  239.           ENDIF
  240.       ENDIF
  241.    50 CONTINUE
  242.       CALL SHAPE(XYZ(1),XYZ(2),XYZ(3),ITYP,SHP,IRET)
  243.       DO 20 J=1,NOEL
  244.           NN=(IPT1.NUM(J,IELL)-1)*(IDIM+1)
  245.           DO 25 K=1,NDIM
  246.               XYL(K)=XYL(K)+XCOOR(NN+K)*SHP(1,J)
  247.    25     CONTINUE
  248.    20 CONTINUE
  249.       DO 2 I=1,NDIM
  250.         IZL.COORPA(I,1)=XYL(I)+DTI*UE(I)
  251. C       write(6,*)' trjrat xyl ',IZL.COORPA(I,1),XYL(I),DTI,UE(I)
  252.    2  CONTINUE
  253.       CALL TRJPEL(IZL,IZN,MELEME,IZVIT,IZCOU,IZCENT,IELTFA,IZSH,TTEMP)
  254.       DO 3 I=1,NDIM
  255.         XARI(I)=IZN.COORPA(I,1)
  256.    3  CONTINUE
  257.       IEL2=IZN.NLEPA(1)
  258.       IF(IEL2.NE.0)THEN
  259. C                 ON VERIFIE QUE IEL2 A AU MOINS UN POINT COMMUN
  260. C       AVEC IEL1  SINON ON DIVISE LE PAS DE TEMPS PAR 2 (4 FOIS )
  261.       CALL MELNEL(IEL1,MELEME,IPT2,NEL2,1)
  262.       CALL MELNEL(IEL2,MELEME,IPT3,NEL3,1)
  263.       NBNN2=IPT2.NUM(/1)
  264.       NBNN3=IPT3.NUM(/1)
  265.       DO 80 I=1,NBNN2
  266.            DO 85 J=1,NBNN3
  267.                IF(IPT2.NUM(I,IEL1).EQ.IPT3.NUM(J,IEL2))GO TO 90
  268.    85      CONTINUE
  269.    80 CONTINUE
  270.       DT=DT*0.5
  271.       GO TO 19
  272.    90 CONTINUE
  273.       ELSE
  274.       NBNUL= NBNUL+1
  275.       IF(NBNUL.LE.4)THEN
  276.             DT=DT*0.5
  277.             GO TO 19
  278.       ELSE
  279.             INTERR(2)=IEL1
  280.             CALL ERREUR(-318)
  281.       ENDIF
  282.       ENDIF
  283.       SEGSUP IZN,IZL
  284. C
  285.       END
  286.  
  287.  
  288.  
  289.  
  290.  
  291.  

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