C FACED     SOURCE    PV        18/06/15    21:15:01     9859           
C TRACE D'UN OBJET PAR FACE EN COMMENCANT PAR CELLES DE DERRIERE
C
C  imod=0 trace en couleur d'effacement
C  imod=1 trace en couleir normale
C  SG 2016/07/18 Programmation comme envvo2 pour gerer les faces TRI7/QUA9
C
      SUBROUTINE FACED(MELEME,XPROJ,ICPR,IVU,MCOUP,KON,LNDEGR,imod)
      IMPLICIT INTEGER(I-N)

-INC PPARAM
-INC CCOPTIO
-INC SMELEME
-INC CCGEOME
-INC CCREEL
C+PP DEGRADE OU NON
      LOGICAL lndegr
C+PP
      SEGMENT XPROJ(3,1)
      SEGMENT ICPR(1)
      SEGMENT IVU(ITE)
      SEGMENT MCOUP(0)
      SEGMENT KON(3,NBCON,NMAX)
*
*     Type de faces prises en compte: T3, Q4, T6, Q8, POLY, T7, Q9
*     Numero dans KDFAC               1   2   3   4     6    7   8
* Pour ne pas se prendre la tête,  on numerote pareil que KDFAC
* Pour les 5 (non utilisé), 6 (polygone) et >8, ca restera à 0
*  NTYFAC=20, codé en dur dans CCGEOME pour KDFAC
      PARAMETER (NTYFAC=20)
* Nb de faces de chaque type, sert également de compteur
      SEGMENT NBFAC(NTYFAC)
* Un segment pointant sur les IFACI
      SEGMENT IPOFAC(NTYFAC)
*     Segment IFACI contenant les noeuds, la couleur et si la face d'un
*     type donné est vue ou non
      SEGMENT IFACI(NNODE+2,NFACI)
*     Nombre de noeuds max pour les polygones
      PARAMETER (NNOMAX=14)
*
      SEGMENT NSOMP(NFACP)
      SEGMENT IFACOL(NFAC)
      SEGMENT TFAC(NFAC)
      SEGMENT KFAK(NFAC)
      SEGMENT NAUX(max(2,NFAC))
      DIMENSION XTR(NNOMAX),YTR(NNOMAX),ZTR(NNOMAX),NPTR(NNOMAX)
      DIMENSION XTR2(3),YTR2(3),ZTR2(3)
*
*dbg      write(ioimp,*) 'coucou faced lndegr=',LNDEGR,' imod=',imod

      SEGACT MELEME
      MELSAU=MELEME
      IPT1=MELEME
*     IF (MCOUP.NE.0) THEN
*      NBNN=0
*      NBELEM=0
*      NBREF=0
*      NBSOUS=LISOUS(/1)-1
*      SEGINI IPT1
*      DO 5 I=1,NBSOUS
*      IPT1.LISOUS(I)=LISOUS(I)
*   5  CONTINUE
*      ISCOUP=LISOUS(NBSOUS+1)
*     ENDIF
*     call ecmail(ipt1,0)
      CALL ECROBJ('MAILLAGE',IPT1)
      CALL ENVELO
*      IF (MCOUP.NE.0) SEGSUP IPT1
      CALL LIROBJ('MAILLAGE',MELEME,1,IRETOU)
      IF (IERR.NE.0) RETURN
      SEGACT MELEME
*      IF (MCOUP.NE.0) THEN
*        IF (LISOUS(/1).NE.0) THEN
*         NBSOUS=LISOUS(/1)+1
*         SEGADJ MELEME
*         LISOUS(NBSOUS)=ISCOUP
*        ELSE
*         NBSOUS=2
*         SEGINI IPT1
*         IPT1.LISOUS(1)=MELEME
*         IPT1.LISOUS(2)=ISCOUP
*         MELEME=IPT1
*        ENDIF
*      ENDIF
      SEGACT XPROJ,ICPR
      SEGACT KON*MOD
      NBCON=KON(/2)
      NBCONR=NBCON-1
      NBPOIN=XPROJ(/2)
      TMIN=xsgran
      TMAX=-xsgran
      DO 1 I=1,NBPOIN
         TMIN=MIN(TMIN,XPROJ(3,I))
         TMAX=MAX(TMAX,XPROJ(3,I))
 1    CONTINUE
      TDIST=TMAX-TMIN
