C TETRA     SOURCE    JC220346  16/11/29    21:15:38     9221           
C---------------------------------------------------------------------|
C                                                                     |
       SUBROUTINE TETRA(II,JJ,IF1,IF2,IGAGNE,IBOUT)
C                                                                     |
C      CETTE SUBROUTINE TENTE DE CREER UN TETRAEDRE A PARTIR          |
C      DES 2 TRIANGLES IF1 ET IF2.                                    |
C                                                                     |
C      -1- RECHERCHE DES FACETTES IF3 ET IF4                          |
C      -2- CREATION DES FACETTES INEXISTENTES                         |
C      -3- MISE A JOUR DU MAILLAGE DE SURFACE                         |
C      -4- TEST DU VOLUME ELEMENTAIRE CREE                            |
C      -5- MEMORISATION DU VOLUME CREE                                |
C                                                                     |
C      - IGAGNE=1 EN CAS DE SUCCES                                    |
C      - IGAGNE=0 EN CAS D'ECHEC                                      |
C                                                                     |
C---------------------------------------------------------------------|
C
      IMPLICIT INTEGER(I-N)
      IMPLICIT REAL*8(A-H,O-Z)
-INC TDEMAIT

-INC PPARAM
-INC CCOPTIO
       LOGICAL REPONS,FACET,SOLTET,DIAGO,IN2,VERDIV,IN
       dimension ipts(10),jpts(10),qualti(10),qualtj(10)
C
       inouv=0
       idepl=0
       if (nfc(4,if1).ne.0) return
       if (nfc(4,if2).ne.0) return
C
*   on se refuse a mailler un tetraedre trop ouvert
*       if (TETA(IF1,IF2,II,JJ).lt.-100) then
*        write (6,*) 'if1 if2 ii jj TETA ',IF1,IF2,TETA(IF1,IF2,II,JJ)
*        goto 155
*       endif
       nfcini=nfcmax
       nptini=nptmax
       ip1=0
       ip2=0
       N3=0
       N4=0
       ICTF=0
       ICTV=0
       ipas=0
       IP=IPRED(IF1,II)
       JP=ISUCC(IF2,II)
*
       if (TETA(IF1,IF2,II,JJ).lt.-100D0) goto 250
C
C      -1- RECHERCHE DES FACETTES IF3 ET IF4
C      -------------------------------------
C      RECHERCHE D'UNE 3EME FACETTE IF3
C      --------------------------------
C
       IF3=IFACE3(IP,II,JP)
         IF (IF3.NE.0) THEN
*  si if3 dans le mauvais sens rien a faire
       if (isucc(if3,ii).ne.ip) then
        IF (IVERB.EQ.1) write (6,*) ' tetra face a l''envers if3 '
        return
       endif
*  SI IF3 QUADRANGULAIRE RENVOYER SUR PYRA1
           IF (NFC(4,IF3).NE.0) THEN
              CALL PYRA1(IP,II,IF1,IF3,IGAGNE,IBOUT)
              RETURN
           ENDIF
           N3=IF3
*         WRITE(6,1010)IF3
1010      FORMAT(' ** IF3=',I3)
         ENDIF
C
C
C
C      RECHERCHE D'UNE 4EME FACETTE IF4
C      --------------------------------
C
       IF4=IFACE3(IP,JJ,JP)
         IF (IF4.NE.0) THEN
*  si if4 dans le mauvais sens rien a faire
       if (isucc(if4,jj).ne.jp) then
        IF (IVERB.EQ.1) write (6,*) ' tetra face a l''envers if4 '
        return
       endif
           IF (NFC(4,IF4).NE.0) THEN
*  SI IF4 QUADRANGULAIRE RENVOYER SUR PYRA1
              CALL PYRA1(JJ,IP,IF1,IF4,IGAGNE,IBOUT)
              RETURN
           ENDIF
           N4=IF4
*          WRITE(6,1020)IF4
1020      FORMAT(' ** IF4=',I3)
         ENDIF
C
C      -2- CREATION DES FACETTES INEXISTANTES
C      --------------------------------------
       IF (N3.EQ.0) THEN
C
C      CREATION DE LA FACETTE IF3
C      --------------------------
       NFCMAX=NFCMAX+1
       IF3=NFCMAX
       ICTF=ICTF+1
C
       NFC(1,IF3)=II
       NFC(2,IF3)=JP
       NFC(3,IF3)=IP
       NFC(4,IF3)=0
C
       ENDIF
C
       IF (N4.EQ.0) THEN
