Sources Esope | Cast3M

trjav2
C TRJAV2    SOURCE    CB215821  20/11/25    13:41:26     10792          
      SUBROUTINE TRJAV2(IZVIT,IZPART,IZN3,IPART,TMIN,TMAX,
     *      MELEME,IELTFA,IZCENT,IFACEL,DELTAT,IZSH)
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C        Issu du sp TRJAVA  cas du calcul analytique
C
C    FAIT AVANCER UNE PARTICULE ( COORDONNEES DE REFERENCES )
C                 JUSQU'A L'INSTANT TMAX OU A L'INSTANT DE SORTIE
C    IZVIT SEGMENT DECRIVANT LES VITESSES ( &lt;--- TRJVIT TRJFLU)
C    IZPART POSITIONS INITIALES DES PARTICULES
C    IZN3  SEGMENT RESULTANT (AJUSTE ICI)
C    IPART NUMERO DE LA PARTICULE
C    TMIN  INSTANT DE DEPART
C    TMAX  INSTANT A NE PAS DEPASSER
C
C    MELEME POINTEUR DU MAILLAGE
C    IELTFA POINTEUR DE DOMAINE.ELTFA
C    IZCENT POINTEUR DE DOMAINE.CENTRE
C    IFACEL POINTEUR DE DOMAINE.FACEL
C    DELTAT PAS DE TEMPS AVEC LEQUEL ON CONSERVE LES RESULTATS
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
C
 
-INC PPARAM
-INC CCOPTIO
-INC SMCOORD
-INC SMELEME
-INC SMCHPOI
C
      POINTEUR IZCENT.MELEME,IELTFA.MELEME,IZFAC1.MELEME
      POINTEUR IFACEL.MELEME
C
      SEGMENT IZPART
           INTEGER NLEPA(NPART),NUMPA(NPART)
           REAL*8  COORPA(NDIM,NPART)
      ENDSEGMENT
C
      SEGMENT  IZN3
           INTEGER NAPAR3(NPOS),NUM3(NPOS)
           REAL*8 CREF3(NDIM,NPOS),TPAR(NPOS)
      ENDSEGMENT
C
      SEGMENT IZSH
            REAL*8 SHP(6,MNO9),SHY(12,MNO9),XYZL(3,MNO9)
      ENDSEGMENT
C
      SEGMENT IZVIT
           REAL*8 TEMTRA(NVIPT)
           INTEGER IPUN(NBS),IDUN(NBS),IPVPT(NVIPT),IFORML
      ENDSEGMENT
C             IDUN(I) nombre d elements avant le sous maillage I
C             IPVPT  pointeurs de izvpt pour chaque pas de temps
      SEGMENT IZVPT
           INTEGER IPUN1(NBS),IPUMAX
      ENDSEGMENT
C
      SEGMENT IZUN
           REAL*8 UN(I1,I2,I3)
      ENDSEGMENT
      POINTEUR IZUN1.IZUN,IZUN2.IZUN
C
      DIMENSION UELEM(3),XARI(3),XDEP(3),XINT(3)
C
C***
C
      NDIM=IDIM
      NPOS=50
      SEGINI IZN3
      SEGACT IZVIT
      IFORMU=IFORML
      IEL1=NLEPA(IPART)
C
      DO 2 ID1=1,NDIM
        XDEP(ID1)=COORPA(ID1,IPART)
 2    CONTINUE
C
      ITER=2
      NAPAR3(1)=IEL1
C
      DO 3 ID1=1,NDIM
         CREF3(ID1,1)=XDEP(ID1)
 3    CONTINUE
C
      TCOUR=TMIN
      TPAR(1)=TMIN
      NUAPAR=IEL1
C
      SEGACT IELTFA,IZCENT,IFACEL
      SEGACT IZSH
      NVIPT=TEMTRA(/1)
   1  CONTINUE
C       IVPT VAUT 1 EN PERMANENT
      IVPT=1
C     MELNEL PERMET D'AVOIR LE NUMERO DE L'ELEMENT DANS UN SOUS MAILLAGE
      CALL MELNEL(IEL1,MELEME,IPT1,NEL0,0)
      NOEL1=IPT1.NUM(/1)
      IELL=IEL1-NEL0
C     ON RECUPERE LES COORDONEES DES  NOEUDS DE L'ELEMENT
      CALL DOXE(XCOOR,IDIM,NOEL1,IPT1.NUM,IELL,XYZL)
      ITY1=IPT1.ITYPEL
      TTEMP=TCOUR
C            ON ACTIVE LE SEGMENT DES FLUX AUX FACES
      CALL TRJVEL(IZVIT,IZUN,IEL1,IVPT,TTEMP)
      ITYG=NUMGEO(ITY1)
C   ON DETERMINE LA TRAJECTOIRE LE TEMPS DE PARCOURS ET LA FACE DE SORTIE
      CALL TRJMA2(XARI,XDEP,UN(1,1,IELL),DTINT,NDIM,ICONT
     * ,ITYG,IART,INOEU)
 4    CONTINUE
C
C***   ON CALCULE LE JACOBIEN DE LA TRANSFORMATION EN XDEP
C
      CALL DETJAC(XYZL,XDEP,ITY1,NDIM,NOEL1,IFORMU,SHP,DET)
C
C***  POUR LES ELEMENTS NON TRIANGULAIRE LE JACOBIEN N'EST PAS CONSTANT
C
      IF(((ITYG-4)*(ITYG-23)).NE.0)THEN
C
C***  ON CALCULE PLUSIEURS JACOBIENS ET ON FAIT LA MOYENNE
C
         CALL DETJAC(XYZL,XARI,ITY1,NDIM,NOEL1,IFORMU,SHP,DET1)
         DTS=DTINT/10
         DO 12 ID1=1,9