c
c     on compte le nombre d elements dont les faces sont de type 1 2 3 4
c     6 7 8 dans NBFAC, attention à 6 : gestion des polygones
      SEGINI NBFAC
      IPT1=MELEME
      SEGACT MELEME
      DO 10 IOB=1,MAX(1,LISOUS(/1))
         IF (LISOUS(/1).NE.0) THEN
            IPT1=LISOUS(IOB)
            SEGACT IPT1
         ENDIF
         NBELEM=IPT1.NUM(/2)
         ILTEL=LTEL(1,IPT1.ITYPEL)
         IF (ILTEL.EQ.0) GOTO 12
         ILTAD=LTEL(2,IPT1.ITYPEL)
         DO 13 IF=1,ILTEL
            IFT=LDEL(1,ILTAD+IF-1)
            NBFAC(IFT)=NBFAC(IFT)+NBELEM
 13      CONTINUE
 12      CONTINUE
 10   CONTINUE
*
      NFACP=NBFAC(6)
      SEGINI NSOMP
c==== CREATION DES FACES ==============================================
* Initialisation des IFACI
      SEGINI IPOFAC
      DO ITYFAC=1,NTYFAC
         NNODE=KDFAC(1,ITYFAC)
* Polygone
         IF (ITYFAC.EQ.6) NNODE=NNOMAX
         IF (NNODE.GT.0) THEN
            NFACI=NBFAC(ITYFAC)
            SEGINI IFACI
            IPOFAC(ITYFAC)=IFACI
         ENDIF
      ENDDO
c     NBFAC sert maintenant de compteur
      DO ITYFAC=1,NTYFAC
         NBFAC(ITYFAC)=0
      ENDDO
      ICOUPE=0
      DO 50 IOB=1,MAX(1,LISOUS(/1))
         IF (LISOUS(/1).NE.0) THEN
            IPT1=LISOUS(IOB)
            ICOUPE=0
            IF (IOB.EQ.LISOUS(/1).AND.MCOUP.NE.0) ICOUPE=1
         ENDIF
         NBELEM=IPT1.NUM(/2)
         ILTEL=LTEL(1,IPT1.ITYPEL)
         IF (ILTEL.EQ.0) GOTO 52
         ILTAD=LTEL(2,IPT1.ITYPEL)
         DO 60 IF=1,ILTEL
            ITYFAC=LDEL(1,ILTAD+IF-1)
            IAD=LDEL(2,ILTAD+IF-1)
            NNODE=KDFAC(1,ITYFAC)
            NNODF=NNODE
* Polygone
            IF (ITYFAC.EQ.6) THEN
               NBNN=IPT1.NUM(/1)
*  23 1
*Erreur dans le module de trace
               IF (NBNN.GT.NNOMAX) THEN
                  CALL ERREUR(23)
                  RETURN
               ENDIF
               NNODE=NBNN
               NNODF=NNOMAX
            ENDIF
            IF (NNODE.GT.0) THEN 
               IFACI=IPOFAC(ITYFAC)
               DO 70 IEL=1,NBELEM
                  NBFAC(ITYFAC)=NBFAC(ITYFAC)+1
                  j=NBFAC(ITYFAC)
* Polygone
                  IF (ITYFAC.EQ.6) NSOMP(j)=NNODE
                  IFACI(NNODF+2,j)=1
                  DO i=1,NNODE
                     IFACI(i,j)=ICPR(IPT1.NUM(LFAC(IAD+i-1),IEL))
                     IF (IVU(IFACI(i,j)).NE.1) IFACI(NNODF+2,j)=0
                  ENDDO
                  IFACI(NNODF+1,j)=IPT1.ICOLOR(IEL)