C
C      CREATION DE LA FACETTE IF4
C      --------------------------
       NFCMAX=NFCMAX+1
       IF4=NFCMAX
       ICTF=ICTF+1
C
       NFC(1,IF4)=JJ
       NFC(2,IF4)=IP
       NFC(3,IF4)=JP
       NFC(4,IF4)=0
C
       ENDIF
C
C  VERIFICATION VALIDITE ARETE CREE
       IF (n3.EQ.0.AND.n4.EQ.0) THEN
       IF (DIAGO(IP,JP,diacri)) then
        iff=ifdiag
        lpts=0
        nfcmax=nfcini
        if (iff.ne.0) then
        if (nfc(1,iff).eq.ip.or.nfc(2,iff).eq.ip.or.nfc(3,iff).eq.ip)
     >     then
            ipt1=isucc(iff,ip)
            if (ipt1.ne.ii.and.ipt1.ne.jj) then
              lpts=lpts+1
              ipts(lpts)=ipt1
            endif
            ipt1=ipred(iff,ip)
            if (ipt1.ne.ii.and.ipt1.ne.jj) then
              lpts=lpts+1
              ipts(lpts)=ipt1
            endif
         endif
        if (nfc(1,iff).eq.jp.or.nfc(2,iff).eq.jp.or.nfc(3,iff).eq.jp)
     >     then
            ipt1=isucc(iff,jp)
            if (ipt1.ne.ii.and.ipt1.ne.jj) then
              lpts=lpts+1
              ipts(lpts)=ipt1
            endif
            ipt1=ipred(iff,jp)
            if (ipt1.ne.ii.and.ipt1.ne.jj) then
              lpts=lpts+1
              ipts(lpts)=ipt1
            endif
         endif
         lpts1=0
         do 700 lpt=1,lpts
          do 701 lpt1=1,lpts1
            if (ipts(lpt).eq.jpts(lpt1)) goto 700
 701      continue
          qualti(lpt)=min(qualt(ii,jj,ip,ipts(lpt)),
     >                    qualt(jj,ii,jp,ipts(lpt)))
          if (qualti(lpt).le.0.001d0) goto 700
          lpts1=lpts1+1
          jpts(lpts1)=ipts(lpt)
          qualtj(lpts1)=qualti(lpt)
 700     continue
         lpts=lpts1
         lpts1=0
         do 702 lpt=1,lpts
          do 703 l=1,lpts1
           if (qualtj(lpt).gt.qualti(l)) goto 704
 703      continue
          l=lpts1+1
 704      continue
          lpts1=lpts1+1
          do 705 l1=lpts1,l+1,-1
            qualti(l1)=qualti(l1-1)
            ipts(l1)=ipts(l1-1)
 705      continue
            qualti(l)=qualtj(lpt)
            ipts(l)=jpts(lpt)
 702     continue
         do 706 lpt=1,lpts
*          write (6,*) 'tetra appel tetra2 type 1 avec ',
*    >     ipts(lpt),iff,qualti(lpt)
            call tetra2(ii,jj,if1,if2,igagne,ipts(lpt))
            if (igagne.eq.1) return
 706     continue
         if (npf(5,ii).ne.0.or.npf(5,jj).ne.0) goto 250
*        write (6,*) ' tetra diago 4-4 => point milieu',ip,jp
         goto 350
        endif
*        write (6,*) ' tetra echec 1 lpts',lpts
         return
         ENDIF
*  test in2
*       IF (.NOT.IN2(ip,jp,nfcini)) THEN
*         write (6,*) 'tetra test in2 echoue avec ',ip,jp
*         nfcmax=nfcini
*         nptmax=nptini
*         return
*       ENDIF
       ENDIF
C
C      -3- MISE A JOUR DU MAILLAGE DE SURFACE
C      --------------------------------------
       CALL REPSUB(IF1)
       CALL REPSUB(IF2)
       CALL REPSUB(IF3)
       CALL REPSUB(IF4)
C
C      -4- TEST DU VOLUME ELEMENTAIRE CREE
C      -----------------------------------
       IF (N3.EQ.0) THEN
         IF (.NOT.FACET(IF3))     then
*          write(6,*) ' tetra facet if3 invalide'
           GOTO 150
         endif
       ENDIF
       IF (N4.EQ.0) THEN
         IF (.NOT.FACET(IF4))    then
*          write(6,*) ' tetra facet if4 invalide'
           GOTO 150
         endif
       ENDIF
