Sources Esope | Cast3M

factra
C FACTRA    SOURCE    CHAT      05/01/12    23:56:02     5004
      SUBROUTINE FACTRA(NDIM,ITYG,XDEP2,XARI2,IZSH,IFACE,XINT,XN,NBFAC)
*************************************************************************
*** SP 'FACTRA' : par convection ou diffusion explicite (où traj=droite),
*** donne n° faces traversees, pts intersections associés, normales du
*** repere local à la face associées.
***
*** APPELES 1 = aucun
*** APPELES 2 = 'FACTR1', 'FACTR2'
***
*** E = 'NDIM' dimension de l'espace
***     'ITYG' entier caracterisant la geometrie de l'element
***     'XDEP2' coordonnees reelles de depart de la particule
***     'XARI2' coordonnees reelles d'arrivee de la particule
***     'IZSH' segment content coords reelles des noeuds de l'elemt considere
***
*** S = 'IFACE' n° des faces traversees par particule
***     'XINT' pts intersection trajectoire particule-faces traversees
***     'XN' vecteurs unitaires normaux aux faces traversees
***     'NBFAC' nbre de faces de l'element considere traversees par particule
***
***       Auteur Cyril Nou
*************************************************************************
 
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8 (A-H,O-Z)
      SEGMENT IZSH
          REAL*8 SHP(6,MNO9),SHY(12,MNO9),XYZL(3,MNO9)
      ENDSEGMENT
      DIMENSION XDEP2(3),XARI2(3),IFACE(6),XINT(3,6),XN(3,6)
      DIMENSION XN1(3),XN2(3),XN3(3),XINTER(3)
      DIMENSION PT1(3),PT2(3),PT3(3),PT4(3)
      DIMENSION ITRIAN(4),ICARRE(5),ICUBE(4,6),IPRISM(4,5),ITETRA(4,4)
***   donnees ordonnées specifiquement pour chaque type d'element afin de
***   parcourir les faces dans l'ordre croissant dans les diverses boucles
      DATA ITRIAN/1,2,3,1/
      DATA ICARRE/1,2,3,4,1/
      DATA ICUBE/1,2,3,4, 5,6,7,8, 1,2,6,5, 2,3,7,6, 3,4,8,7, 1,4,8,5/
      DATA IPRISM/1,2,3,3, 4,5,6,6, 1,2,5,4, 2,3,6,5, 1,3,6,4/
      DATA ITETRA/1,2,4,4, 1,2,3,3, 2,3,4,4, 1,3,4,4/
***   initialisation des arguments de sortie à 0
      DO 10 I=1,6
         IFACE(I)=0
         DO 20 J=1,NDIM
            XINT(J,I)=0.D0
            XN(J,I)=0.D0
 20      CONTINUE
 10   CONTINUE
      NBFAC=0
 
***   cas TRI3 (triangles), boucle sur les différentes faces
      IF (ITYG.EQ.4) THEN
         DO 30 I=1,3
***         recuperation des pts definissant la face
            DO 35 J=1,NDIM
               PT1(J)=XYZL(J,ITRIAN(I))
               PT2(J)=XYZL(J,ITRIAN(I+1))
 35         CONTINUE
***         recherche des normales et du pt d'intersection
            CALL FACTR1(NDIM,PT1,PT2,PT3,XDEP2,XARI2,XN1,XN2,XN3
     $           ,XINTER,ITEST)
***         'ITEST'=1 s'il y a intersection avec plan associé à la face
            IF (ITEST.EQ.1) CALL FACTR2(NDIM,ITYG,PT1,PT2,PT3,PT4
     $           ,XDEP2,XARI2,XN1,XN2,XN3,XINTER,ITEST)
***         'ITEST'=1 s'il y a intersection avec la face
            IF (ITEST.EQ.1) THEN
***            recuperation des infos de la face effectivement traversée
               NBFAC=NBFAC+1
               IFACE(NBFAC)=I
               DO 40 K=1,NDIM
                  XINT(K,NBFAC)=XINTER(K)
                  XN(K,NBFAC)=XN1(K)
 40            CONTINUE
            ENDIF
 30      CONTINUE
