C LIEUPT    SOURCE    CHAT      05/01/13    01:17:33     5004
      SUBROUTINE LIEUPT(EPSILO,NDIM,ITYG,XREEL,IZSH,JTEST)

****************************************************************************
*** SP 'LIEUPT' : permet de determiner si pt considere est à l'interieur
*** ou à l'exterieur de l'element considere à epsilon pres
***
*** APPELES 1 = aucun
*** APPELES 2 = 'LIEUP2'
***
*** E = 'EPSILO' marge relative acceptée position pt % element
***     'NDIM' dimension de l'espace
***     'ITYG' entier caracterisant la geometrie de l'element considere
***     'XREEL' pt considere en coordonnees reelles
***     'IZSH' segment content coord reelles noeuds, fcts forme et base
***
*** S = 'JTEST' vaut 1 si pt à l'interieur de l'elemt a epsilon pres, 0 sinon
***
***              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 PT1(3),PT2(3),PT3(3),PT4(3),XREEL(3)
      DIMENSION ITRIAN(5),ICARRE(6),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,2/
      DATA ICARRE/1,2,3,4,1,2/
      DATA ICUBE/1,2,3,5, 5,6,7,1, 1,2,6,3, 2,3,7,1, 3,4,8,1, 1,4,8,2/
      DATA IPRISM/1,2,3,4, 4,5,6,1, 1,2,5,3, 2,3,6,1, 1,3,6,2/
      DATA ITETRA/1,2,4,3, 1,2,3,4, 2,3,4,1, 1,3,4,2/
***   on postule initialement 'XREEL' appartient à element
      JTEST=1

**************
*** CAS 2D ***
**************

***   cas TRI3 (triangles), boucle sur les différentes faces
      IF (ITYG.EQ.4) THEN
         DO 10 I=1,3
***         recuperation des pts definissant la face + l'origine
            DO 20 J=1,NDIM
               PT1(J)=XYZL(J,ITRIAN(I))
               PT2(J)=XYZL(J,ITRIAN(I+1))
               PT3(J)=XYZL(J,ITRIAN(I+2))
 20         CONTINUE
***         teste la position de 'XREEL' par rapport à la face
            CALL LIEUP2(NDIM,PT1,PT2,PT3,PT4,XREEL,TEST)
***          si 'XREEL' est du mauvais cote de la face
            IF (TEST.LT.(-EPSILO)) THEN
***            'JTEST'=0 car 'XREEL' est hors du TRI3
               JTEST=0
               RETURN
            ENDIF
 10      CONTINUE
***   cas QUA4 (quadrangles)
      ELSEIF (ITYG.EQ.8) THEN
         DO 30 I=1,4
            DO 40 J=1,NDIM
               PT1(J)=XYZL(J,ICARRE(I))
               PT2(J)=XYZL(J,ICARRE(I+1))
               PT3(J)=XYZL(J,ICARRE(I+2))
 40         CONTINUE
            CALL LIEUP2(NDIM,PT1,PT2,PT3,PT4,XREEL,TEST)
            IF (TEST.LT.(-EPSILO)) THEN
               JTEST=0
               RETURN
            ENDIF
 30      CONTINUE

**************
*** CAS 3D ***
**************

***   cas CUB8 (cubes)
      ELSEIF (ITYG.EQ.14) THEN
         DO 50 I=1,6
            DO 60 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))
 60         CONTINUE
            CALL LIEUP2(NDIM,PT1,PT2,PT3,PT4,XREEL,TEST)
            IF (TEST.LT.(-EPSILO)) THEN
               JTEST=0
               RETURN
            ENDIF
 50      CONTINUE
***   cas PRI6 (prismes)
      ELSEIF (ITYG.EQ.16) THEN
         DO 70 I=1,5
            DO 80 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))
 80         CONTINUE
            CALL LIEUP2(NDIM,PT1,PT2,PT3,PT4,XREEL,TEST)
            IF (TEST.LT.(-EPSILO)) THEN
               JTEST=0
               RETURN
            ENDIF
 70      CONTINUE
***   cas TET4 (tetraedres)
      ELSEIF (ITYG.EQ.23) THEN
         DO 90 I=1,4
            DO 100 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))
 100        CONTINUE
            CALL LIEUP2(NDIM,PT1,PT2,PT3,PT4,XREEL,TEST)
            IF (TEST.LT.(-EPSILO)) THEN
               JTEST=0
               RETURN
            ENDIF
 90      CONTINUE
      ENDIF
      RETURN
      END