*
*  verification longueur de l'arete a creer
*
         DNORM=(XYZ(1,IP)-XYZ(1,JP))**2
     #        +(XYZ(2,IP)-XYZ(2,JP))**2
     #        +(XYZ(3,IP)-XYZ(3,JP))**2
         DTEST=tetrl*XYZ(4,IP)*XYZ(4,JP)
*        write (6,*) ' tetra dnorm dtest ',dnorm,dtest
         IF (DNORM.GT.DTEST) THEN
           CALL REPSUB(IF4)
           CALL REPSUB(IF3)
           CALL REPSUB(IF2)
           CALL REPSUB(IF1)
           NFCMAX=NFCini
           GOTO 250
         ENDIF
*  verification validite du tetraedre
       IF (.NOT.SOLTET(IF1,IF2,IF3,IF4,IPTT)) GOTO 151
C
C      LE VOLUME CREE EST VALIDE.
C      --------------------------
C      -5- MEMORISATION DU VOLUME CREE
C      -------------------------------
*       write (6,*) ' tetra creation tetraedre '
       NVOL=NVOL+1
       IF (NFV(1,IF1).EQ.0) NFV(1,IF1)=NVOL
       IF (NFV(1,IF1).NE.NVOL) NFV(2,IF1)=NVOL
       IF (NFV(1,IF2).EQ.0) NFV(1,IF2)=NVOL
       IF (NFV(1,IF2).NE.NVOL) NFV(2,IF2)=NVOL
       IF (NFV(1,IF3).EQ.0) NFV(1,IF3)=NVOL
       IF (NFV(1,IF3).NE.NVOL) NFV(2,IF3)=NVOL
       IF (NFV(1,IF4).EQ.0) NFV(1,IF4)=NVOL
       IF (NFV(1,IF4).NE.NVOL) NFV(2,IF4)=NVOL
       IVOL(9,NVOL)=25
       DO 130 I=1,3
          IVOL(I,NVOL)=NFC(I,IF1)
130    CONTINUE
       IVOL(4,NVOL)=JP
C
*      WRITE(6,1070)NVOL,(IVOL(I,NVOL),I=1,9)
*1070   FORMAT(I3,4X,14I4)
       do 961 npfi=40,1,-1
        if (npf(npfi,ii).ne.0) goto 962
 961   continue
 962   continue
       do 963 npfj=40,1,-1
        if (npf(npfj,jj).ne.0) goto 964
 963   continue
 964   continue
       qual=qualt(ii,jj,ip,jp)
       if (iimpi.eq.1) write (6,1100) nvol,nfcmax,nfacet,qual,
     >   (ivol(i,nvol),i=1,4)
 1100  format (' TETRA  ',i5,2i6,f8.4,4i6)
C
*      DO 140 J=1,NPTMAX
*          WRITE(6,1080)J,(NPF(I,J),I=1,40)
*1080      FORMAT(I4,4X,40I3)
* 140    CONTINUE
C
       IGAGNE=1
C
       RETURN
C
350    CONTINUE
*  IP JP est invalide (diago)
*  IP JP II JJ ont 4 facettes
*  on va essayer de mettre un point au centre de
*  gravite de ii jj ip jp et des 2 noeuds supplementaires
*  des facettes touchant IP JP II JJ
       IF3=NOISIN(IP,II,IF1)
       if (nfc(4,if3).ne.0) return
       IP1=ISUCC(IF3,IP)
       IF5=NOISIN(II,JP,IF2)
       if (nfc(4,if5).ne.0) return
       IP3=ISUCC(IF5,II)
       IF4=NOISIN(JP,JJ,IF2)
       if (nfc(4,if4).ne.0) return
       IP2=ISUCC(IF4,JP)
       IF6=NOISIN(JJ,IP,IF1)
       if (nfc(4,if6).ne.0) return
       IP4=ISUCC(IF6,JJ)
*      write (6,*) ' II JJ IP JP IP1 IP2 IP3 IP4',
*    >   II,JJ,IP,JP,IP1,IP2,IP3,IP4
       NFCMAX=NFCini
       idepl=3
       inouv=1
       DO 351 I=1,4
          XYZ(I,NPTMAX+1)=(XYZ(I,II)+XYZ(I,JJ)+XYZ(I,IP)+
     >                   XYZ(I,JP)+(XYZ(I,IP1)+XYZ(I,IP3))/2+
     >                   (XYZ(I,IP2)+XYZ(I,IP4))/2)/6
 351   CONTINUE