* TRI3 cas des coupes
                  IF (ITYFAC.EQ.1) THEN
                     IF (ICOUPE.EQ.1) THEN
                        IF (MCOUP(IEL)/8.EQ.1) IFACI(NNODF+2,j)=2
                        IF (MCOUP(IEL)/16.EQ.1) IFACI(NNODF+2,j)=3
                     ENDIF
                  ENDIF
 70            CONTINUE
            ENDIF
 60      CONTINUE
 52      CONTINUE
         IF (LISOUS(/1).NE.0) SEGDES IPT1
 50   CONTINUE
      SEGDES MELEME
C  IF FAUT MAINTENANT RETASSER ET CLASSER LES TABLEAUX DES FACES
C  PROBLEME ELLES NE SONT PAS FORCEMENT DECRITES DE LA MEME FACON
      NFAC=0
      DO ITYFAC=1,NTYFAC
         NFAC=NFAC+NBFAC(ITYFAC)
      ENDDO
      SEGINI TFAC,KFAK,IFACOL
      IFAC=0
      DO ITYFAC=1,NTYFAC
         NNODE=KDFAC(1,ITYFAC)
         IF (ITYFAC.EQ.6) THEN
            NNODF=NNOMAX
         ELSE
            NNODF=NNODE
         ENDIF
         IF (NNODE.GT.0.OR.ITYFAC.EQ.6) THEN
            IFACI=IPOFAC(ITYFAC)
            DO I=1,NBFAC(ITYFAC)
               IFAC=IFAC+1
* Polygone
               IF (ITYFAC.EQ.6) NNODE=NSOMP(I)
               XXXX = 0.
               DO J=1,NNODE
* ligne suivante : erreur de faced pour poly ?
*                  XXXX = XXXX + XPROJ(3,ICPR(IFACI(J,I)))
                  XXXX = XXXX + XPROJ(3,IFACI(J,I))
               ENDDO
               XXXX=XXXX/NNODE
               TFAC(IFAC)=XXXX
               IF (IFACI(NNODF+2,I).EQ.1) TFAC(IFAC)=TFAC(IFAC)-TDIST
               IFACOL(IFAC)=IFACI(NNODF+1,I)
               KFAK(IFAC)=I+((ITYFAC-1)*NFAC)
* TRI3/coupe
               IF (ITYFAC.EQ.1) THEN
                  IF (IFACI(NNODF+2,I).EQ.2) TFAC(I)=TFAC(IFAC)-2*TDIST
                  IF (IFACI(NNODF+2,I).EQ.3) KFAK(IFAC)=0
               ENDIF
            ENDDO
         ENDIF
      ENDDO
C  IL N'Y A PLUS QU'A TRIER ET RETASSER KFAK SUIVANT TFAC
      SEGINI NAUX
      IF (IDIM.EQ.2) GOTO 209
      NAUX(1)=1
      NAUX(2)=NFAC
      IZ=2
 208  IZ=IZ-1
      IF (IZ.LE.0) GOTO 209
      IPB=NAUX(IZ*2-1)
      IPH=NAUX(IZ*2)
      IF(IPB.GE.IPH) GOTO 208
      JPB=IPB-1
      JPH=IPH+1
C  CALCUL DU PIVOT
      PV=0.