***   cas QUA4 (quadrangles)
      ELSEIF (ITYG.EQ.8) THEN
         DO 50 I=1,4
            DO 55 J=1,NDIM
               PT1(J)=XYZL(J,ICARRE(I))
               PT2(J)=XYZL(J,ICARRE(I+1))
 55         CONTINUE
            CALL FACTR1(NDIM,PT1,PT2,PT3,XDEP2,XARI2,XN1,XN2,XN3
     $           ,XINTER,ITEST)
            IF (ITEST.EQ.1) CALL FACTR2(NDIM,ITYG,PT1,PT2,PT3,PT4
     $           ,XDEP2,XARI2,XN1,XN2,XN3,XINTER,ITEST)
            IF (ITEST.EQ.1) THEN
               NBFAC=NBFAC+1
               IFACE(NBFAC)=I
               DO 60 K=1,NDIM
                  XINT(K,NBFAC)=XINTER(K)
                  XN(K,NBFAC)=XN1(K)
 60            CONTINUE
            ENDIF
 50      CONTINUE
***   cas CUB8 (cubes)
      ELSEIF (ITYG.EQ.14) THEN
         DO 70 I=1,6
            DO 75 J=1,NDIM
               PT1(J)=XYZL(J,ICUBE(1,I))
               PT2(J)=XYZL(J,ICUBE(2,I))
               PT3(J)=XYZL(J,ICUBE(3,I))
               PT4(J)=XYZL(J,ICUBE(4,I))
 75         CONTINUE
            CALL FACTR1(NDIM,PT1,PT2,PT3,XDEP2,XARI2,XN1,XN2,XN3
     $           ,XINTER,ITEST)
            IF (ITEST.EQ.1) CALL FACTR2(NDIM,ITYG,PT1,PT2,PT3,PT4
     $           ,XDEP2,XARI2,XN1,XN2,XN3,XINTER,ITEST)
            IF (ITEST.EQ.1) THEN
               NBFAC=NBFAC+1
               IFACE(NBFAC)=I
               DO 80 K=1,NDIM
                  XINT(K,NBFAC)=XINTER(K)
                  XN(K,NBFAC)=XN1(K)
 80            CONTINUE
            ENDIF
 70      CONTINUE
***   cas PRI6 (prismes)
      ELSEIF (ITYG.EQ.16) THEN
         DO 90 I=1,5
            DO 95 J=1,NDIM
               PT1(J)=XYZL(J,IPRISM(1,I))
               PT2(J)=XYZL(J,IPRISM(2,I))
               PT3(J)=XYZL(J,IPRISM(3,I))
               PT4(J)=XYZL(J,IPRISM(4,I))
 95         CONTINUE
            CALL FACTR1(NDIM,PT1,PT2,PT3,XDEP2,XARI2,XN1,XN2,XN3
     $           ,XINTER,ITEST)
            IF (ITEST.EQ.1) CALL FACTR2(NDIM,ITYG,PT1,PT2,PT3,PT4
     $           ,XDEP2,XARI2,XN1,XN2,XN3,XINTER,ITEST)
            IF (ITEST.EQ.1) THEN
               NBFAC=NBFAC+1
               IFACE(NBFAC)=I
               DO 100 K=1,NDIM
                  XINT(K,NBFAC)=XINTER(K)
                  XN(K,NBFAC)=XN1(K)
 100           CONTINUE
            ENDIF
 90      CONTINUE
***   cas TET4 (tetraedres)
      ELSEIF (ITYG.EQ.23) THEN
         DO 110 I=1,4
            DO 115 J=1,NDIM
               PT1(J)=XYZL(J,ITETRA(1,I))
               PT2(J)=XYZL(J,ITETRA(2,I))
               PT3(J)=XYZL(J,ITETRA(3,I))
               PT4(J)=XYZL(J,ITETRA(4,I))
 115        CONTINUE
            CALL FACTR1(NDIM,PT1,PT2,PT3,XDEP2,XARI2,XN1,XN2,XN3
     $           ,XINTER,ITEST)
            IF (ITEST.EQ.1) CALL FACTR2(NDIM,ITYG,PT1,PT2,PT3,PT4
     $           ,XDEP2,XARI2,XN1,XN2,XN3,XINTER,ITEST)
            IF (ITEST.EQ.1) THEN
               NBFAC=1
               IFACE(NBFAC)=I
               DO 120 K=1,NDIM
                  XINT(K,NBFAC)=XINTER(K)
                  XN(K,NBFAC)=XN1(K)
 120           CONTINUE
            ENDIF
 110     CONTINUE
      ENDIF
      RETURN
      END