*      call vcrit(nptmax+1)
 356   CONTINUE
*  verif de validite de la position du point
       iechec=0
       v1=0.
       v2=0.
       v3=0.
       v4=0.
       v5=0.
       v6=0.
       v7=0.
       v8=0.
       v9=0.
       v10=0.
       v11=0.
       v12=0.
       v1=vol(nptmax+1,ii,ip,ip1)
       if (v1.gt.0.) then
       iechec=1
*       write (6,*) ' tetra vol 1 positif ',v1
       endif
       v2=vol(nptmax+1,ii,ip1,jp)
       if (v2.gt.0.) then
       iechec=1
*       write (6,*) ' tetra vol 2 positif ',v2
       endif
       if (ip1.ne.ip3) then
        v3=vol(nptmax+1,ii,ip,ip3)
        if (v3.gt.0.) then
       iechec=1
*        write (6,*) ' tetra vol 3 positif ',v3
        endif
        v4=vol(nptmax+1,ii,ip3,jp)
        if (v4.gt.0.) then
       iechec=1
*        write (6,*) ' tetra vol 4 positif ',v4
        endif
       endif
       v5=vol(nptmax+1,ip2,jj,jp)
       if (v5.gt.0.) then
       iechec=1
*       write (6,*) ' tetra vol 5 positif ',v5
       endif
       v6=vol(nptmax+1,ip2,ip,jj)
       if (v6.gt.0.) then
       iechec=1
*        write (6,*) ' tetra vol 6 positif ',v6
        endif
        if (ip2.ne.ip4) then
        v7=vol(nptmax+1,ip4,jj,jp)
        if (v7.gt.0.) then
       iechec=1
*        write (6,*) ' tetra vol 7 positif ',v7
        endif
        v8=vol(nptmax+1,ip4,ip,jj)
        if (v8.gt.0.) then
       iechec=1
*        write (6,*) ' tetra vol 8 positif ',v8
        endif
       endif
        if (ip2.eq.ip4.and.ip1.eq.ip3.and.ip1.ne.ip2) then
        v9=vol(nptmax+1,ip2,ip1,ip)
        if (v9.gt.0.) then
       iechec=1
*        write (6,*) ' tetra vol 9 positif ',v9
        endif
        v10=vol(nptmax+1,ip1,ip2,jp)
        if (v10.gt.0.) then
       iechec=1
*        write (6,*) ' tetra vol 10 positif ',v10
        endif
       endif
        v11=vol(nptmax+1,ii,jj,ip)
        if (v11.gt.0.) then
       iechec=1
*        write (6,*) ' tetra vol 11 positif ',v11
        endif
        v12=vol(nptmax+1,jj,ii,jp)
        if (v12.gt.0.) then
       iechec=1
*        write (6,*) ' tetra vol 12 positif ',v12
        endif
       vtot=v1+v2+v3+v4+v5+v6+v7+v8+v9+v10+v11+v12
       if (iechec.eq.1) then
*  deplacement du point
        if (idepl.ne.0) then
         xyz(1,nptmax+1)=1.10*xyz(1,nptmax+1)-0.10*(xyz(1,ii)+xyz(1,jj))/2
         xyz(2,nptmax+1)=1.10*xyz(2,nptmax+1)-0.10*(xyz(2,ii)+xyz(2,jj))/2
         xyz(3,nptmax+1)=1.10*xyz(3,nptmax+1)-0.10*(xyz(3,ii)+xyz(3,jj))/2
         xyz(4,nptmax+1)=1.10*xyz(4,nptmax+1)-0.10*(xyz(4,ii)+xyz(4,jj))/2
*        call vcrit(nptmax+1)
*        write (6,*) ' tetra 1 deplacement du point ',nptmax+1
         idepl=idepl-1
         goto 356
        endif
*        goto 250
*           write (6,*) ' tetra echec 3 pt plus deplacable'
         return
       endif
*       write (6,*) ' tetra v* ',vtot,v1,v2,v3,v4,v5,v6,v7,v8,v9,v10
       if (v1.gt.vtot/1000) return
       if (v2.gt.vtot/1000) return
       if (ip1.ne.ip3.and.v3.gt.vtot/1000) return
       if (ip1.ne.ip3.and.v4.gt.vtot/1000) return
       if (v5.gt.vtot/1000) return
       if (v6.gt.vtot/1000) return
       if (ip2.ne.ip4.and.v7.gt.vtot/1000) return
       if (ip2.ne.ip4.and.v8.gt.vtot/1000) return
       if (ip2.eq.ip4.and.ip1.eq.ip3.and.ip1.ne.ip2.and.
     >      v9.gt.vtot/1000) return
       if (ip2.eq.ip4.and.ip1.eq.ip3.and.ip1.ne.ip2.and.
     >      v10.gt.vtot/1000) return
       if (v11.gt.vtot/1000) return
       if (v12.gt.vtot/1000) return
       NPTMAX=NPTMAX+1
       GOTO 352