*     DO 207 J=IPB,IPH
*     PV=PV+TFAC(J)
*207  CONTINUE
*     PV=PV/(IPH-IPB+1)
      PV=(TFAC(IPB)+TFAC(IPH))/2.
 242  JPB=JPB+1
      IF (JPH.EQ.JPB) GOTO 245
      IF (TFAC(JPB).LT.PV) GOTO 243
      GOTO 242
 243  JPH=JPH-1
      IF (JPH.EQ.JPB) GOTO 245
      IF (TFAC(JPH).GT.PV) GOTO 244
      GOTO 243
 244  IAUX=KFAK(JPB)
      IAUXX=IFACOL(JPB)
      KFAK(JPB)=KFAK(JPH)
      IFACOL(JPB)=IFACOL(JPH)
      KFAK(JPH)=IAUX
      IFACOL(JPH)=IAUXX
      TAUX=TFAC(JPB)
      TFAC(JPB)=TFAC(JPH)
      TFAC(JPH)=TAUX
      GOTO 242
 245  IF (JPB.GE.IPB) GOTO 247
      JPB=JPB+1
      JPH=JPH+2
      GOTO 248
 247  IF (JPH.LE.IPH) GOTO 249
      JPB=JPB-2
      JPH=JPH-1
      GOTO 248
 249  IF (TFAC(JPB).LE.PV) GOTO 250
      IF (JPH.EQ.IPH) GOTO 251
 252  JPH=JPH+1
      GOTO 248
 250  IF (JPB.EQ.IPB) GOTO 252
 251  JPB=JPB-1
 248  IF (JPB.EQ.IPB) GOTO 253
      NAUX(2*IZ)=JPB
      IZ=IZ+1
 253  IF (JPH.EQ.IPH) GOTO 208
      NAUX(2*IZ)=IPH
      NAUX(2*IZ-1)=JPH
      IZ=IZ+1
      GOTO 208
 209  CONTINUE
C  LES FACES SONT CLASSEES DANS KFAK LES FACES ON ETE ELIMINEES PAR
C  ENVELO .   IL NE RESTE PLUS QU'A TRACER
      DO 300 I=1,NFAC
         IFF=KFAK(I)
         IF (IFF.EQ.0) GOTO 300
         IT=(IFF-1)/NFAC
         IF=IFF-IT*NFAC
         IT=IT+1
         IOK=0
*     
         ITYFAC=IT
         NNODE=KDFAC(1,ITYFAC)
         IF (ITYFAC.EQ.6) NNODE=NSOMP(IF)
         IF (NNODE.GT.0) THEN
            IFACI=IPOFAC(ITYFAC)
            DO IP=1,NNODE
               XTR(IP)=XPROJ(1,IFACI(IP,IF))
               YTR(IP)=XPROJ(2,IFACI(IP,IF))
               ZTR(IP)=XPROJ(3,IFACI(IP,IF))
               NPTR(IP)=IFACI(IP,IF)
               IF (IVU(IFACI(IP,IF)).EQ.1) IOK=1
*!!               write(ioimp,*) 'ip=',ip,' x y z n iok',xtr(ip),ytr(ip),
*!!     $              ztr(ip), iok
            ENDDO
         ENDIF
         IF (IOK.EQ.0) GOTO 300
         NP=NNODE
* SG Pour les TRI7,QUA9 seuls les points du contour sont pris en compte
*    pour le calcul de la normale et la suppression des points du contour.
         IF (ITYFAC.EQ.7.OR.ITYFAC.EQ.8) NP=NP-1
*!!         write(ioimp,*) 'np=',np
C  DETERMINATION DE LA COULEUR D'APRES L'ORIENTATION DE LA NORMALE
         XN=0.
         YN=0.
         ZN=0.
         XDE=XTR(NP)
         YDE=YTR(NP)
         ZDE=ZTR(NP)
         XW1=XTR(1)-XDE
         YW1=YTR(1)-YDE
         ZW1=ZTR(1)-ZDE
         DO 900 I2=2,NP-1
            XW2=XTR(I2)-XDE
            YW2=YTR(I2)-YDE
            ZW2=ZTR(I2)-ZDE
            XN=XN+(YW1*ZW2-ZW1*YW2)
            YN=YN+(ZW1*XW2-XW1*ZW2)
            ZN=ZN+(XW1*YW2-YW1*XW2)
            XW1=XW2
            YW1=YW2
            ZW1=ZW2
 900     CONTINUE
         DN=SQRT(XN**2+YN**2+ZN**2)
         IF (DN.EQ.0.) DN=1.
         ZN=ASIN(ABS(ZN/DN))