C              CETTE ROUTINE PEUT ETRE FUSIONNEE AVEC TRJLDC
            CALL TRJSAU(XINT,XDEP,UN(1,1,IELL),DTS,NDIM,ITYG)
            CALL DETJAC(XYZL,XINT,ITY1,NDIM,NOEL1,IFORMU,SHP,DETINT)
            DET=DET+2*DETINT
            DO 14 ID2=1,NDIM
               XDEP(ID2)=XINT(ID2)
 14         CONTINUE
 12      CONTINUE
         DET=(DET+DET1)/20
       ENDIF
       DT=DTINT*DET
C
C***                 TRAITEMENT DES CAS PARTICULIERS
C
      IF(INOEU.NE.0)THEN
C
C***             RATRAPAGE A UN NOEUD DU MAILLLAGE
C
         TCOUR=TCOUR+DT
         CALL TRJRND(IZN3,IZVIT,MELEME,TCOUR,INOEU,IEL1,NDIM,ITER,
     *               IEL2,XARI,DTINT,ICONT,INOEU)
C         INOEU=INOE
         IEL1=IEL2
C       CAS OU ON SORT PAR UN NOEUD EN BORD DE MAILLAGE
         IF(DTINT.EQ.0.D0)THEN
            DT=DTINT
C        ON VA EN 10 STOCKER LES RESULTATS ET ON SORT
            GOTO 10
         ENDIF
         GOTO 4
C
C***                  RATRAPAGE A UNE ARETE
C
      ELSEIF(IART.NE.0)THEN
C
C***      ON CHERCHE LES FACES QUI SE COUPENT EN FORMANT L'ARETE
C
         CALL TRJARF(IART,NF1,NF2,ITYG)
C         ON CHOISIT UNE FACE AYANT UN FLUX SORTANT
         IF(UN(1,NF1,IELL).GT.0.D0)THEN
            ICONT=NF1
            GOTO 7
         ELSEIF(UN(1,NF2,IELL).GT.0.D0)THEN
            ICONT=NF2
            GOTO 7
C
C***     SINON ON EST EN BORD DE DOMAINE ON STOCKE LES RESULTATS EN 10
C
         ELSE
            DT=0.D0
            GOTO 10
         ENDIF
      ENDIF
C
C************* RECHERCHE DES VOISINS ET SCKAGE DES RESULTATS*********
C
C***  LA PARTICULE SORT DE L'ELEMENT COURANT ON CHERCHE LE VOISIN
C
 7    CONTINUE
        CALL TRJIEL(IEL1,IEL2,ICONT,NF,IFACEL,IZCENT,IELTFA)
C
C***   SI ON EST EN BORD DE DOMAINE ON ARRETE LE CALCUL
C
        IF(IEL2.EQ.0)GO TO 10
C
C***  ON ARRETE LE CALCUL PAR DEPASSEMENT DU TEMPS IMPOSE
C
       IF((TMAX-TCOUR).LE.0.D0)GOTO 21
C
C***  ON N'EST PAS EN BORD DE DOMAINE ET ON A TROUVE UN VOISIN
C
C     ON COMMENCE PAR LOCALISER LE NOUVEL ELEMENT
      CALL MELNEL(IEL2,MELEME,IPT2,NEL2,1)
      IELL2=IEL2-NEL2
      JTY=IPT2.ITYPEL
C     ON CHERCHE LA FACE COMMUNE ENTRE IEL1 ET IEL2
      CALL TRJFAC(IEL2,NF,JFA,IELTFA)
      IOR=1
      IF(NDIM.EQ.3)CALL TRJIOR(IELL,IELL2,IPT1,IPT2,ICONT,JFA,IOR)
      ITYG=NUMGEO(ITY1)
      JTYG=NUMGEO(JTY)
C     ON TRANSFORME LES COOR DE SORTIE DE IEL1 EN COOR D'ENTREE DANS IEL2
      CALL TRJTRJ(XDEP,XARI,ITYG,JTYG,JFA,ICONT,IOR)
      TCOUR=TCOUR+DT
      NUAPAR=IEL2
  13  CONTINUE
C
C***                 ON STOCKE LES RESULTATS
C
      DO 22 ID1=1,NDIM
        CREF3(ID1,ITER)=XDEP(ID1)
  22  CONTINUE
      TPAR(ITER)=TCOUR
      NAPAR3(ITER)=NUAPAR
      IEL1=NUAPAR
      IF ((NPOS-ITER).LE.1) THEN
        NPOS=NPOS+50
        SEGADJ IZN3
      ENDIF
      ITER=ITER+1
C
C***    ON RETOURNE EN 1 POUR ITERER LE PROCEDER AVEC UN NOUVEL ELEMENT
C
      GO TO 1
C
C***  LA PARTICULE I EST SORTIE DU DOMAINE DE CALCUL
C              ON STOCKE LES RESULTATS ET ON SORT
  10  CONTINUE
      DO 23 ID1=1,NDIM
        XDEP(ID1)=XARI(ID1)
  23  CONTINUE
      TCOUR=TCOUR+DT
      NUAPAR=IEL1
  21  CONTINUE
      DO 20 ID1=1,NDIM
        CREF3(ID1,ITER)=XDEP(ID1)
  20  CONTINUE
      TPAR(ITER)=TCOUR
      NAPAR3(ITER)=NUAPAR
      NPOS=ITER
      SEGADJ IZN3
C
c 999  FORMAT('PARTICULE',I4,': PROBLEMES DANS LE COIN D''UN ELEMENT !!',
c     *        2I7 )
C
      END