250    CONTINUE
*   IL FAUT FAIRE DEUX ELEMENTS
       idepl=3
       inouv=1
       NFCMAX=NFCini
       NPTMAX=NPTMAX+1
*  CREATION POINT
       xyz(4,nptmax)=(xyz(4,ii)+xyz(4,jj))/2
*  deplacement du point pour l'eloigner de ii jj
       xn1=(xyz(2,jj)-xyz(2,ii))*(xyz(3,ip)-xyz(3,ii))-
     >     (xyz(3,jj)-xyz(3,ii))*(xyz(2,ip)-xyz(2,ii))
       yn1=(xyz(3,jj)-xyz(3,ii))*(xyz(1,ip)-xyz(1,ii))-
     >     (xyz(1,jj)-xyz(1,ii))*(xyz(3,ip)-xyz(3,ii))
       zn1=(xyz(1,jj)-xyz(1,ii))*(xyz(2,ip)-xyz(2,ii))-
     >     (xyz(2,jj)-xyz(2,ii))*(xyz(1,ip)-xyz(1,ii))
       sn1=sqrt(xn1**2+yn1**2+zn1**2)
       xn1=xn1/sn1
       yn1=yn1/sn1
       zn1=zn1/sn1
       xn2=(xyz(2,jp)-xyz(2,ii))*(xyz(3,jj)-xyz(3,ii))-
     >     (xyz(3,jp)-xyz(3,ii))*(xyz(2,jj)-xyz(2,ii))
       yn2=(xyz(3,jp)-xyz(3,ii))*(xyz(1,jj)-xyz(1,ii))-
     >     (xyz(1,jp)-xyz(1,ii))*(xyz(3,jj)-xyz(3,ii))
       zn2=(xyz(1,jp)-xyz(1,ii))*(xyz(2,jj)-xyz(2,ii))-
     >     (xyz(2,jp)-xyz(2,ii))*(xyz(1,jj)-xyz(1,ii))
       sn2=sqrt(xn2**2+yn2**2+zn2**2)
       xn2=xn2/sn2
       yn2=yn2/sn2
       zn2=zn2/sn2
*       xn=(xn1+xn2)/2
*       yn=(yn1+yn2)/2
*       zn=(zn1+zn2)/2
       xn=(xn1+xn2)/2
       yn=(yn1+yn2)/2
       zn=(zn1+zn2)/2
       sn=sqrt(xn**2+yn**2+zn**2)
       xn=xn/sn
       yn=yn/sn
       zn=zn/sn
*      xmil=(xyz(1,ii)+xyz(1,jj))/2
*      ymil=(xyz(2,ii)+xyz(2,jj))/2
*      zmil=(xyz(3,ii)+xyz(3,jj))/2
       xv=xyz(1,jj)-xyz(1,ii)
       yv=xyz(2,jj)-xyz(2,ii)
       zv=xyz(3,jj)-xyz(3,ii)
       sv=sqrt(xv**2+yv**2+zv**2)
       xv=xv/sv
       yv=yv/sv
       zv=zv/sv
       xli=xv*(xyz(1,ip)-xyz(1,ii))+yv*(xyz(2,ip)-xyz(2,ii))+
     >     zv*(xyz(3,ip)-xyz(3,ii))
       xlj=xv*(xyz(1,jp)-xyz(1,ii))+yv*(xyz(2,jp)-xyz(2,ii))+
     >     zv*(xyz(3,jp)-xyz(3,ii))
*      xl=(xli+xlj+2*sv+2*0)/6
       xl=0.5*sv
*      xl=(xli+xlj)/2
       xmil=xyz(1,ii)+xl*xv
       ymil=xyz(2,ii)+xl*yv
       zmil=xyz(3,ii)+xl*zv
       expf = xyz(4,nptmax)
       xyz(1,nptmax)=xmil+xn*expf*expfac
       xyz(2,nptmax)=ymil+yn*expf*expfac
       xyz(3,nptmax)=zmil+zn*expf*expfac