C+PP DEGRADE OU NON
         IF (lndegr) ZN=REAL(XPI/2.D0)
C+PP
         IFACO=IFACOL(I)
         if (imod.eq.0) ifaco=8
* SG Pour les TRI7 et QUA9, on choisit de decouper en TRI3 :
*    idealement, on aurait mis le calcul de la normale dedans
         IF (ITYFAC.EQ.7.OR.ITYFAC.EQ.8) THEN
            XTR2(1)=XTR(NP+1)
            YTR2(1)=YTR(NP+1)
            ZTR2(1)=ZTR(NP+1)
            DO IPO=1,NP
               IPOP=IPO+1
               IF (IPOP.GT.NP) IPOP=1
               XTR2(2)=XTR(IPO)
               YTR2(2)=YTR(IPO)
               ZTR2(2)=ZTR(IPO)
               XTR2(3)=XTR(IPOP)
               YTR2(3)=YTR(IPOP)
               ZTR2(3)=ZTR(IPOP)
               CALL TRFACE(3,XTR2,YTR2,ZTR2,ZN,IFACO,IEFF)
               if (imod.eq.0) ieff=0
            ENDDO
         ELSE
            CALL TRFACE(NP,XTR,YTR,ZTR,ZN,IFACO,IEFF)
         ENDIF
         

         if (imod.eq.0) ieff=0
C  SI IEFF <> 0 IL FAUT METTRE LES TRAITS EN EFFACEMENT
         N1=NPTR(NP)
         DO 450 IIP=1,NP
            N2=NPTR(IIP)
            NI=N1
            NJ=N2
* 459        CONTINUE
 457        DO 454 K=1,NBCONR
               IF (KON(1,K,NI).NE.NJ) GOTO 454
C  8 = EFFACEMENT
               KON(2,K,NI)=8
               IF (IEFF.EQ.0) KON(2,K,NI)=IFACOL(I)
               GOTO 456
 454        CONTINUE
            IF (KON(1,NBCON,NI).EQ.0) GOTO 456
            NI=KON(1,NBCON,NI)
            GOTO 457
 456        CONTINUE
            NI=N2
            NJ=N1
* 469        CONTINUE
 467        DO 464 K=1,NBCONR
               IF (KON(1,K,NI).NE.NJ) GOTO 464
C     8 = EFFACEMENT
               KON(2,K,NI)=8
               IF (IEFF.EQ.0) KON(2,K,NI)=IFACOL(I)
               GOTO 466
 464        CONTINUE
            IF (KON(1,NBCON,NI).EQ.0) GOTO 466
            NI=KON(1,NBCON,NI)
            GOTO 467
 466        CONTINUE
            N1=N2
 450     CONTINUE
 300  CONTINUE
C'EST FINI
      SEGACT MELEME
      IF (LISOUS(/1).NE.0) THEN
         NBSOUS=LISOUS(/1)
*     IF (MCOUP.NE.0) NBSOUS=NBSOUS-1
         DO 490 IO=1,NBSOUS
            IPT2=LISOUS(IO)
            segact ipt2
            if (ipt2.itypel.gt.3.AND.ipt2.itypel.NE.32) SEGSUP IPT2
 490     CONTINUE
      ENDIF
      if ((itypel.eq.0) .AND. (MELEME .NE. MELSAU)) SEGSUP MELEME
      MELEME=MELSAU
      SEGDES KON
      DO ITYFAC=1,NTYFAC
         IFACI=IPOFAC(ITYFAC)
         IF (IFACI.NE.0) THEN
            SEGSUP IFACI
         ENDIF
      ENDDO
      SEGSUP IPOFAC
      SEGSUP NBFAC
      SEGSUP NAUX,TFAC,KFAK,IFACOL,NSOMP
      RETURN
      END