C CRCOU2    SOURCE    PV        21/11/04    21:15:06     11174          
C  CALCUL D'UN OBJET COUPE VIRTUEL  (CAS DE DEFORMEES)
C
C  CRCOU2 REMPLACE MELEME PAR L'OBJET COUPE
C  MCOUP ET LE DERNIER COMPOSANT DE MELEME REMPLACE DECRIVENT
C  L'INTERSECTION
C  IL FAUDRA METTRE UNE INDICATION AUX TRIANGLES SE TROUVANT DANS LE PLA
C  DE COUPE POUR REUSSIR LES PARTIES VUES PARTIES CACHEES
C
C SG 2016/07/18 : traitement du cas des faces tri7/qua9
C
      SUBROUTINE CRCOU2(JOEIL,ICOUP1,ICOUP2,ICOUP3,MELEME,MCOUP,VCHCA,
     # SXCORD,ICPR,MELEM2,MCOU2,ITE,IVU,mchamor,isect)
      IMPLICIT INTEGER(I-N)
-INC CCGEOME
-INC CCREEL
-INC PPARAM
-INC CCOPTIO
-INC SMCOORD
-INC SMELEME
-INC SMCHAML
      SEGMENT SXCORD
       real*8   XCORD(IDIM,ITE)
      ENDSEGMENT
      REAL*8 XN,YN,ZN,SIGA,SIGB,SIGC,SIGRAV
* sg
      REAL*8 xo,yo,zo,xct,yct,zct,xc2,yc2,zc2,xc3,yc3,zc3
      REAL*8 xv,yv,zv,xw,yw,zw 
      REAL*8 rn,sig,chgrav
* sg
      SEGMENT IVU(0)
      SEGMENT VCHCA(0)
      SEGMENT ITR(0)
      SEGMENT ITRT(0)
      SEGMENT XTR(0)
      SEGMENT XTRT(0)
      SEGMENT IQUAT(0)
      SEGMENT XQUAT(0)
      SEGMENT ISEGM(0)
      SEGMENT XSEGM(NBPOIN)
      SEGMENT MCOUP(0)
      SEGMENT ICPR(0)
      SEGMENT IJGRAV(2,100)
      DIMENSION JTT(9),ISEG(9)
      LOGICAL VOLUM
      logical lquaf
*
*dbg      write(ioimp,*) 'coucou crcou2 isect,vchca,xcord=',isect,vchca
*dbg     $     ,xcord
      n2ptel=0
      n2el=0
      mcham=mchamor
      MEDSAU=0
      SEGINI ITR,ISEGM,MCOUP,IJGRAV,ITRT,IQUAT
      SEGACT SXCORD,MCOORD
      NBPOIN=nbpts
      if (mcham.ne.0) segini xtr,xtrt,xquat,xsegm
      ITE=XCORD(/2)
      IREF=(JOEIL-1)*4
      XO=XCOOR(IREF+1)
      YO=XCOOR(IREF+2)
      ZO=XCOOR(IREF+3)
*dbg      WRITE(IOIMP,*) 'joeil,x,y,z=',joeil,xo,yo,zo
      IREF=(ICOUP1-1)*4
      XCT=XCOOR(IREF+1)
      YCT=XCOOR(IREF+2)
      ZCT=XCOOR(IREF+3)
*dbg      WRITE(IOIMP,*) 'icoup1,x,y,z=',icoup1,xct,yct,zct
      IREF=(ICOUP2-1)*4
      XC2=XCOOR(IREF+1)
      YC2=XCOOR(IREF+2)
      ZC2=XCOOR(IREF+3)
*dbg      WRITE(IOIMP,*) 'icoup2,x,y,z=',icoup2,xc2,yc2,zc2
      IREF=(ICOUP3-1)*4
      XC3=XCOOR(IREF+1)
      YC3=XCOOR(IREF+2)
      ZC3=XCOOR(IREF+3)
*dbg      WRITE(IOIMP,*) 'icoup3,x,y,z=',icoup3,xc3,yc3,zc3
C  NORMALE AU PLAN :
      XV=XC2-XCT
      YV=YC2-YCT
      ZV=ZC2-ZCT
      XW=XC3-XCT
      YW=YC3-YCT
      ZW=ZC3-ZCT
      XN=YV*ZW-ZV*YW
      YN=ZV*XW-XV*ZW
      ZN=XV*YW-YV*XW
      RN=XN**2+YN**2+ZN**2
*      IF (RN.EQ.0.) CALL ERREUR(21)
      IF (ABS(RN).LT.(XPETIT**0.25D0)) CALL ERREUR(21)
      IF (IERR.NE.0) RETURN
      RN=SQRT(RN)
      XN=XN/RN
      YN=YN/RN
      ZN=ZN/RN
      SIG=(XO-XCT)*XN+(YO-YCT)*YN+(ZO-ZCT)*ZN
      IF (SIG.LE.0.) THEN
         XN=-XN
         YN=-YN
         ZN=-ZN
      ENDIF
*  DANS IVU ON NOTE LES POINTS DU MAUVAIS COTE DU PLAN DE COUPE
      DO 60 I=1,XCORD(/2)
         XP=XCORD(1,I)
         YP=XCORD(2,I)
         ZP=XCORD(3,I)
         SIG=(XP-XCT)*XN+(YP-YCT)*YN+(ZP-ZCT)*ZN
         IF (SIG.GT.0.) IVU(I)=0
*dbg         WRITE(IOIMP,*) 'I,X,Y,Z,S,IVU=',I,XP,YP,ZP,SIG,IVU(I)
 60   CONTINUE
 5000 CONTINUE
      SEGACT MELEME
      NBSOU=LISOUS(/1)
      NBSOUS=MAX(2,NBSOU+1)
      NBREF=0
      if (mcham.ne.0) nbref=nbsous
      NBNN=0
      NBELEM=0
      SEGINI IPT8
      IPT1=MELEME
      ISU=0
*  prob optimiseur
      melva1=meleme
      melva2=meleme
      DO 10 ISOUS=1,MAX(1,NBSOU)
         IF (NBSOU.NE.0) THEN
            IPT1=LISOUS(ISOUS)
            if (mcham.ne.0) then
               melva1=lisref(isous)
               segact melva1
               lva1=melva1.velche(/1)
            endif
            SEGACT IPT1
         ENDIF
         NBREF=0
         NBSOUS=0
         NBNN=IPT1.NUM(/1)
         NBELEM=IPT1.NUM(/2)
         SEGINI IPT2
         IPT2.ITYPEL=IPT1.ITYPEL
         if (mcham.ne.0) then
            n1ptel=nbnn
            n1el=nbelem
            segini melva2
         endif
         JELN=0
         DO 20 JEL=1,NBELEM
            DO 25 INOEU=1,NBNN
               IPT=IPT1.NUM(INOEU,JEL)
               IF (IVU(ICPR(IPT)).EQ.0) GOTO 30
 25         CONTINUE
            if (isect.ne.0) goto 20
            JELN=JELN+1
*     write (6,*) ' mcham melva1 melva2 meleme ',mcham,melva1,melva2,
*    >    meleme
            DO INOEU=1,NBNN
               IPT2.NUM(INOEU,JELN)=IPT1.NUM(INOEU,JEL)
            ENDDO
            if (mcham.ne.0) then
               DO INOEU=1,NBNN
                  melva2.velche(INOEU,JELN)=
     $                 melva1.velche(min(INOEU,lva1),JEL)
               ENDDO
            ENDIF
            IPT2.ICOLOR(JELN)=IPT1.ICOLOR(JEL)
            GOTO 20
 30         CONTINUE
            VOLUM=KSURF(IPT1.ITYPEL).NE.IPT1.ITYPEL
*     IF (VOLUM) THEN
            IPT=IPT1.NUM(1,JEL)
            XGRAV=XCORD(1,ICPR(IPT))
            YGRAV=XCORD(2,ICPR(IPT))
            ZGRAV=XCORD(3,ICPR(IPT))
            IF (VCHCA.NE.0) VGRAV=VCHCA(IPT)
            if (mcham.ne.0) chgrav=melva1.velche(1,JEL)
            DO 21 INOEU=2,NBNN
               IPT=IPT1.NUM(INOEU,JEL)
               XGRAV=XGRAV+XCORD(1,ICPR(IPT))
               YGRAV=YGRAV+XCORD(2,ICPR(IPT))
               ZGRAV=ZGRAV+XCORD(3,ICPR(IPT))
               IF (VCHCA.NE.0) VGRAV=VGRAV+VCHCA(IPT)
               if (mcham.ne.0) chgrav=chgrav+
     $              melva1.velche(min(INOEU,lva1),JEL)
 21         CONTINUE
            XGRAV=XGRAV/NBNN
            YGRAV=YGRAV/NBNN
            ZGRAV=ZGRAV/NBNN
C  ON PREND COMME POINT DU PLAN LA PROJECTION DU CENTRE GRAVITE
            SIGRAV=(XGRAV-XCT)*XN+(YGRAV-YCT)*YN+(ZGRAV-ZCT)*ZN
            XCT=XGRAV-SIGRAV*XN
            YCT=YGRAV-SIGRAV*YN
            ZCT=ZGRAV-SIGRAV*ZN
            SIGRAV=(XGRAV-XCT)*XN+(YGRAV-YCT)*YN+(ZGRAV-ZCT)*ZN
            IF (ABS(SIGRAV).LE.1D-30) SIGRAV=1D-30
            NJGRAV=0
            IF (VCHCA.NE.0) VGRAV=VGRAV/NBNN
            if (mcham.ne.0) CHGRAV=CHGRAV/NBNN
*     ENDIF
C  DECOMPOSITION EN FACE PUIS EN TRIANGLE
            NBFAC=LTEL(1,IPT1.ITYPEL)
            IAD=LTEL(2,IPT1.ITYPEL)-1
            IF (NBFAC.EQ.0) GOTO 20
            DO 161 IFAC=1,NBFAC
               ITYP=LDEL(1,IAD+IFAC)
* SG  Pour les TRI7/QUA9, il y a un traitement particulier pour ne pas
*     tracer les segments relies au noeud central. Il prend en compte le
*     fait que le decoupage est barycentrique
         lquaf=(ityp.eq.7.OR.ityp.eq.8) 
         JAD=LDEL(2,IAD+IFAC)-1
         NPFAC=KDFAC(1,ITYP)
         IDEP=KDFAC(2,ITYP)
         IFEP=IDEP+3*(KDFAC(3,ITYP)-1)
         DO 160 ITRIAN=IDEP,IFEP,3
            kfa=kfac(itrian)
            kfb=kfac(itrian+1)
            kfc=kfac(itrian+2)
            JFA=JAD+KFAC(ITRIAN)
            JFB=JAD+KFAC(ITRIAN+1)
            JFC=JAD+KFAC(ITRIAN+2)
            ISEG(1)=0
            ISEG(2)=0
            ISEG(3)=0
            ISEG(4)=0
            ISEG1=0
            ISEG2=0
            ISEG3=0
            if (lquaf) then
               if (KFA.NE.NPFAC.AND.KFB.NE.NPFAC) ISEG1=1
               if (KFB.NE.NPFAC.AND.KFC.NE.NPFAC) ISEG2=1
               if (KFA.NE.NPFAC.AND.KFC.NE.NPFAC) ISEG3=1
            else
               IF (MOD(JFB-JFA+NPFAC,NPFAC).EQ.1) ISEG1=1
               IF (MOD(JFB-JFA+NPFAC,NPFAC).EQ.NPFAC-1) ISEG1=1
               IF (MOD(JFC-JFB+NPFAC,NPFAC).EQ.1) ISEG2=1
               IF (MOD(JFC-JFB+NPFAC,NPFAC).EQ.NPFAC-1) ISEG2=1
               IF (MOD(JFA-JFC+NPFAC,NPFAC).EQ.1) ISEG3=1
               IF (MOD(JFA-JFC+NPFAC,NPFAC).EQ.NPFAC-1) ISEG3=1
            endif
            IAFA=LFAC(JFA)
            IBFA=LFAC(JFB)
            ICFA=LFAC(JFC)
            IA=IPT1.NUM(IAFA,JEL)
            IB=IPT1.NUM(IBFA,JEL)
            IC=IPT1.NUM(ICFA,JEL)
            XA=XCORD(1,ICPR(IA))
            YA=XCORD(2,ICPR(IA))
            ZA=XCORD(3,ICPR(IA))
            XB=XCORD(1,ICPR(IB))
            YB=XCORD(2,ICPR(IB))
            ZB=XCORD(3,ICPR(IB))
            XC=XCORD(1,ICPR(IC))
            YC=XCORD(2,ICPR(IC))
            ZC=XCORD(3,ICPR(IC))
C  TESTONS SI LES POINTS SONT DU BON COTE
            SIGA=(XA-XCT)*XN+(YA-YCT)*YN+(ZA-ZCT)*ZN
            IF (ABS(SIGA).LE.1D-30) SIGA=1D-30
            SIGB=(XB-XCT)*XN+(YB-YCT)*YN+(ZB-ZCT)*ZN
            IF (ABS(SIGB).LE.1D-30) SIGB=1D-30
            SIGC=(XC-XCT)*XN+(YC-YCT)*YN+(ZC-ZCT)*ZN
            IF (ABS(SIGC).LE.1D-30) SIGC=1D-30
            if (isect.eq.1) goto 350
            IF (SIGA.GT.0..OR.SIGB.GT.0..OR.SIGC.GT.0.) GOTO 103
            IF (VOLUM) THEN
               IA1=MIN(IA,IB,IC)
               IB1=MAX(IA,IB,IC)
               IC1=IA+IB+IC
               DO 102 IDJF=ITRT(/1)-4,1,-5
                  ITR1=MIN(ITRT(IDJF),ITRT(IDJF+1),ITRT(IDJF+2))
                  ITR2=MAX(ITRT(IDJF),ITRT(IDJF+1),ITRT(IDJF+2))
                  ITR3=ITRT(IDJF)+ITRT(IDJF+1)+ITRT(IDJF+2)
                  IF (IA1.NE.ITR1) GOTO 102
                  IF (IB1.NE.ITR2) GOTO 102
                  IF (IC1.NE.ITR3) GOTO 102
*      WRITE (6,*) ' TRIANGLE EN DOUBLE CAS 1 ',IA,IB,IC
C  LE TRIANGLE EST DECLARE NON TRACABLE
                  ITRT(IDJF+4)=MOD(ITRT(IDJF+4),16)+16
                  GOTO 301
 102           CONTINUE
            ENDIF
*     WRITE (6,*) ' NOUVEAU TRIANGLE ',IA,IB,IC
            ITRT(**)=IA
            ITRT(**)=IB
            ITRT(**)=IC
            ITRT(**)=IPT1.ICOLOR(JEL)
            ITRT(**)=ISEG1+2*ISEG2+4*ISEG3
            if (mcham.ne.0) then
               xtrt(**)=melva1.velche(min(lva1,iafa),jel)
               xtrt(**)=melva1.velche(min(lva1,ibfa),jel)
               xtrt(**)=melva1.velche(min(lva1,icfa),jel)
            endif
 301        CONTINUE
            GOTO 101
 103        CONTINUE
      IF (SIGA.GT.0..AND.SIGB.GT.0..AND.SIGC.GT.0.) GOTO 101
      IPOS=0
      IF (SIGA.LE.0.) THEN
        IPOS=IPOS+1
        JTT(IPOS)=IA
        ISEG(IPOS)=ISEG1
      if (mcham.ne.0) xsegm(ia)=
     >           melva1.velche(min(lva1,iafa),jel)
      ENDIF
      IF (SIGA*SIGB.LT.0..AND.(.NOT.(SIGA.LE.0..AND.SIGB.LE.0.))) THEN
        IA1=MIN(IA,IB)
        IB1=MAX(IA,IB)
        DO 104 NSEGM=1,ISEGM(/1),3
        IF (IA1.NE.ISEGM(NSEGM)) GOTO 104
        IF (IB1.NE.ISEGM(NSEGM+1)) GOTO 104
        ISUP=ISEGM(NSEGM+2)
        if (mcham.ne.0)
     >    xsegm(isup)=(melva1.velche(min(lva1,iafa),jel)*sigb-
     >               melva1.velche(min(lva1,ibfa),jel)*siga)/(sigb-siga)
*       WRITE (6,*) ' POINT EN DOUBLE ',IA1,IB1,ISUP
        GOTO 105
 104    CONTINUE
        ITE=ITE+1
        SEGADJ SXCORD
        XCORD(1,ITE)=(XA*SIGB-XB*SIGA)/(SIGB-SIGA)
        XCORD(2,ITE)=(YA*SIGB-YB*SIGA)/(SIGB-SIGA)
        XCORD(3,ITE)=(ZA*SIGB-ZB*SIGA)/(SIGB-SIGA)
        IF (VCHCA.NE.0)
     #   VCHCA(**)=(VCHCA(IA)*SIGB-VCHCA(IB)*SIGA)/(SIGB-SIGA)
        ICPR(**)=ITE
        ISUP=ICPR(/1)
        ISEGM(**)=IA1
        ISEGM(**)=IB1
        ISEGM(**)=ISUP
        if (mcham.ne.0)
     >    xsegm(**)=(melva1.velche(min(lva1,iafa),jel)*sigb-
     >               melva1.velche(min(lva1,ibfa),jel)*siga)/(sigb-siga)
 105    CONTINUE
        IPOS=IPOS+1
        JTT(IPOS)=ISUP
        IF (SIGB.LE.0.) THEN
          ISEG(IPOS)=ISEG1
        ELSE
          ISEG(IPOS)=1
        ENDIF
      ENDIF
      IF (SIGB.LE.0.) THEN
        IPOS=IPOS+1
        JTT(IPOS)=IB
        ISEG(IPOS)=ISEG2
        if (mcham.ne.0) xsegm(ib)=
     >      melva1.velche(min(lva1,ibfa),jel)
      ENDIF
      IF (SIGB*SIGC.LT.0..AND.(.NOT.(SIGB.LE.0..AND.SIGC.LE.0.))) THEN
        IB1=MIN(IB,IC)
        IC1=MAX(IB,IC)
        DO 106 NSEGM=1,ISEGM(/1),3
        IF (IB1.NE.ISEGM(NSEGM)) GOTO 106
        IF (IC1.NE.ISEGM(NSEGM+1)) GOTO 106
        ISUP=ISEGM(NSEGM+2)
        if (mcham.ne.0)
     >    xsegm(isup)=(melva1.velche(min(lva1,ibfa),jel)*sigc-
     >               melva1.velche(min(lva1,icfa),jel)*sigb)/(sigc-sigb)
*       WRITE (6,*) ' POINT EN DOUBLE ',IB1,IC1,ISUP
        GOTO 107
 106    CONTINUE
        ITE=ITE+1
        SEGADJ SXCORD
        XCORD(1,ITE)=(XB*SIGC-XC*SIGB)/(SIGC-SIGB)
        XCORD(2,ITE)=(YB*SIGC-YC*SIGB)/(SIGC-SIGB)
        XCORD(3,ITE)=(ZB*SIGC-ZC*SIGB)/(SIGC-SIGB)
        IF (VCHCA.NE.0)
     #   VCHCA(**)=(VCHCA(IB)*SIGC-VCHCA(IC)*SIGB)/(SIGC-SIGB)
        ICPR(**)=ITE
        ISUP=ICPR(/1)
        ISEGM(**)=IB1
        ISEGM(**)=IC1
        ISEGM(**)=ISUP
        if (mcham.ne.0)
     >    xsegm(**)=(melva1.velche(min(lva1,ibfa),jel)*sigc-
     >               melva1.velche(min(lva1,icfa),jel)*sigb)/(sigc-sigb)
 107    CONTINUE
        IPOS=IPOS+1
        JTT(IPOS)=ISUP
        IF (SIGC.LE.0.) THEN
          ISEG(IPOS)=ISEG2
        ELSE
          ISEG(IPOS)=1
        ENDIF
      ENDIF
      IF (SIGC.LE.0.) THEN
        IPOS=IPOS+1
        JTT(IPOS)=IC
        ISEG(IPOS)=ISEG3
        if (mcham.ne.0) xsegm(ic)=
     >        melva1.velche(min(lva1,icfa),jel)
      ENDIF
      IF (SIGC*SIGA.LT.0..AND.(.NOT.(SIGC.LE.0..AND.SIGA.LE.0.))) THEN
        IC1=MIN(IC,IA)
        IA1=MAX(IC,IA)
        DO 108 NSEGM=1,ISEGM(/1),3
        IF (IC1.NE.ISEGM(NSEGM)) GOTO 108
        IF (IA1.NE.ISEGM(NSEGM+1)) GOTO 108
        ISUP=ISEGM(NSEGM+2)
        if (mcham.ne.0)
     >    xsegm(isup)=(melva1.velche(min(lva1,icfa),jel)*siga-
     >              melva1.velche(min(lva1,iafa),jel)*sigc)/(siga-sigc)
*       WRITE (6,*) ' POINT EN DOUBLE ',IC1,IA1,ISUP
        GOTO 109
 108    CONTINUE
        ITE=ITE+1
        SEGADJ SXCORD
        XCORD(1,ITE)=(XC*SIGA-XA*SIGC)/(SIGA-SIGC)
        XCORD(2,ITE)=(YC*SIGA-YA*SIGC)/(SIGA-SIGC)
        XCORD(3,ITE)=(ZC*SIGA-ZA*SIGC)/(SIGA-SIGC)
        IF (VCHCA.NE.0)
     #   VCHCA(**)=(VCHCA(IC)*SIGA-VCHCA(IA)*SIGC)/(SIGA-SIGC)
        ICPR(**)=ITE
        ISUP=ICPR(/1)
        ISEGM(**)=IC1
        ISEGM(**)=IA1
        ISEGM(**)=ISUP
        if (mcham.ne.0)
     >    xsegm(**)=(melva1.velche(min(lva1,icfa),jel)*siga-
     >              melva1.velche(min(lva1,iafa),jel)*sigc)/(siga-sigc)
 109    CONTINUE
        IPOS=IPOS+1
        JTT(IPOS)=ISUP
        IF (SIGA.LE.0.) THEN
          ISEG(IPOS)=ISEG3
        ELSE
          ISEG(IPOS)=1
        ENDIF
      ENDIF
      IF (IPOS.GT.4) WRITE (6,*) ' ANOMALIE IPOS SURFACE ',IPOS
      IF (IPOS.LT.3) WRITE (6,*) ' ANOMALIE IPOS SURFACE ',IPOS
      IF (IPOS.LT.3) GOTO 101
      IF (IPOS.EQ.3) THEN
      IF (VOLUM) THEN
       IA1=MIN(JTT(1),JTT(2),JTT(3))
       IB1=MAX(JTT(1),JTT(2),JTT(3))
       IC1=JTT(1)+JTT(2)+JTT(3)
       DO 110 IDJF=ITRT(/1)-4,1,-5
       ITR1=MIN(ITRT(IDJF),ITRT(IDJF+1),ITRT(IDJF+2))
       ITR2=MAX(ITRT(IDJF),ITRT(IDJF+1),ITRT(IDJF+2))
       ITR3=ITRT(IDJF)+ITRT(IDJF+1)+ITRT(IDJF+2)
        IF (IA1.NE.ITR1) GOTO 110
        IF (IB1.NE.ITR2) GOTO 110
        IF (IC1.NE.ITR3) GOTO 110
*      WRITE (6,*) ' TRIANGLE EN DOUBLE CAS 2 ',JTT(1),JTT(2),JTT(3)
C  LE TRIANGLE EST DECLARE NON TRACABLE
                     ITRT(IDJF+4)=MOD(ITRT(IDJF+4),16)+16
                     GOTO 303
 110              CONTINUE
               ENDIF
*     WRITE (6,*) ' NOUVEAU TRIANGLE ',JTT(1),JTT(2),JTT(3)
       ITRT(**)=JTT(1)
       ITRT(**)=JTT(2)
       ITRT(**)=JTT(3)
       ITRT(**)=IPT1.ICOLOR(JEL)
       ITRT(**)=ISEG(1)+2*ISEG(2)+4*ISEG(3)
      if (mcham.ne.0) then
       xtrt(**)=xsegm(jtt(1))
       xtrt(**)=xsegm(jtt(2))
       xtrt(**)=xsegm(jtt(3))
      endif
      ENDIF
 303  CONTINUE
      IF (IPOS.EQ.4) THEN
      IF (VOLUM) THEN
      IA1=MIN(JTT(1),JTT(2),JTT(3),JTT(4))
      IB1=MAX(JTT(1),JTT(2),JTT(3),JTT(4))
      IC1=JTT(1)+JTT(2)+JTT(3)+JTT(4)
      ID1=JTT(1)*JTT(2)*JTT(3)*JTT(4)
      DO 111 IDJF=IQUAT(/1)-5,1,-6
      ITR1=MIN(IQUAT(IDJF),IQUAT(IDJF+1),IQUAT(IDJF+2),IQUAT(IDJF+3))
      ITR2=MAX(IQUAT(IDJF),IQUAT(IDJF+1),IQUAT(IDJF+2),IQUAT(IDJF+3))
      ITR3=IQUAT(IDJF)+IQUAT(IDJF+1)+IQUAT(IDJF+2)+IQUAT(IDJF+3)
      ITR4=IQUAT(IDJF)*IQUAT(IDJF+1)*IQUAT(IDJF+2)*IQUAT(IDJF+3)
       IF (IA1.NE.ITR1) GOTO 111
       IF (IB1.NE.ITR2) GOTO 111
       IF (IC1.NE.ITR3) GOTO 111
       IF (ID1.NE.ITR4) GOTO 111
*      WRITE (6,*) ' QUADRILA EN DOUBLE ',JTT(1),JTT(2),JTT(3),JTT(4)
C  LE QUADRILATERE EST DECLARE NON TRACABLE
       IQUAT(IDJF+5)=MOD(IQUAT(IDJF+5),16)+16
       GOTO 302
 111  CONTINUE
      ENDIF
*     WRITE (6,*) ' NOUVEAU QUADRILATERE ',JTT(1),JTT(2),JTT(3),JTT(4)
      IQUAT(**)=JTT(1)
      IQUAT(**)=JTT(2)
      IQUAT(**)=JTT(3)
      IQUAT(**)=JTT(4)
      IQUAT(**)=IPT1.ICOLOR(JEL)
      IQUAT(**)=ISEG(1)+2*ISEG(2)+4*ISEG(3)+8*ISEG(4)
      if (mcham.ne.0) then
       xquat(**)=xsegm(jtt(1))
       xquat(**)=xsegm(jtt(2))
       xquat(**)=xsegm(jtt(3))
       xquat(**)=xsegm(jtt(4))
      endif
      ENDIF
 302  CONTINUE
      GOTO 101
 101  CONTINUE
 350  CONTINUE
      IF (.NOT.VOLUM) GOTO 160
C  ETUDE POUR LE CAS DES VOLUMES DU TETRAEDRE S'APPUYANT SUR LE
C  TRIANGLE ET LE CENTRE DE GRAVITE
C  TETRAEDRE
      IF (SIGA.LT.0..AND.SIGB.LT.0..AND.SIGC.LT.0..AND.SIGRAV
     #  .LT.0.) GOTO 160
      IF (SIGA.GT.0..AND.SIGB.GT.0..AND.SIGC.GT.0..AND.SIGRAV
     #  .GT.0.) GOTO 160
      IPOS=0
      ISEG(1)=0
      ISEG(2)=0
      ISEG(3)=0
      ISEG(4)=0
      IF (SIGRAV*SIGA.LE.0.) THEN
       IF (SIGA.EQ.0.) THEN
        ISUP=IA
        if (mcham.ne.0)
     >    xsegm(isup)=melva1.velche(min(lva1,iafa),jel)
        GOTO 359
       ENDIF
        DO 358 NSEGM=1,NJGRAV
        IF (IA.NE.IJGRAV(1,NSEGM)) GOTO 358
        ISUP=IJGRAV(2,NSEGM)
        if (mcham.ne.0)
     >    xsegm(isup)=(melva1.velche(min(lva1,iafa),jel)*sigrav-
     >              chgrav*siga)/(sigrav-siga)
        GOTO 359
 358    CONTINUE
        ITE=ITE+1
        SEGADJ SXCORD
        XCORD(1,ITE)=(XA*SIGRAV-XGRAV*SIGA)/(SIGRAV-SIGA)
        XCORD(2,ITE)=(YA*SIGRAV-YGRAV*SIGA)/(SIGRAV-SIGA)
        XCORD(3,ITE)=(ZA*SIGRAV-ZGRAV*SIGA)/(SIGRAV-SIGA)
        IF (VCHCA.NE.0)
     #   VCHCA(**)=(VCHCA(IA)*SIGRAV-VGRAV*SIGA)/(SIGRAV-SIGA)
        ICPR(**)=ITE
        ISUP=ICPR(/1)
       if (mcham.ne.0)
     >    xsegm(**)=(melva1.velche(min(lva1,iafa),jel)*sigrav-
     >              chgrav*siga)/(sigrav-siga)
        NJGRAV=NJGRAV+1
        IJGRAV(1,NJGRAV)=IA
        IJGRAV(2,NJGRAV)=ISUP
 359    CONTINUE
        IPOS=IPOS+1
        JTT(IPOS)=ISUP
      ENDIF
      IF (SIGA*SIGB.LT.0..AND.
     # (.NOT.(SIGRAV*SIGA.LE.0..AND.SIGRAV*SIGB.LE.0.))) THEN
        IA1=MIN(IA,IB)
        IB1=MAX(IA,IB)
        DO 404 NSEGM=1,ISEGM(/1),3
        IF (IA1.NE.ISEGM(NSEGM)) GOTO 404
        IF (IB1.NE.ISEGM(NSEGM+1)) GOTO 404
        ISUP=ISEGM(NSEGM+2)
*       WRITE (6,*) ' POINT EN DOUBLE ',IA1,IB1,ISUP
        if (mcham.ne.0)
     >    xsegm(isup)=(melva1.velche(min(lva1,iafa),jel)*sigb-
     >               melva1.velche(min(lva1,ibfa),jel)*siga)/(sigb-siga)
        GOTO 405
 404    CONTINUE
        ITE=ITE+1
        SEGADJ SXCORD
        XCORD(1,ITE)=(XA*SIGB-XB*SIGA)/(SIGB-SIGA)
        XCORD(2,ITE)=(YA*SIGB-YB*SIGA)/(SIGB-SIGA)
        XCORD(3,ITE)=(ZA*SIGB-ZB*SIGA)/(SIGB-SIGA)
        IF (VCHCA.NE.0)
     #   VCHCA(**)=(VCHCA(IA)*SIGB-VCHCA(IB)*SIGA)/(SIGB-SIGA)
        ICPR(**)=ITE
        ISUP=ICPR(/1)
        ISEGM(**)=IA1
        ISEGM(**)=IB1
        ISEGM(**)=ISUP
        if (mcham.ne.0)
     >    xsegm(**)=(melva1.velche(min(lva1,iafa),jel)*sigb-
     >               melva1.velche(min(lva1,ibfa),jel)*siga)/(sigb-siga)
 405    CONTINUE
        IPOS=IPOS+1
        JTT(IPOS)=ISUP
      IF (SIGRAV*SIGB.GT.0.) THEN
        ISEG(IPOS)=1
      ENDIF
      ENDIF
      IF (SIGRAV*SIGB.LE.0.) THEN
       IF (SIGB.EQ.0.) THEN
        ISUP=IB
        if (mcham.ne.0)
     >    xsegm(isup)=melva1.velche(min(lva1,ibfa),jel)
        GOTO 369
       ENDIF
        DO 368 NSEGM=1,NJGRAV
        IF (IB.NE.IJGRAV(1,NSEGM)) GOTO 368
        ISUP=IJGRAV(2,NSEGM)
        GOTO 369
 368    CONTINUE
        ITE=ITE+1
        SEGADJ SXCORD
        XCORD(1,ITE)=(XB*SIGRAV-XGRAV*SIGB)/(SIGRAV-SIGB)
        XCORD(2,ITE)=(YB*SIGRAV-YGRAV*SIGB)/(SIGRAV-SIGB)
        XCORD(3,ITE)=(ZB*SIGRAV-ZGRAV*SIGB)/(SIGRAV-SIGB)
        IF (VCHCA.NE.0)
     #   VCHCA(**)=(VCHCA(IB)*SIGRAV-VGRAV*SIGB)/(SIGRAV-SIGB)
        ICPR(**)=ITE
        ISUP=ICPR(/1)
        if (mcham.ne.0)
     >    xsegm(**)=(melva1.velche(min(lva1,ibfa),jel)*sigrav-
     >               chgrav*sigb)/(sigrav-sigb)
        NJGRAV=NJGRAV+1
        IJGRAV(1,NJGRAV)=IB
        IJGRAV(2,NJGRAV)=ISUP
 369    CONTINUE
        IPOS=IPOS+1
        JTT(IPOS)=ISUP
        if (mcham.ne.0)
     >    xsegm(isup)=(melva1.velche(min(lva1,ibfa),jel)*sigrav-
     >               chgrav*sigb)/(sigrav-sigb)
      ENDIF
      IF (SIGB*SIGC.LT.0..AND.
     # (.NOT.(SIGRAV*SIGB.LE.0..AND.SIGRAV*SIGC.LE.0.))) THEN
        IB1=MIN(IB,IC)
        IC1=MAX(IB,IC)
        DO 406 NSEGM=1,ISEGM(/1),3
        IF (IB1.NE.ISEGM(NSEGM)) GOTO 406
        IF (IC1.NE.ISEGM(NSEGM+1)) GOTO 406
        ISUP=ISEGM(NSEGM+2)
*       WRITE (6,*) ' POINT EN DOUBLE ',IB1,IC1,ISUP
        if (mcham.ne.0)
     >    xsegm(isup)=(melva1.velche(min(lva1,ibfa),jel)*sigc-
     >               melva1.velche(min(lva1,icfa),jel)*sigb)/(sigc-sigb)
        GOTO 407
 406    CONTINUE
        ITE=ITE+1
        SEGADJ SXCORD
        XCORD(1,ITE)=(XB*SIGC-XC*SIGB)/(SIGC-SIGB)
        XCORD(2,ITE)=(YB*SIGC-YC*SIGB)/(SIGC-SIGB)
        XCORD(3,ITE)=(ZB*SIGC-ZC*SIGB)/(SIGC-SIGB)
        IF (VCHCA.NE.0)
     #   VCHCA(**)=(VCHCA(IB)*SIGC-VCHCA(IC)*SIGB)/(SIGC-SIGB)
        ICPR(**)=ITE
        ISUP=ICPR(/1)
        ISEGM(**)=IB1
        ISEGM(**)=IC1
        ISEGM(**)=ISUP
        if (mcham.ne.0)
     >    xsegm(**)=(melva1.velche(min(lva1,ibfa),jel)*sigc-
     >               melva1.velche(min(lva1,icfa),jel)*sigb)/(sigc-sigb)
 407    CONTINUE
        IPOS=IPOS+1
        JTT(IPOS)=ISUP
      IF (SIGRAV*SIGC.GT.0.) THEN
        ISEG(IPOS)=1
      ENDIF
      ENDIF
      IF (SIGRAV*SIGC.LE.0.) THEN
       IF (SIGC.EQ.0.) THEN
        ISUP=IC
       if (mcham.ne.0)
     >    xsegm(isup)=melva1.velche(min(lva1,icfa),jel)
        GOTO 379
       ENDIF
        DO 378 NSEGM=1,NJGRAV
        IF (IC.NE.IJGRAV(1,NSEGM)) GOTO 378
        ISUP=IJGRAV(2,NSEGM)
        if (mcham.ne.0)
     >    xsegm(isup)=(melva1.velche(min(lva1,icfa),jel)*sigrav-
     >               chgrav*sigc)/(sigrav-sigc)
        GOTO 379
 378    CONTINUE
        ITE=ITE+1
        SEGADJ SXCORD
        XCORD(1,ITE)=(XC*SIGRAV-XGRAV*SIGC)/(SIGRAV-SIGC)
        XCORD(2,ITE)=(YC*SIGRAV-YGRAV*SIGC)/(SIGRAV-SIGC)
        XCORD(3,ITE)=(ZC*SIGRAV-ZGRAV*SIGC)/(SIGRAV-SIGC)
        IF (VCHCA.NE.0)
     #   VCHCA(**)=(VCHCA(IC)*SIGRAV-VGRAV*SIGC)/(SIGRAV-SIGC)
        ICPR(**)=ITE
        ISUP=ICPR(/1)
        if (mcham.ne.0)
     >    xsegm(**)=(melva1.velche(min(lva1,icfa),jel)*sigrav-
     >               chgrav*sigc)/(sigrav-sigc)
        NJGRAV=NJGRAV+1
        IJGRAV(1,NJGRAV)=IC
        IJGRAV(2,NJGRAV)=ISUP
 379    CONTINUE
        IPOS=IPOS+1
        JTT(IPOS)=ISUP
      ENDIF
      IF (SIGC*SIGA.LT.0..AND.
     # (.NOT.(SIGRAV*SIGC.LE.0..AND.SIGRAV*SIGA.LE.0.))) THEN
        IC1=MIN(IC,IA)
        IA1=MAX(IC,IA)
        DO 408 NSEGM=1,ISEGM(/1),3
        IF (IC1.NE.ISEGM(NSEGM)) GOTO 408
        IF (IA1.NE.ISEGM(NSEGM+1)) GOTO 408
        ISUP=ISEGM(NSEGM+2)
*       WRITE (6,*) ' POINT EN DOUBLE ',IC1,IA1,ISUP
       if (mcham.ne.0)
     >    xsegm(isup)=(melva1.velche(min(lva1,icfa),jel)*siga-
     >               melva1.velche(min(lva1,iafa),jel)*sigc)/(siga-sigc)
        GOTO 409
 408    CONTINUE
        ITE=ITE+1
        SEGADJ SXCORD
        XCORD(1,ITE)=(XC*SIGA-XA*SIGC)/(SIGA-SIGC)
        XCORD(2,ITE)=(YC*SIGA-YA*SIGC)/(SIGA-SIGC)
        XCORD(3,ITE)=(ZC*SIGA-ZA*SIGC)/(SIGA-SIGC)
        IF (VCHCA.NE.0)
     #   VCHCA(**)=(VCHCA(IC)*SIGA-VCHCA(IA)*SIGC)/(SIGA-SIGC)
        ICPR(**)=ITE
        ISUP=ICPR(/1)
        ISEGM(**)=IC1
        ISEGM(**)=IA1
        ISEGM(**)=ISUP
        if (mcham.ne.0)
     >    xsegm(**)=(melva1.velche(min(lva1,icfa),jel)*siga-
     >               melva1.velche(min(lva1,iafa),jel)*sigc)/(siga-sigc)
 409    CONTINUE
        IPOS=IPOS+1
        JTT(IPOS)=ISUP
      IF (SIGRAV*SIGA.GT.0.) THEN
        ISEG(IPOS)=1
      ENDIF
      ENDIF
      IF (IPOS.GT.4) WRITE (6,*) ' ANOMALIE IPOS VOLUME ',IPOS
      IF (IPOS.LT.3) WRITE (6,*) ' ANOMALIE IPOS VOLUME ',IPOS
      IF (IPOS.LT.3) GOTO 160
      ITR(**)=JTT(1)
      ITR(**)=JTT(2)
      ITR(**)=JTT(3)
      ITR(**)=IPT1.ICOLOR(JEL)
      if (mcham.ne.0) then
       xtr(**)=xsegm(jtt(1))
       xtr(**)=xsegm(jtt(2))
       xtr(**)=xsegm(jtt(3))
      endif
*     WRITE (6,*) ' TRIANGLE (CAS 1) DU PLAN DE COUPE '
*     WRITE (6,*) JTT(1),JTT(2),JTT(3)
C  TRIANGLE DU PLAN DE COUPE
      IF (IPOS.EQ.3) THEN
       MCOUP(**)=ISEG(1)+2*ISEG(2)+4*ISEG(3)+8
      ELSE
       MCOUP(**)=ISEG(1)+2*ISEG(2)+4*0      +8
      ENDIF
 403  CONTINUE
      IF (IPOS.GE.4) THEN
      ITR(**)=JTT(1)
      ITR(**)=JTT(3)
      ITR(**)=JTT(4)
      ITR(**)=IPT1.ICOLOR(JEL)
      if (mcham.ne.0) then
       xtr(**)=xsegm(jtt(1))
       xtr(**)=xsegm(jtt(3))
       xtr(**)=xsegm(jtt(4))
      endif
*     WRITE (6,*) ' TRIANGLE (CAS 2) DU PLAN DE COUPE '
*     WRITE (6,*) JTT(1),JTT(3),JTT(4)
C  TRIANGLE DU PLAN DE COUPE
       MCOUP(**)=0      +2*ISEG(3)+4*ISEG(4)+8
      ENDIF
 402  CONTINUE
 160  CONTINUE
 161  CONTINUE
  20  CONTINUE
      NBELEM=JELN
      SEGADJ IPT2
      ISU=ISU+1
      IPT8.LISOUS(ISU)=IPT2
      if (mcham.ne.0) then
       n1el=nbelem
       segadj melva2
       ipt8.lisref(isu)=melva2
      endif
      SEGDES IPT1
  10  CONTINUE
      IF (NBSOU.NE.0) SEGDES MELEME
  12  CONTINUE
      NBSOUS=ISU+1
*     IF (ITR(/1).EQ.0) NBSOUS=ISU
      NBELEM=0
      NBNN=0
      NBREF=0
      if (mcham.ne.0) nbref=nbsous
      SEGADJ IPT8
      MELEME=IPT8
*     IF (ITR(/1).EQ.0) RETURN
*      WRITE(6,*) ' NOMBRE DE TRIANGLE ENGENDREES DANS LE PLAN DE COUPE',
*     $     ITR(/1)/4
*      WRITE(6,*) ' NOMBRE DE TRIANGLE ENGENDREES HORS LE PLAN DE COUPE',
*     $     ITRT(/1)/5
*      WRITE(6,*) ' NOMBRE DE QUADRILA ENGENDREES HORS LE PLAN DE COUPE',
*     $     IQUAT(/1)/6
      NBELEM=ITR(/1)/4+ITRT(/1)/5+(IQUAT(/1)/6)*2
      NBNN=3
      NBREF=0
      NBSOUS=0
      SEGINI IPT4
      IPT4.ITYPEL=4
      if (mcham.ne.0) then
       n1el=nbelem
       n1ptel=nbnn
       segini melva4
      endif
      ITRL1=ITR(/1)/4
      DO 200 IEL=1,ITRL1
         DO NBNN=1,3
            IPT4.NUM(NBNN,IEL)=ITR((IEL-1)*4+NBNN)
         ENDDO
         if (mcham.ne.0) then
            DO NBNN=1,3
               melva4.velche(nbnn,iel)=xtr((IEL-1)*3+nbnn)
            ENDDO
         ENDIF
         IPT4.ICOLOR(IEL)=ITR(IEL*4)
 200  CONTINUE
      ITRL2=ITRT(/1)/5
      DO 202 IEL=1,ITRL2
         DO 203 NBNN=1,3
            IPT4.NUM(NBNN,IEL+ITRL1)=ITRT((IEL-1)*5+NBNN)
            if (mcham.ne.0)
     >           melva4.velche(nbnn,iel+itrl1)=xtrt((IEL-1)*3+nbnn)
 203     CONTINUE
         IPT4.ICOLOR(IEL+ITRL1)=ITRT(IEL*5-1)
         MCOUP(**)=ITRT(IEL*5)
 202  CONTINUE
      ITRL3=(IQUAT(/1)/6)
      DO 204 IEL=1,ITRL3
         IEL1=2*IEL-1
         DO 205 NBNN=1,3
            IPT4.NUM(NBNN,IEL1+ITRL1+ITRL2)=IQUAT((IEL-1)*6+NBNN)
            if (mcham.ne.0)
     >           melva4.velche(nbnn,iel1+itrl1+itrl2)=xquat((IEL-1)*4
     $           +nbnn)
 205     CONTINUE
         IPT4.ICOLOR(IEL1+ITRL1+ITRL2)=IQUAT(IEL*6-1)
         MCOUP(**)=MOD(IQUAT(IEL*6),4)+16*(IQUAT(IEL*6)/16)
         IEL2=2*IEL
         IPT4.NUM(1,IEL2+ITRL1+ITRL2)=IQUAT((IEL-1)*6+1)
         IPT4.NUM(2,IEL2+ITRL1+ITRL2)=IQUAT((IEL-1)*6+3)
         IPT4.NUM(3,IEL2+ITRL1+ITRL2)=IQUAT((IEL-1)*6+4)
         IPT4.ICOLOR(IEL2+ITRL1+ITRL2)=IQUAT(IEL*6-1)
         MCOUP(**)=0+MOD(IQUAT(IEL*6)/4,2)*2+MOD(IQUAT(IEL*6)/8,2)*4
     $        +16*(IQUAT(IEL*6)/16)
         if (mcham.ne.0) then
            melva4.velche(1,iel2+itrl1+itrl2)=xquat((IEL-1)*4+1)
            melva4.velche(2,iel2+itrl1+itrl2)=xquat((IEL-1)*4+3)
            melva4.velche(3,iel2+itrl1+itrl2)=xquat((IEL-1)*4+4)
         endif
 204  CONTINUE
*     WRITE (6,*) ' LISTE DES ELEMENTS RAJOUTES '
*     DO 9834 J=1,NBELEM
*     WRITE (6,*) (IPT4.NUM(I,J),I=1,3)
*9834 CONTINUE
      LISOUS(ISU+1)=IPT4
      if (mcham.ne.0) lisref(ISU+1)=melva4
      IF (MELEM2.NE.0) THEN
         IF (MEDSAU.EQ.0) THEN
C           ON COUPE MELEM2
            mcham=0
            MEDSAU=MELEME
            MEDCOU=MCOUP
            MELEME=MELEM2
            SEGSUP ITR,ITRT,IQUAT
            SEGINI MCOUP,ITR,ITRT,IQUAT
            GOTO 5000
         ELSE
            MELEM2=MELEME
            MELEME=MEDSAU
            MCOU2=MCOUP
            MCOUP=MEDCOU
         ENDIF
      ENDIF
      SEGSUP ITR,ISEGM,IJGRAV,ITRT,IQUAT
      if (mchamor.ne.0) segsup xtr,xtrt,xquat,xsegm
*dbg      write(ioimp,*) 'sortie crcou2' 
      RETURN
      END