*      call vcrit(nptmax)
       IP1=0
       IP2=0
       IP3=0
       IP4=0
       goto 352
 355   CONTINUE
*  deplacement du point
         xyz(1,nptmax)=0.70*xyz(1,nptmax)+0.30*(xyz(1,ii)+xyz(1,jj))/2
         xyz(2,nptmax)=0.70*xyz(2,nptmax)+0.30*(xyz(2,ii)+xyz(2,jj))/2
         xyz(3,nptmax)=0.70*xyz(3,nptmax)+0.30*(xyz(3,ii)+xyz(3,jj))/2
         xyz(4,nptmax)=0.70*xyz(4,nptmax)+0.30*(xyz(4,ii)+xyz(4,jj))/2
*        call vcrit(nptmax)
*        write (6,*) ' tetra 2 deplacement du point ',nptmax
         idepl=idepl-1
 352   CONTINUE
       IPTT=NPTMAX
       CALL DIST(nptmax,nptaux,GL,IOK,ii,jj,ip,jp,ip1,ip2,ip3,ip4,0,0)
       IF (IOK.EQ.0) THEN
*        WRITE (6,*) ' tetra DIST ',nptaux
         if (nptaux.eq.0) then
           if (idepl.ne.0) goto 355
           NPTMAX=nptini
           NFCMAX=NFCini
*           write (6,*) ' tetra echec 4 dist et pt non deplacable'
           return
         endif
          if (idepl.ne.0) goto 355
         NPTMAX=nptini
         NFCMAX=NFCini
         iptt=nptaux
         idepl=0
         inouv=0
*        return
       ELSEIF (gl.lt.xyz(4,nptmax)/4) then
*        WRITE (6,*) ' tetra GL '
          if (idepl.ne.0) goto 355
         NPTMAX=nptini
         NFCMAX=NFCini
*           write (6,*) ' tetra echec 5 gl et pt non deplacable'
         return
       ENDIF
 354   continue
*  verification diago avant l'appel a tetra2
        iff=0
        lpts=0
        if (diago(iptt,ii,diacrd)) iff=ifdiag
        if (diago(iptt,jj,diacrd)) iff=ifdiag
        if (diago(iptt,ip,diacrd)) iff=ifdiag
        if (diago(iptt,jp,diacrd)) iff=ifdiag
*        write (6,*) ' tetra iff ',iff
        if (iff.ne.0) then
           nptmaa=nptmax
           nptmax=nptini
           if (nfc(4,iff).ne.0) return
        if (nfc(1,iff).eq.ip.or.nfc(2,iff).eq.ip.or.nfc(3,iff).eq.ip)
     >     then
            ipt1=isucc(iff,ip)
            if (ipt1.ne.ii.and.ipt1.ne.jj.and.
     >          ipt1.ne.jp.and.ipt1.ne.iptt) then
                  lpts=lpts+1
                  ipts(lpts)=ipt1
            endif
            ipt1=ipred(iff,ip)
            if (ipt1.ne.ii.and.ipt1.ne.jj.and.
     >          ipt1.ne.jp.and.ipt1.ne.iptt) then
                  lpts=lpts+1
                  ipts(lpts)=ipt1
            endif
         endif
        if (nfc(1,iff).eq.jp.or.nfc(2,iff).eq.jp.or.nfc(3,iff).eq.jp)
     >     then
            ipt1=isucc(iff,jp)
            if (ipt1.ne.ii.and.ipt1.ne.jj.and.
     >          ipt1.ne.ip.and.ipt1.ne.iptt) then
                  lpts=lpts+1
                  ipts(lpts)=ipt1
            endif
            ipt1=ipred(iff,jp)
            if (ipt1.ne.ii.and.ipt1.ne.jj.and.
     >          ipt1.ne.ip.and.ipt1.ne.iptt) then
                  lpts=lpts+1
                  ipts(lpts)=ipt1
            endif
         endif
        if (nfc(1,iff).eq.ii.or.nfc(2,iff).eq.ii.or.nfc(3,iff).eq.ii)
     >     then
            ipt1=isucc(iff,ii)
            if (ipt1.ne.ip.and.ipt1.ne.jp.and.
     >          ipt1.ne.jj.and.ipt1.ne.iptt) then
                  lpts=lpts+1
                  ipts(lpts)=ipt1
            endif
            ipt1=ipred(iff,ii)
            if (ipt1.ne.ip.and.ipt1.ne.jp.and.
     >          ipt1.ne.jj.and.ipt1.ne.iptt) then
                  lpts=lpts+1
                  ipts(lpts)=ipt1
            endif
         endif
        if (nfc(1,iff).eq.jj.or.nfc(2,iff).eq.jj.or.nfc(3,iff).eq.jj)
     >     then
            ipt1=isucc(iff,jj)
            if (ipt1.ne.ip.and.ipt1.ne.jp.and.
     >          ipt1.ne.ii.and.ipt1.ne.iptt) then
                  lpts=lpts+1
                  ipts(lpts)=ipt1
            endif
            ipt1=ipred(iff,jj)
            if (ipt1.ne.ip.and.ipt1.ne.jp.and.
     >          ipt1.ne.ii.and.ipt1.ne.iptt) then
                  lpts=lpts+1
                  ipts(lpts)=ipt1
            endif
         endif
        if (nfc(1,iff).eq.iptt.or.nfc(2,iff).eq.iptt.or.
     >                            nfc(3,iff).eq.iptt)  then
            ipt1=isucc(iff,iptt)
            if (ipt1.ne.ip.and.ipt1.ne.jp.and.
     >          ipt1.ne.ii.and.ipt1.ne.jj) then
                  lpts=lpts+1
                  ipts(lpts)=ipt1
            endif
            ipt1=ipred(iff,iptt)
            if (ipt1.ne.ip.and.ipt1.ne.jp.and.
     >          ipt1.ne.ii.and.ipt1.ne.jj) then
                  lpts=lpts+1
                  ipts(lpts)=ipt1
            endif
         endif
         lpts1=0
         lpts2=lpts
         do 710 lpt=1,lpts
          do 711 lpt1=1,lpts1
            if (ipts(lpt).eq.jpts(lpt1)) goto 710
 711      continue
          qualti(lpt)=min(qualt(ii,jj,ip,ipts(lpt)),
     >                    qualt(jj,ii,jp,ipts(lpt)))
          if (qualti(lpt).le.0.001d0) goto 710
          lpts1=lpts1+1
          jpts(lpts1)=ipts(lpt)
          qualtj(lpts1)=qualti(lpt)
 710     continue
         lpts=lpts1
         lpts1=0
         do 712 lpt=1,lpts
          do 713 l=1,lpts1
           if (qualtj(lpt).gt.qualti(l)) goto 714
 713      continue
          l=lpts1+1
 714      continue
          lpts1=lpts1+1
          do 715 l1=lpts1,l+1,-1
            qualti(l1)=qualti(l1-1)
            ipts(l1)=ipts(l1-1)
 715      continue
            qualti(l)=qualtj(lpt)
            ipts(l)=jpts(lpt)
 712     continue
         do 716 lpt=1,lpts
*          write (6,*) 'tetra appel tetra2 type 2 avec ',
*    >     ipts(lpt),iff,qualti(lpt)
            call tetra2(ii,jj,if1,if2,igagne,ipts(lpt))
            if (igagne.eq.1) return
 716     continue
*           write (6,*) ' tetra echec 6 lpts2 lpts',lpts2,lpts,ip,jp
*         return
          nptmax=nptmaa
        endif
*  test supplementaire
      if (inouv.eq.1) then
       xyz(1,nptmax+1)=1.90*xyz(1,nptmax)-0.90*(xyz(1,ii)+xyz(1,jj))/2
       xyz(2,nptmax+1)=1.90*xyz(2,nptmax)-0.90*(xyz(2,ii)+xyz(2,jj))/2
       xyz(3,nptmax+1)=1.90*xyz(3,nptmax)-0.90*(xyz(3,ii)+xyz(3,jj))/2
       xyz(4,nptmax+1)=1.90*xyz(4,nptmax)-0.90*(xyz(4,ii)+xyz(4,jj))/2
*      call vcrit(nptmax+1)
       IF (.NOT.IN2(ii,nptmax+1,nfcini)) THEN
*  deplacement du point
        if (idepl.ne.0) goto 355
        NPTMAX=nptini
        NFCMAX=NFCini
*       write (6,*) ' tetra 2eme test in2 echoue ',ii,jj,iptt
        RETURN
       ENDIF
      endif
            quala=qualt(ii,jj,ip,iptt)
            qualb=qualt(jj,ii,jp,iptt)
*           write (6,*) ' tetra quala qualb ',quala,qualb
       if (quala.gt.0.001.and.qualb.gt.0.001)
     >       call tetra2(ii,jj,if1,if2,igagne,iptt)
       if (igagne.eq.1) return
*  deplacement du point
        if (idepl.ne.0) goto 355
       nptmax=nptini
       idjf=1
       goto 152
150    CONTINUE
       idjf=0
C
C      LE VOLUME CREE EST INVALIDE: IL FAUT DONC DETRUIRE LES FACETTES
C      CREEES. -------------------------------------------------------
       CALL REPSUB(IF4)
       CALL REPSUB(IF3)
       CALL REPSUB(IF2)
       CALL REPSUB(IF1)
C
       NFCMAX=NFCini
 152   CONTINUE
C
*   on essaye en cas d'intersection de facette d'utiliser cette information
*      call facetx(k1,k2,k3,ifr)
       k1=iirfac
       k2=j1rfac
       k3=j2rfac
       if (k1.gt.nptmax) k1=0
       if (k2.gt.nptmax) k2=0
       if (k3.gt.nptmax) k3=0
       ifr=ifrfac
*       write (6,*) ' tetra retour de facetx ',k1,k2,k3,ifr
       if (k1.eq.0) then
*          write (6,*) ' tetra echec 7 facetx => 0 '
          if (ifdiag.eq.0) return
          k1=nfc(1,ifdiag)
          k2=nfc(2,ifdiag)
          k3=nfc(3,ifdiag)
       endif
       if (ip.eq.k1.or.jp.eq.k1) then
         if (ifr.eq.if3.or.ifr.eq.0) then
         iptt=0
         if (jj.eq.k2) iptt=k3
         if (jj.eq.k3) iptt=k2
         if (iptt.ne.0) call tetra2(ii,jj,if1,if2,igagne,iptt)
         if (igagne.eq.1) return
         endif
         if (ifr.eq.if4.or.ifr.eq.0) then
         iptt=0
         if (ii.eq.k2) iptt=k3
         if (ii.eq.k3) iptt=k2
         if (iptt.ne.0) call tetra2(ii,jj,if1,if2,igagne,iptt)
         if (igagne.eq.1) return
         endif
       endif
       if (ip.eq.k2.or.jp.eq.k2) then
         if (ifr.eq.if3.or.ifr.eq.0) then
         iptt=0
         if (jj.eq.k1) iptt=k3
         if (jj.eq.k3) iptt=k1
         if (iptt.ne.0) call tetra2(ii,jj,if1,if2,igagne,iptt)
         if (igagne.eq.1) return
         endif
         if (ifr.eq.if4.or.ifr.eq.0) then
         iptt=0
         if (ii.eq.k1) iptt=k3
         if (ii.eq.k3) iptt=k1
         if (iptt.ne.0) call tetra2(ii,jj,if1,if2,igagne,iptt)
         if (igagne.eq.1) return
         endif
       endif
       if (ip.eq.k3.or.jp.eq.k3) then
         if (ifr.eq.if3.or.ifr.eq.0) then
         iptt=0
         if (jj.eq.k1) iptt=k2
         if (jj.eq.k2) iptt=k1
         if (iptt.ne.0) call tetra2(ii,jj,if1,if2,igagne,iptt)
         if (igagne.eq.1) return
         endif
         if (ifr.eq.if4.or.ifr.eq.0) then
         iptt=0
         if (ii.eq.k1) iptt=k2
         if (ii.eq.k2) iptt=k1
         if (iptt.ne.0) call tetra2(ii,jj,if1,if2,igagne,iptt)
         if (igagne.eq.1) return
         endif
       endif
         if (idjf.eq.1) then
*           write (6,*) ' tetra echec 8 apres facetx',k1,k2,k3,ifr
            return
         endif
         if (npf(5,ii).ne.0.or.npf(5,jj).ne.0) goto 250
*        write (6,*) ' tetra diago 4-4 => point milieu',ip,jp
          goto 350
       RETURN
151    CONTINUE
       CALL REPSUB(IF4)
       CALL REPSUB(IF3)
       CALL REPSUB(IF2)
       CALL REPSUB(IF1)
       NFCMAX=NFCini
         if (iptt.ne.0) call tetra2(ii,jj,if1,if2,igagne,iptt)
         if (igagne.eq.1) return
         if (npf(5,ii).ne.0.or.npf(5,jj).ne.0) goto 250
*        write (6,*) ' tetra diago 4-4 => point milieu',ip,jp
          goto 350
       return
 155   CONTINUE
*       write (6,*) ' tetra probleme avec nouvelle arrete '
*      CALL COM33(II,JJ,IF1,IF2,IGAGNE)
C
*       write (6,*) ' tetra echec 9'
       RETURN
       END