      subroutine mshptg(cr,h,c,nu,nbs,nbsmx,tri,arete,nba,sd,nbsd,reft
     &                 ,nbt,coef,puis,iop,err,ilang)
c-----------------------------------------------------------------------
c      but:  construire une triangulation a partir d'un ensemble de
c             points et d'un maillage frontalier
c-----------------------------------------------------------------------
c     entre :
c     -------
c           cr(2,nbsmx)  tableau des coordonnees des nbs points donnes 
c
c           h (nbsmx)    tableau du h local voulu autour de chaque point
c                          donnes 
c                                
c           nbs          nombre de points donnes
c           nbsmx        nombre de points maximal a cree 
c                        si nbs  = nbsmx ont ne cree pas de points
c                        si nbs  < nbsmx => erreur
c
c           arete(2,nba) tableau des aretes du maillage a forcer
c                          exemple :la frontiere 
c
c           nba          le nombre d'aretes du maillage
c
c           sd(2,nbsd)   tableau definisant les nbsd  sous domaine
c                          (reference des triangles gerener)
c                          abs(sd(1,i)) =  numero d'une l'arete 
c                          si sd(1,i) est positive alors le sous domaine
c                          est a gauche de l'arete sinon il est a droite
c                          sd(2,i) donne le numero du sous doimaine 
c
c           iop          option de trace et d'impression
c                        1) impre = mod(iop,10)
c                           impre = 0 => pas impression
c                           impre > 4 => le resultat est imprime
c                           impre > 5 => debug
c                        2) iop  <10 ==> pas de graphique
c                            iop > 9  on trace le  mailage construit
c                            iop > 29  + numero des points
c                            iop > 49  on  trace  tout les mouvement 
c                                      des triangles
c                            iop > 99  un (cr) est attendu apres tracer
c                                      de triangle
c                            le graphique utilise fortran-3d
c                            -------------------------------
c
c           puis         coefficent de generation des points 
c                        .1  => on propage plus loin les rafinement 
c                               donnes par h 
c                        .25 => valeur conseillee
c
c           coef         coefficent sur le test arret 
c                          le valeur conseillee est .75
c                          remarque le nombre de point genere est en 
c                          O(coef**2)
c
c        tableaux de travail:
c        --------------------
c
c           c(2,nbsmx)    tableau d'entiers (copie de coordonnees)
c           tri(ltri)     tableau d'entiers 
c        
c        out :
c        -----
c
c         nbs         nombre de points   donnes + generes 
c         nbt         nombre de triangles generes 
c         cr(1:2,nbs) coordonnees des sommets donnes + generes
c         nu(1:3,nbt) sommets des triangles (tableau des connections)
c                       telle que les sommets tourne dans le sens direct
c         reft(1:nbt) numero de sous domaine de chaque triangle 
c
c         err    si err = 0 alors pas de probleme
c                sinon nbt = 0 et pas de triangulation
c
c     dimension des tableaux
c     ----------------------
c     definition des parameters
c     nbtmx = 2*(nbs-1) ,  ltri = max(4*nbs+2*nbsd,nba)
c
c     integer : nu(6*nbtmx) , reft(nbtmx) , c(2*nbsmx) , tri(ltri) 
c     integer : arete(2,nba), sd(2,nbsd)
c     real    : cr(2*nbsmx) , h(nbsmx)
c
c----------------------------------------------------------------------
c     programmeur F.Hecht, INRIA Rocquencourt, 78153 le chesnay, France 
c           version 1.0  mars 1986     
c-----------------------------------------------------------------------
c    DCL
      integer nbs,nbt,nba,nbsmx,nbsd
      integer c(2,nbsmx),tri(*),arete(2,nba),sd(2,nbsd)
      integer nu(12*(nbsmx-1)),reft(2*(nbsmx-1))
      integer iop,err
      real cr(2,nbsmx),h(nbsmx),coef,puis
c...... dcl des variables locales ..........................................
      integer i,j,k,t,tete,impre,nbtgrn,nbsgrn
      real trfri(4)
c...... dcl des parametres de la regularisations` 
      integer itermx
      real omega,eps
Cfj      parameter (itermx=20,omega = 1.4,eps = .005)
      parameter (itermx=20,omega = 1.4,eps = 1.e-7)
c...... fin dens dcl....................................................
      impre=mod(iop,10)
      err = 0
      nbt = 0
      if(nbs.lt.3.or.nbsmx.lt.nbs) then
        err = 1
        print *,'*** fatal error mshpts : le nombre de points ',nbs
     &         ,' est < 3 ou > ',nbsmx,' nb de maximal '
        return
      endif
c-------------------------
c preparation des donnees
c-------------------------
      call mshtri (cr,c,nbs,tri,tri(nbs+1),trfri,iop,err,ilang)
      if (err.ne.0) then
        print*,'*** ERREUR mshtri'
        return
      endif

c--------------------------------
c maillage de l enveloppe convexe
c--------------------------------
      call mshcxi (c,nu,tri,nbs,tete,iop,err)
c-----------------------------------------------------------------------
c     definition de tableau nu(1:6,2*nbs-2)
c-----------------------------------------------------------------------
c     nu(*,ie) definit soit un element ,soit un sommet frontiere
c     si nu(5:6,ie) = (0,0) alors ie est un sommet frontiere
c     avec nu(1,ie) = numero du sommet
c          nu(2,ie) = 8*t + a
c                     ou t est le numero du triangle ayant l'arete 
c                     frontiere (a) dont le premier sommet est nu(1,ie)
c          nu(3,ie) = pointeur dans nu sur sommet frontiere precedent
c          nu(4,ie) = pointeur dans nu sur sommet frontiere suivant
c
c     sinon ie est un element :
c          nu(1:3,ie) numero des 3 sommets du triangle ie tournant dans
c                     le sens direct
c          nu(4:6,ie) = (d4,d5,d6) donnee des 3 aretes ai
c           ai est forme des sommets nu(i-3,ie),nu(mod(i,3)+1,ie)
c           si di < 0 alors arete i est frontiere et -di est pointeur
c             sur 1er sommet frontiere de i
c           sinon arete est interne et di = 8*ta + ata
c              ou ta est le numero du triangle adjacent a l'arete
c              et ata est le numero de l'arete dans ta
c------------------------------------------------------------------------
      if (err.ne.0) then
        print*,'*** ERREUR mshcxi'
        return
      endif
c
      do 10 i=1,nbs
        tri(i)=0
10    continue
c      if(impre.gt.4)  print *,'frontiere convexe tete =',tete
      i=tete
20    continue
      if(impre.gt.4) then
c        print *,i,' s = ',nu(6*(i-1)+1),' t = ',nu(6*(i-1)+2)/8
c     &         ,' a= ',mod(6*(i-1)+2,8),' precedent =',nu(6*(i-1)+3)
      endif
      j=nu(6*(i-1)+4)
      tri(nu(6*(i-1)+1))=nu(6*(j-1)+1)
      i=j
      if(i.ne.tete) goto 20
c-----------------------------
c traitement frontiere
c-----------------------------
c      if(impre.ge.4)  print *,' les elements '
      k=0
      call mshfrt(c,nu,nbs,arete,nba,sd,nbsd,reft,tri,iop,err)
      if (err.ne.0) then
        print*,'*** ERREUR mshfrt'
        return
      endif
c-------------------------------------------------------------------
c       on a modifie nu les sommets frontiere n'ont plus de sens
c       ainsi que les pointeurs sur ces elements
c-------------------------------------------------------------------
      nbsgrn = nbs
      call mshgpt (c,cr,nu,h,reft,nbsgrn,nbsmx,nbtgrn,coef,puis,trfri
     &            ,iop,err)
      if (err.ne.0) then
        print*,'*** ERREUR mshgpt'
        return
      endif
c
c     construction du tableau nu(1:3,1:nbt)
c------------------------------------------
      nbt=0
      k = 0
      j = 1
      do 200 t=1,nbtgrn
       if(nu(j+5).ne.0) then
         nbt=nbt + 1
         reft(nbt) = reft(t)
         do 190 i=0,2
          k=k+1
          nu(k)=nu(j+i)
190      continue
         if(iop.ge.10) call mshdrw (c,nu,3,nbt,iop)
c         if(impre.gt.4) print '(7i12)',t,(nu(j+i),i=0,5)
       endif
       j = j + 6
200   continue
c     dans nu il y a (s1(t),s2(t),s3(t),t=1,nbt)
c     ou s1 s2 s3 sont les 3 sommets de t
c------------------------------------------------
      do 210 i =1,nbs
       tri(i)=1
210   continue
      do 220 i=nbs+1,nbsgrn
       tri(i)=0
220   continue
      call mshvoi (nu,tri(nbsgrn+1),nu(nbt*3+1),nbt,nbsgrn)
      call mshrgl (cr,tri,nbsgrn,nu,tri(nbsgrn+1),nu(nbt*3+1),nbt
     &            ,omega,itermx,eps,iop)

      nbs = nbsgrn
      end
c**********************************************************************
      subroutine mshvoi (nu,w1,w,nbt,nbs)
c
c RECHERCHE DU VOISINAGE
C-----------------------
      integer nbt,nbs,nu(3*nbt),w1(0:nbs),w(3*nbt),i,is
      do 10 i=1,nbs
       w1(i)=0
10    continue
      do 30 i=1,3*nbt
        w1(nu(i)) = w1(nu(i)) + 1
30    continue
      w1(0) = 0
      do 40 i=1,nbs
        w1(i)= w1(i-1) + w1(i)
40    continue
      do 60 i=1,nbt*3
         is = nu(i) -1
         w1(is)    = w1(is) + 1
         w(w1(is)) = i
60    continue
      do 70 i = nbs,1,-1
        w1(i)=w1(i-1)
70    continue
      w1(0) = 0
      end
c**********************************************************************
      subroutine mshrgl (c,nrfs,nbs,nu,w1,w,nbt,omega,itermx,eps,iop)
C
C REGULARISATION PAR MOYENNE BARYCENTRIQUE
C-----------------------------------------
      integer itermx
      integer nbs,nbt,nrfs(nbs),nu(3*nbt),w1(0:nbs),w(3*nbt),iop
      real bx,by,err,omega,depx,depy,eps,dx,xmin,xmax,ymin,ymax
      integer iter,i1,i2,i,is,impre,k,ic
      real c(2,nbs)
cc      impre = mod(iop,10)
      xmin=c(1,1)
      ymin=c(2,1)
      xmax=c(1,1)
      ymax=c(2,1)
      do 10 ic=1,nbs
       xmin=min(c(1,ic),xmin)
       ymin=min(c(2,ic),ymin)
       xmax=max(c(1,ic),xmax)
       ymax=max(c(2,ic),ymax)
10    continue
      dx = max(xmax-xmin,ymax-ymin)
      do 50 iter=1,itermx
       err = 0
       i2 = w1(0)
       do 40 is=1,nbs
        i1 = i2 + 1
        i2 = w1(is)
        if(i2.ge.i1.and.nrfs(is).eq.0) then
          bx = 0
          by = 0
          do 20 i=i1,i2
           if(mod(w(i),3).eq.0) then
            k = w(i) - 2
           else
            k = w(i) + 1
           endif
           bx = bx + c(1,nu(k))
           by = by + c(2,nu(k))
20        continue
          bx = bx / (i2-i1+1)
          by = by / (i2-i1+1)
          depx=omega*(c(1,is)-bx)
          depy=omega*(c(2,is)-by)
          c(1,is)=c(1,is)-depx
          c(2,is)=c(2,is)-depy
          err = max(err,abs(depx),abs(depy))
        endif
40     continue
c       if(impre.ne.0) print *,' iteration ',iter
c     &                       ,' deplacement maximal : ',err,err/dx
c--------------------------------     
       if(err.le.eps*dx) then
c------------graphique------trace a la fin de la regularisation---
      return     
      end if
50    continue
c      print *,' warning mshrlg : on a pas converge en ',itermx
c     &       ,' iterations'
      end
c**********************************************************************
Cfj      subroutine mshdr1 (c,nu,i6,t,iop)
CfjC
CfjC ROUTINE DE TRACE D UN TRIANGLE
CfjC-------------------------------
Cfj      integer i6,nu(i6,*),t,i,i1,i2,i3,iop
Cfj      real x,y,c(2,*)
Cfj      i1=nu(1,t)
Cfj      i2=nu(2,t)
Cfj      i3=nu(3,t)
Cfj      if(iop.ge.100) pause
Cfj      end
c**********************************************************************
      subroutine mshgpt (c,cr,nu,h,reft,nbs,nbsmx,nbt,coef,puis,trfri
     &                  ,iop,err)
      integer nbs,nbsmx,nbt,iop,err
      integer c(2,nbsmx),nu(6,2*(nbsmx-1)),reft(2*(nbsmx-1))
      real cr(2,nbsmx),h(nbsmx),coef,puis,trfri(4)
      integer tete,t,s1,s2,s3,itera,nbsold,ix,iy
      real x,y,aire,hs,det,h1,h2,h3,pui
c      logical impre,impr1,impr2
c      impre = mod(iop,10).ne.0
c      impr1 = mod(iop,10).ge.2
c      impr2 = mod(iop,10).ge.8
      pui = puis
c      if(impre) print*,'mshgpt:  nbs = ',nbs,' nbsmx = ',nbsmx
      nbt = 2*nbs-2
      if(nbs.ge.nbsmx) return
      tete = 0
c     initialisation de la liste des triangles libre      
      do 10 t=1,nbt
       if(nu(6,t).eq.0) then
        nu(1,t)=tete
        tete = t
       endif
10    continue
      itera = 0
20    continue
      itera = itera + 1
      nbsold = nbs 
      do 100 t=1,nbt
       if(nu(6,t).ne.0) then
         s1=nu(1,t)
         s2=nu(2,t)
         s3=nu(3,t)
c        calcul de 2 fois l'aire du triangle 
         det =   ( cr(1,s2) - cr(1,s1) ) * ( cr(2,s3) - cr(2,s1) )
     &         - ( cr(2,s2) - cr(2,s1) ) * ( cr(1,s3) - cr(1,s1) )
         aire = det*coef
         if(puis.gt.0.) then
           hs = ((h(s1)**pui + h(s2)**pui + h(s3)**pui )/3.)**(1./pui )
         elseif(puis.gt.-1.)then
           hs = (h(s1)*h(s2)*h(s3))**(1./3.) 
         elseif(puis.gt.-2) then
           hs = 3.*h(s1)*h(s2)*h(s3)
     &        / ( h(s1)*h(s2) + h(s1)*h(s3) + h(s2)*h(s3) )
         else
           hs = sqrt(3.)*(h(s1)*h(s2)*h(s3)
     &        / sqrt(h(s1)*h(s2))**2+(h(s1)*h(s3))**2+(h(s2)*h(s3))**2)
         endif
         if(aire.gt.hs*hs) then
           h1 = 1.
           h2 = 1.
           h3 = 1.
           x = (cr(1,s1)*h1 + cr(1,s2)*h2 + cr(1,s3)*h3 )/(h1+h2+h3)
           y = (cr(2,s1)*h1 + cr(2,s2)*h2 + cr(2,s3)*h3 )/(h1+h2+h3)

           ix = nint( trfri(1) * ( x -trfri(2) )            )
           iy = nint( trfri(1) * ( y -trfri(3) ) - trfri(4) )

           if(     (c(1,s2)-     ix)*(c(2,s3)-iy     ) 
     &            -(c(2,s2)-     iy)*(c(1,s3)-ix     ).le.0 
     &        .or. (ix     -c(1,s1))*(c(2,s3)-c(2,s1)) 
     &            -(iy     -c(2,s1))*(c(1,s3)-c(1,s1)).le.0 
     &        .or. (c(1,s2)-c(1,s1))*(iy     -c(2,s1)) 
     &            -(c(2,s2)-c(2,s1))*(ix     -c(1,s1)).le.0  ) then
c             if(impr1) then
c               print *,'warning mshgpt: le point genere n''est pas dans'
c     &                ,' le triangle',nbs,' h =',h1,h2,h3,' h(s) =',hs
c               print *,' s1  ',c(1,s1),c(2,s1),' s2  ',c(1,s2),c(2,s2)
c     &                ,' s3  ',c(1,s3),c(2,s3),' nbs ',ix,iy
c             endif
           else
            if(nbs.ge.nbsmx) then
c             print *,' warning mshgpt :on ne peut plus cree de points '
c             print *,' nb de point = ',nbs,' nb de point maximal=',nbsmx
             return
            endif
            nbs = nbs +1
            c(1,nbs) = ix
            c(2,nbs) = iy
            cr(1,nbs) =  ix            /trfri(1) + trfri(2)
            cr(2,nbs) = (iy + trfri(4))/trfri(1) + trfri(3)
            h(nbs)    = hs
c            if(impr2) print *,'mhsgpt: triangle ',t
c     &               ,' generation du point',nbs,cr(1,nbs),cr(2,nbs)
             call msha1p (t,nbs,c,nu,reft,tete,nbt,nbs,iop,err)
             if (err.ne.0) then
               print*,'*** ERREUR msha1p'
               return
             endif
           endif
         endif
       endif
100   continue
c      if(impr1) then 
c        print *,' iterations de generation des point',itera
c     &         ,' nb de points ',nbs
c      endif
      if(nbsold.ne.nbs) goto 20
c      if(impre) then
c        print *,'mshptg : nb d''iteration:',itera,'et nb de points',nbs
c        print *,'------------------------------------------------------'
c      endif
      end
c**********************************************************************
      subroutine msha1p (t,s,c,nu,reft,tete,nbt,nbs,iop,err)
      integer t,s,nbt,nbs,tete,iop,err
      integer c(2,nbs),nu(6,2*(nbs-1)),reft(2*(nbs-1))
      integer t1,t2,t3,ta2,ta3,ia2,ia3,tta,ia,i
      t1 = t
c      if(mod(iop,10).ge.9) then
c        print *,' t old',t,(nu(i,t1),i=1,6)
c      endif
      if(tete.eq.0) then
       nbt= nbt+1
       t2=nbt
      else
       t2=tete
       tete=nu(1,tete)
      endif
      if(tete.eq.0) then
       nbt= nbt+1
       t3=nbt
      else
       t3=tete
       tete=nu(1,tete)
      endif
      nu(1,t2) = s
      nu(2,t2) = nu(2,t)
      nu(3,t2) = nu(3,t)
      nu(4,t2) = 5 + (2**3)*t1
      nu(5,t2) = nu(5,t)
      nu(6,t2) = 5 + (2**3)*t3

      nu(1,t3) = nu(1,t)
      nu(2,t3) = s
      nu(3,t3) = nu(3,t)
      nu(4,t3) = 6 + (2**3)*t1
      nu(5,t3) = 6 + (2**3)*t2
      nu(6,t3) = nu(6,t)

      tta = nu(5,t)
      if(tta.gt.0) then
       ta2 = tta/(2**3)
       ia2 = tta - (2**3)*ta2
       nu(ia2,ta2) = 5 + (2**3)*t2
c       if(mod(iop,10).ge.9)print *,' ta2 ',ta2,ia2,(nu(i,ta2),i=1,6)
      endif
      tta = nu(6,t)
      if(tta.gt.0) then
       ta3 = tta/(2**3)
       ia3 = tta - (2**3)*ta3
       nu(ia3,ta3) = 6 + (2**3)*t3
c       if(mod(iop,10).ge.9)print *,' ta3 ',ta3,ia3,(nu(i,ta3),i=1,6)
      endif

      nu(3,t1) = s
      nu(5,t1) = 4 + (2**3)*t2
      nu(6,t1) = 4 + (2**3)*t3

      reft(t2) = reft(t)
      reft(t3) = reft(t)

      if(iop.ge.50) call mshdrw (c,nu,6,t1,iop)
      if(iop.ge.50) call mshdrw (c,nu,6,t2,iop)
      if(iop.ge.50) call mshdrw (c,nu,6,t3,iop)
c      if(mod(iop,10).ge.9) then
c        print *,' t1 ',t1,(nu(i,t1),i=1,6)
c        print *,' t2 ',t2,(nu(i,t2),i=1,6)
c        print *,' t3 ',t3,(nu(i,t3),i=1,6)
c      endif
      call mshopt (c,nu,t1,4,nbs,iop,err)
      if (err.ne.0) then
        print*,'*** ERREUR mshopt'
        return
      endif
      call mshopt (c,nu,t2,5,nbs,iop,err)
      if (err.ne.0) then
        print*,'*** ERREUR mshopt'
        return
      endif
      call mshopt (c,nu,t3,6,nbs,iop,err)
      if (err.ne.0) then
        print*,'*** ERREUR mshopt'
        return
      endif
      if(err.ne.0) return
      end
c**********************************************************************
      integer function mshlcl(c,nu,tete,s,nbs)
      integer nbs,c(2,nbs),nu(6,nbs+nbs-2),tete,s
      integer x,y,sf,pt,ppt,det
      logical init
      x=c(1,s)
      y=c(2,s)
      init=.true.
      pt=tete
10    continue
      ppt=pt
      pt=nu(4,pt)
      if(pt.ne.tete) then
        det=x*c(2,nu(1,pt)) -y*c(1,nu(1,pt))
        if(det.lt.0) then
          init=.false.
          goto 10
        elseif(init.and.det.eq.0) then
          goto 10
        endif
      endif
      mshlcl=ppt
      end
c**********************************************************************
      subroutine mshdrw (c,nu,i6,t,iop)
      integer i6,c(2,*),nu(i6,*),t,i,i1,i2,i3,iop
      real x,y
      if(iop.ge.100) pause
      end
c**********************************************************************
      subroutine mshtri (cr,c,nbs,tri,nu,trfri,iop,err,ilang)
      integer nbs,c(1:2,1:nbs),tri(1:nbs),nu(1:nbs),iop,err
      real cr(1:2,1:nbs),trfri(4)
      integer iii,ic,xx,ip,i,j,jc,k,trik,tri3,det,ierr
      real aa1,aa2,xmin,xmax,ymin,ymax,xx1,xx2,yy1,yy2,dxx
      real precis
      parameter (precis=2.**15-1.)
      ierr = 0
      iii=1
      xmin=cr(1,1)
      ymin=cr(2,1)
      xmax=cr(1,1)
      ymax=cr(2,1)
      do 10 ic=1,nbs
       xmin=min(cr(1,ic),xmin)
       ymin=min(cr(2,ic),ymin)
       xmax=max(cr(1,ic),xmax)
       ymax=max(cr(2,ic),ymax)
       tri(ic)=ic
       if(cr(1,ic).lt.cr(1,iii)) then
        iii=ic
       endif
10    continue
      aa1 = precis/(xmax-xmin)
      aa2 = precis/(ymax-ymin)
      aa1 = min(aa1,aa2)
      aa2 = aa1*(cr(2,iii)-ymin)
      trfri(1)=  aa1
      trfri(2) = cr(1,iii)
      trfri(3) = ymin
      trfri(4) = aa2
      do 20 ic=1,nbs
       c(1,ic) = nint(aa1*(cr(1,ic)-cr(1,iii)))
       c(2,ic) = nint(aa1*(cr(2,ic)-ymin)-aa2)
       nu(ic)= c(1,ic)**2 + c(2,ic)**2
20    continue
c----------------------------------------------------------
      call mshtr1 (nu,tri,nbs)
      ip = 1
      xx=nu(ip)
      do 30 jc=1,nbs
       if(nu(jc).gt.xx)then
        call mshtr1 (nu(ip),tri(ip),jc-ip)
        do 25 i=ip,jc-2
          if (nu(i).eq.nu(i+1)) then
            ierr=ierr+1
            if (ilang.eq.0) then
              print*,'*** points confondus :',tri(i),tri(i+1)
            else
              print*,'*** merged points:',tri(i),tri(i+1)
            endif
          endif
 25     continue
        xx=nu(jc)
        ip=jc
       endif 
       ic=tri(jc)
       nu(jc)=c(2,ic)
30    continue        
      call mshtr1 (nu(ip),tri(ip),nbs-ip)
      do 35 i=ip,jc-2
        if (nu(i).eq.nu(i+1)) then
          ierr=ierr+1
          if (ilang.eq.0) then
            print*,'*** points confondus :',tri(i),tri(i+1)
          else
            print*,'*** merged points:',tri(i),tri(i+1)
          endif
        endif
 35   continue
      if (ierr.ne.0) then
        if (ilang.eq.0) then
          print*,'*** ERREUR: il y a des points confondus (',ierr,')'
        else
          print*,'*** ERROR: there are some merged points (',ierr,')'
        endif
        err = 2
        return
      endif
      k=2
50    continue
      if(k.le.nbs) then
        k=k+1
        det = c(1,tri(2))*c(2,tri(k)) - c(2,tri(2))*c(1,tri(k))
        if(det.eq.0) goto 50
      else
c        print *,'fatal error mshtri tous les points sont alignes'
c        print *,'tri =',(tri(k),k=1,nbs)
        err = 3
        return
c        stop 'fatal error'
      endif
c     k est le premier point non aligne        
      trik = tri(k)
      do 60 j=k-1,3,-1
        tri(j+1)=tri(j)
60    continue
      tri(3)=trik
      if(det.lt.0) then
c       on inverse les  points 2 3 tries
        tri3=tri(3)
        tri(3)=tri(2)
        tri(2)=tri3
      endif
      end
c**********************************************************************
      subroutine mshtr1 (criter,record,n)
      integer record(n)
      integer criter(n)
c     trie selon les valeurs de criter croissantes
c     record suit le reordonnancement
c
      integer i,l,r,j,n
      integer rec
      integer crit
c     
      if(n.le.1) return
      l=n/2+1
      r=n
2     if(l.le.1)goto 20
      l=l-1
      rec=record(l)
      crit=criter(l)
      goto 3
20    continue
      rec=record(r)
      crit=criter(r)
      record(r)=record(1)
      criter(r)=criter(1)
      r=r-1
      if(r.eq.1)goto 999
3     j=l
4     i=j
      j=2*j
      if(j-r)5,6,8
5     if(criter(j).lt.criter(j+1))j=j+1
6     if(crit.ge.criter(j))goto 8
      record(i)=record(j)
      criter(i)=criter(j)
      goto 4
8     record(i)=rec
      criter(i)=crit
      goto 2
999   record(1)=rec
      criter(1)=crit
      return
      end
c**********************************************************************
      subroutine mshcvx(direct,c,nu,pfold,nbs,iop,err)
      integer nbs,c(2,nbs),nu(6,nbs+nbs-2),pfold,iop,err
      logical direct
      integer pp,ps,i1,i2,i3,i4,i5,i6
      integer pf,psf,ppf,s1,s2,s3,t,t4,t5,a4,a5,det,tt4,tt5
      if(direct) then
        pp=3
        ps=4
        i1=1
        i2=3
        i3=2
        i4=6
        i5=5
        i6=4
      else
        pp=4
        ps=3
        i1=1
        i2=2
        i3=3
        i4=4
        i5=5
        i6=6
      endif
10    continue
      ppf=pfold
      pf =nu(ps,pfold)
      psf=nu(ps,pf)
      s1=nu(1,ppf)
      s2=nu(1,pf)
      s3=nu(1,psf)
      det =   ( c(1,s2) - c(1,s1) ) * ( c(2,s3) - c(2,s1) )
     &      - ( c(2,s2) - c(2,s1) ) * ( c(1,s3) - c(1,s1) )
c      print *,' mshcvx convexification de ', s1,s2,s3,det,direct
      if(((.not.direct).and.det.gt.0).or.(direct.and.det.lt.0)) then
c       on ajoute un triangle t et on detruit une arete
c       -----------------------------------------------
        if(direct) then
         tt4 = nu(2,ppf)
         tt5 = nu(2,pf)
        else
         tt4 = nu(2,pf)
         tt5 = nu(2,psf)
        endif
        t4 = tt4/(2**3)
        t5 = tt5/(2**3)
        a4 = tt4 -8 * t4
        a5 = tt5 -8 * t5
c       destruction de l'arete frontiere en pf 
c       --------------------------------------
        nu(ps,ppf) = psf
        nu(pp,psf) = ppf
c       on remplace l'arete frontiere par l'element genere 
c       ---------------------------------------------------
        t = pf
c       update de l'arete non detruite
c       ------------------------------
        if(direct) then
          nu(2,ppf) = (2**3) * t + i6
        else
          nu(2,psf) = (2**3) * t + i6
        endif
c       on cree l'element
c       -----------------
        nu(i1,t ) = s1
        nu(i2,t ) = s2
        nu(i3,t ) = s3
        nu(i4,t ) = (2**3) * t4 + a4
        nu(i5,t ) = (2**3) * t5 + a5
        if(direct) then
          nu(i6,t ) = -ppf
        else
          nu(i6,t ) = -psf
        endif
        nu(a4,t4) = (2**3) * t + i4
        nu(a5,t5) = (2**3) * t + i5
        if(iop.ge.50) call mshdrw (c,nu,6,t,iop)
        call mshopt (c,nu,t5,a5,nbs,iop,err)
        if (err.ne.0) then
          print*,'*** ERREUR mshopt'
          return
        endif
        goto 10
      endif
      end
c**********************************************************************
      subroutine mshcxi (c,nu,tri,nbs,tete,iop,err)
      integer nbs,c(2,nbs),nu(6,2*nbs-2),tri(nbs),tete,iop
      integer mshlcl,err
      integer i,j,s,t,pf,ppf,psf,npf,pp,ps,taf,iaf,free,ttaf
      parameter (pp=3,ps=4)
c     initialisation de la free liste dans nu
      do 10 i=1,nbs+nbs-2
       nu(1,i)=i+1
       do 10 j=2,6
        nu(j,i)=0
10    continue
      nu(1,nbs+nbs-2)=0
      free = 1
c     initialisation du premier triangle
      t=free
      free = nu(1,free)
c     initialisation de la liste frontiere
      tete=free
      pf  =free
      do 20 i=1,3
       nu(i  ,t) = tri(i)
       nu(3+i,t) = -pf
       ppf       = pf
       free      = nu(1,pf)
       pf        = free
       if(i.eq.3) pf=tete
       nu(1,ppf) = tri(i)
       nu(2,ppf) = i + 3 + (2**3) * t
       nu(ps,ppf) = pf
       nu(pp,pf ) = ppf
20    continue      
      if(iop.ge.50) call mshdrw (c,nu,6,t,iop)
c      print *,' free =',free,' nu ='
c      print '(6(i12))',((nu(i,j),i=1,6).j=1,4)
      do 30 i=4,nbs
       s=tri(i)
c       print *,' on attaque le sommet ',s
c       print *,'++++++++++++++++++++++++++++++'
       pf=mshlcl(c,nu,tete,s,nbs)
c      creation d'un nouveau triangle et modification de la frontiere
c      --------------------------------------------------------------
       t=free
       free = nu(1,free)
       npf  = free
       free = nu(1,free)
       ppf  = nu(pp,pf)
       psf  = nu(ps,pf)
       ttaf  = nu(2,pf)
       taf  = ttaf / (2**3)
       iaf  = ttaf - (2**3) * taf
c
c                  npf 
c               1  x s               ---
c                 / \                ---
c              4 /   \ 6        ---  vide ---
c               /  t  \              ---
c            2 /   5   \ 3           ---
c------ --<---x---------x---------x- frontiere--<---
c          psf \  iaf  /  pf         ---
c               \ taf /         --- omega ---
c                \   /               ---
c                 \ /                ---
c                  x                 ---
c                                    ---
c     generation  de l'element t
       nu(1,t) = s
       nu(2,t) = nu(1,psf)
       nu(3,t) = nu(1,pf )
       nu(4,t) = -npf
       nu(5,t) = (2**3) * taf + iaf
       nu(6,t) = -pf
       nu(iaf,taf) = (2**3) * t + 5
c      update de la liste frontiere  
       nu(ps,npf) = psf
       nu(ps,pf ) = npf
       nu(pp,npf) = pf
       nu(pp,psf) = npf
       nu(1,npf)  = s
       nu(2,npf)  = (2**3) * t + 4
       nu(2,pf )  = (2**3) * t + 6
       if(iop.ge.50) call mshdrw (c,nu,6,t,iop)
       call mshopt (c,nu,t,5,nbs,iop,err)
       if (err.ne.0) then
         print*,'*** ERREUR mshopt'
         return
       endif
       call mshcvx (.true. ,c,nu,npf,nbs,iop,err)
       if (err.ne.0) then
         print*,'*** ERREUR mshcvx'
         return
       endif
       call mshcvx (.false.,c,nu,npf,nbs,iop,err)
       if (err.ne.0) then
         print*,'*** ERREUR mshcvx'
         return
       endif
30    continue
      end
c**********************************************************************
      subroutine mshopt (c,nu,t,a,nbs,iop,err)
      integer nbs,c(2,nbs),nu(6,nbs+nbs-2),t,a,iop,err
      integer vide
      parameter (vide=-2**30)
      integer mxpile
      parameter (mxpile=1024)
Cfj      parameter (mxpile=256)
      integer pile(2,mxpile)
      integer t1,t2,i,s1,s2,s3,s4,sin1,cos1,sin2,cos2,sgn
      integer tt1,tt,i11,i12,i13,i21,i22,i23,a1,a2,aa,mod3(1:3)
      real reel1,reel2
      real*8 reel8
      data mod3/2,3,1/
      i=1
      pile(1,i) = t 
      pile(2,i) = a
10    continue
      if(i.gt.0) then
        t1=pile(1,i)
        a1=pile(2,i)
        i=i-1
        if(t1.le.0) goto 10
        tt1 = nu(a1,t1)
        if(tt1.le.0) goto 10
        t2 = tt1/(2**3)
        a2 = tt1-t2*(2**3)
c        print *,' mshopt :t1,a1,t2,a2 =',t,a,nu(a,t)/8,mod(nu(a,t),8)
c     &         ,' niveau = ',i
        i11 =   a1 -3
        i12 =   mod3(i11) 
        i13 =   mod3(i12)
        i21 =   a2 -3
        i22 =   mod3(i21)
        i23 =   mod3(i22)
        s1 = nu(i13,t1)
        s2 = nu(i11,t1)
        s3 = nu(i12,t1)
        s4 = nu(i23,t2)
c        print *,i11,i12,i13,nu(i11,t1),nu(i12,t1),nu(i13,t1)
c        print *,i21,i22,i23,nu(i21,t2),nu(i22,t2),nu(i23,t2)
c        print *,s1,s2,s3,s4
c              critere d optimisation du quadrilatere
c----------------------------------------------------
        sin1 =   (c(2,s3)-c(2,s1)) * (c(1,s2)-c(1,s1))
     &         - (c(1,s3)-c(1,s1)) * (c(2,s2)-c(2,s1))
        cos1 =   (c(1,s3)-c(1,s1)) * (c(1,s3)-c(1,s2))
     &         + (c(2,s3)-c(2,s1)) * (c(2,s3)-c(2,s2))
        if(sin1.eq.0..and.cos1.eq.0.) then
c          print *,'fatal error mshopt:'
c     &           ,'3 points confondus ',s1,s2,s3
          err = 20
          return
        end if
c       b est la cotangente de angle (s1,s3,s2)
        sin2  =   (c(1,s4)-c(1,s1)) * (c(2,s2)-c(2,s1))
     &          - (c(2,s4)-c(2,s1)) * (c(1,s2)-c(1,s1))
        cos2  =   (c(1,s4)-c(1,s2)) * (c(1,s4)-c(1,s1))
     &          + (c(2,s4)-c(2,s2)) * (c(2,s4)-c(2,s1))
        reel1=float(cos2)*float(sin1)
        reel2=float(cos1)*float(sin2)
        if(abs(reel1)+abs(reel2).ge.2**30) then
c          print *,'on a un overflow en entier on calcule en reel*8'
          reel8  =    dble(cos2)*dble(sin1)
     &             +  dble(cos1)*dble(sin2)
          reel8=min(max(reel8,-1.d0),1.d0)
          sgn=reel8
        else
          sgn = cos2*sin1 + cos1*sin2
        endif
        if(min(max(sgn,-1),+1)*sin1.ge.0) goto 10
c       on inverse le quadrilatere
c       update des sommets
c-------------------------
        nu(i12,t1) = s4
        nu(i22,t2) = s1
c       update des aretes a1,a2
c-------------------------------
        tt1 = nu(i22+3,t2)
        nu(a1 ,t1) = tt1
        if(tt1.gt.0) then
         tt=tt1/(2**3)
         aa = tt1-(2**3)*tt
         nu(aa,tt)=  a1 +  (2**3) * t1
        elseif(tt1.ne.vide) then
         nu(2,-tt1)= a1 +  (2**3) * t1
        endif
        tt1 = nu(i12+3,t1)
        nu(a2 ,t2) = tt1
        if(tt1.gt.0) then
         tt=tt1/(2**3)
         aa=tt1-(2**3)*tt
         nu(aa,tt)= a2 +  (2**3) * t2
        elseif(tt1.ne.vide) then
         nu(2,-tt1)= a2 +  (2**3) * t2
        endif
        nu(i12+3,t1) =   i22+3 + (2**3)*t2
        nu(i22+3,t2) =   i12+3 + (2**3)*t1
        if(i+4.gt.mxpile) then
          print*,'*** fatal error mshopt la pile est trop petite '
     &         ,mxpile,i+4
          err =21
          return
        endif
        if(iop.ge.50) call mshdrw (c,nu,6,t1,iop)
        if(iop.ge.50) call mshdrw (c,nu,6,t2,iop)
        i=i+1
        pile(1,i)=t1
        pile(2,i)=a1
        i=i+1
        pile(1,i)=t2
        pile(2,i)=a2
        i=i+1
        pile(1,i)=t1
        pile(2,i)=i13+3
        i=i+1
        pile(1,i)=t2
        pile(2,i)=i23+3
        goto 10
      endif
      end
c**********************************************************************
      subroutine mshfrt (c,nu,nbs,arete,nba,sd,nbsd,reft,w,iop,err)
      integer nbs,c(2,nbs),nu(6,nbs+nbs-2),reft(2*(nbs-1))
      integer nba,nbsd,arete(2,nba),sd(2,nbsd),err,w(*),iop
      integer i,j,k,ifrt,nbt,is,ie,nbac,nbacpp,err1,itera
      integer ss1,s1,s2,t,ta,is1,s2t,s3t,det2,det3,a,ap,isd,jsd
      integer p3(1:3),impre
      integer vide
      parameter (vide=-2**30)
      logical fin,sens
      data p3/2,3,1/
      impre=mod(iop,10)
      if(nba.eq.0) return
c      ifrt=0
      nbt = nbs+nbs-2
      do 10 i=1,nbs
       reft(i)=0
10    continue
      do 20 i=1,nba
       reft(arete(1,i))=vide
       reft(arete(2,i))=vide
20    continue
      nbac = 0
      do 30 a=1,nba
       s1=min(arete(1,a),arete(2,a))
       s2=max(arete(1,a),arete(2,a))
       if(s1.eq.s2) then
c         print *,'warning :mshfrt l''arete ',a,' de sommets :',s1,s2
c     &         ,' est degeneree'
         nbac = nbac + 1
       else 
         i = reft(s1)
25       continue
         if(i.ne.vide) then
           if(s2.eq.max(arete(1,i),arete(2,i))) then
c              print *,'warning :mshfrt l''arete ',i,' = (',s1,s2 
c     &               ,')  est egale a arete ',a,' = ('
c     &               ,arete(1,a),arete(2,a),')'
             nbac = nbac + 1
           else
            i = w(i)
            goto 25
           endif
         else
           w(a)=reft(s1)
           reft(s1)=a
         endif
       endif
30    continue
      nbacpp = 1       
      itera = 0
      err1 = 0
50    continue
      itera = itera + 1
      if(err1.ne.0) then
        err=err1
        return
      endif
      if(nbac.lt.nba) then
        if(nbacpp.eq.0) then
c          print *,' fatal error mshfrt :l''algorithme boucle :'
c     &           ,nba,nbac,' iteration =',itera
c          if(impre.ge.9) then
c            print *,' dump '
c            print *,' nbt = ',nbt,' nbs =',nbs
c            print *,' tetes de listes ='
c            print '(6(i5,i9,'';''))',(i,reft(i),i=1,nbs)
c            print *,' chainage  = '
c            print '(4(i5,i6,i6,'';''))',(i,w(i),i=1,nba)
c            print *,' arete = '
c            print '(4(i5,i6,i6,'';''))'
c     &                       ,(i,arete(1,i),arete(2,i),i=1,nba)
c            print *,' i , nu(1:6,i) = '
cc            print '(10(i5,6(i12)/))',(i,(nu(j,i),j=1,6),i=1,nbt)
c          endif
          do 70 i=1,nbs
            a = reft(i)
60          continue
            if(a.gt.0) then
              s1 = arete(1,i)
              s2 = arete(2,i)
c              print *,' arete ',a,'  s1 = ',s1,' s2 = ',s2
              a=w(a)
              goto 60
            endif
70        continue
          err = 7
          return
        endif
c---------------------------------------------------------------------
c     on s'occupe des aretes a forcer
c---------------------------------------------------------------------
c        if(impre.ge.9) print *,'   on s''occupe des aretes a forcer'
        nbacpp = 0
        do 120 ie=1,nbt
         if(nu(5,ie).ne.0) then
           do 110 is=1,3
            s1  =nu(    is ,ie)
            s2t =nu( p3(is),ie)
            ss1 = min(s1,s2t)
            ap = 0
            a = reft(ss1)
80          continue
            if(a.gt.0) then
              s2 = max(arete(1,a),arete(2,a))
              t    = ie
              ta   = 0
              if(s2.eq.max(s1,s2t)) then
                if(iop.ge.70) call mshdrw(c,nu,6,ie,iop)
                if(nu(is+3,ie).gt.0) then
                  ta = nu(is+3,ie) /(2**3)
                  i  = nu(is+3,ie)-(2**3) * ta
                  nu(i,ta)=vide
                endif
                nu(is+3,ie)=vide
                goto 100
              endif
              ap = a
              a  = w(a)
              goto 80
            endif
            if(itera.eq.1) goto 110
            ss1 = s1
            ap = 0
            a = reft(ss1)
90          continue
            if(a.gt.0) then
              s2 = max(arete(1,a),arete(2,a))
              t    = ie
              ta   = 0
c             recherche si l' element coupe l''arete a
              is1  = is
              s3t  = nu(p3(p3(is)),t)
              det2 =  (c(1,s2t)-c(1,s1))*(c(2,s2)-c(2,s1))
     &              - (c(2,s2t)-c(2,s1))*(c(1,s2)-c(1,s1))
              det3 =  (c(1,s3t)-c(1,s1))*(c(2,s2)-c(2,s1))
     &              - (c(2,s3t)-c(2,s1))*(c(1,s2)-c(1,s1))
              if(impre.ge.9) then
c                print *,'t,is,det2,det3,s1,s2,s2t,s3t = '
c     &                 , t,is,det2,det3,s1,s2,s2t,s3t
              endif
              if(det2.gt.0.and.det3.lt.0) then
                call mshfr1 (c,nu,nbs,t,ta,is1,s2,iop,err)
                if (err.ne.0) then
                  print*,'*** ERREUR mshfr1'
                  return
                endif
                goto 100
              elseif(det2.eq.0.and.reft(s2t).eq.0) then
                print *,'*** fatal error mshfrt: le point ',s2t
     &                 ,' qui ne doit pas etre frontiere , l''est'
                err1 = 10
              elseif(det3.eq.0.and.reft(s3t).eq.0) then
                print *,'*** fatal error mshfrt: le point ',s3t
     &                 ,' qui ne doit pas etre frontiere , l''est'
                err1 = 10
              endif
              ap = a
              a  = w(a)
              goto 90
            endif
            goto 110
100         continue
            nbacpp = nbacpp + 1
            if(ap.eq.0) then
              reft(ss1)=w(a)
            else
              w(ap)=w(a)
            endif
            if(nbac+nbacpp.eq.nba) goto 130
110         continue
         endif
120      continue
        nbac = nbac + nbacpp
        goto 50
      endif
130   continue
c-----------------------------------------------------------------------
c     prise en compte des sous domaines 
c-----------------------------------------------------------------------
      do 140 i=1,nbs+nbsd+nbsd
       w(i)=0
140   continue
      do 150 i=1,nbsd
       a = abs(sd(1,i))
       s1=min(arete(1,a),arete(2,a))
       w(i+i)=w(s1+nbsd+nbsd)
       w(s1+nbsd+nbsd)=i
150   continue
      do 180 t=1,nbt  
       reft(t)=vide
       if(nu(6,t).ne.0) then
         do 170 i=1,3
           ss1=min(nu(i,t),nu(p3(i),t))
           jsd = nbsd+nbsd+ss1
160        continue           
           isd=w(jsd)
           if(isd.gt.0) then
             a=sd(1,isd)
             if(a.gt.0) then
               if(nu(i,t).eq.arete(1,a).and.nu(p3(i),t).eq.arete(2,a))
     &                                                 then
                 reft(t)=sd(2,isd)
                 w(isd+isd-1) = t
                 w(jsd)=w(isd+isd)
                 goto 170
               endif
             elseif(a.lt.0) then
               if(nu(i,t).eq.arete(2,-a).and.nu(p3(i),t).eq.arete(1,-a))
     &                                                 then
                 reft(t)=sd(2,isd)
                 w(isd+isd-1) = t
                 w(jsd)=w(isd+isd)
                 goto 170
               endif
             else
c               print *,' fatal erreur sous domaine ',isd,' arete nulle'
               err=11
             endif
             jsd = isd+isd
             goto 160
           endif
170      continue         
       endif
180   continue
      do 190 isd=1,nbsd
       if(w(isd+isd-1).eq.0 ) then
         err= 11
c         print4 *,' fatal erreur mshfrt le sous domaine ',isd
c     &          ,' n''est reference par aucun element '
c         print *,' revoir l''orientation '
c     &          ,' ou la definition des sous domaines'
       else
         w(isd+isd)=3
       endif
190   continue
      if (err.ne.0) then
        print *,
     &     '*** fatal error mshfrt :les sous domaines sont mal definis'
        print *,' l''arete ',a,' est reference plus d''une fois '
        print *,' nombre de sous domaine ',nbsd
        print '(a,/,10(/,i4,a,i10,a,i10))','sous domaine '
     &        ,(i,' arete= ',sd(1,i),' ref= ',sd(2,i),i=1,nbsd)
        return
      endif


      i=nbsd+nbsd
200   continue
      if(i.gt.0) then
        w(i)=w(i)+1
        if(w(i).le.6) then
          ta=nu(w(i),w(i-1))
          if(ta.gt.0) then
            ta = ta / (2**3)
            if(nu(1,ta).gt.0) then
              if(iop.ge.80) call mshdrw(c,nu,6,ta,iop)
              nu(1,ta)=-nu(1,ta)
              if(reft(ta).ne.reft(w(i-1))) then
                if(reft(ta).ne.vide) then
c                  print *,' mshfrt :error sous domaine element ',ta
c     &               ,' ref old = ',reft(ta),' ref new = ',reft(w(i-1))
                else
                  reft(ta)=reft(w(i-1))
                endif
                w(i+1)=ta
                w(i+2)=3
                i=i+2
              endif
            endif
          endif
        else
          i=i-2
        endif
        goto  200
      endif
      do 220 ie=1,nbt
       if(nu(1,ie).lt.0) then
         nu(1,ie)=-nu(1,ie)
       else
         do 210 i=1,6
          nu(i,ie)=0
210      continue
       endif
220   continue
      end
c**********************************************************************
      subroutine mshfr1 (c,nu,nbs,it1,ita,is1,s2,iop,err)
      integer nbs,c(2,nbs),nu(6,nbs+nbs-2),it,is1,s2,err,it1,ita,iop
      integer lstmx
      parameter (lstmx=256)
      integer lst(3,lstmx)
      integer i,j,k,s1,s,s3,x,y,det,det2,det3,nbac,s2t,s3t,a,t,ta
      integer l1,l2,l3,la,p3(1:5)
      data p3 /2,3,1,2,3/
      t = it1
      s1 = nu(is1,t)
      x = c(1,s2)-c(1,s1)
      y = c(2,s2)-c(2,s1)
      nbac = 0
      l1 = is1
      l2 = p3(l1)
      l3 = p3(l2)
c      s2t = nu(l2,t)
c      s3t = nu(l3,t)
      la = l2 + 3
c      print *,'  mshfr1 :',it1,is1,s1,s2
20    continue
      nbac = nbac + 1
      if(nbac.gt.lstmx) then
c        print *,' fatal error mshfr1 : lst trop petit ',nbac,lstmx
        err =8
        return
      endif
      lst(2,nbac) = t
      lst(3,nbac) = la
      if(iop.ge.70) call mshdrw (c,nu,6,t,iop)
      ta = nu(la,t)
      if(ta.le.0) then
c        print *,' fatal error mshfr1:la frontiere est croisee en ',t
        err =9
        return
      endif
      t  = ta/8
      la = ta-8*t 
      s3 = nu(p3(la-2),t)
      if(s3.ne.s2) then
        det = x*(c(2,s3)-c(2,s1))-y*(c(1,s3)-c(1,s1))
c        print *,' s3 = ',s3,det
        if(det.gt.0) then
          la = 3+p3(la-3)
        elseif(det.lt.0) then
          la = 3+p3(la-2)
        else
c          print2 *,' fatal error mshfr1: le point ',s3
c     &           ,' qui ne doit pas etre frontiere , l''est'
          err = 10
          return
        endif
        goto 20
      endif
c     c'est la que l'on peut construire la nouvelle frontiere
c     avec lst,nbac
      if(iop.ge.70) call mshdrw (c,nu,6,t,iop)
      call mshfr2 (c,nu,nbs,lst,nbac,it1,ita,s1,s2,iop,err)
      return
      end
c**********************************************************************
      subroutine mshfr2 (c,nu,nbs,lst,nbac,t,ta,ss1,ss2,iop,err)
      integer nbs,nbac,c(2,nbs),nu(6,nbs+nbs-2),lst(3,nbac)
      integer t,ta,ss1,ss2,iop,err
      integer ptlst,ttlst,pslst,pplst,s1,s2,s3,s4,x41,y41,x,y
      integer i,t1,a1,tt1,t2,a2,tt,i11,i12,i13,i21,i22,i23,aas,aa
      integer det1,det4,det2,det3
      integer mod3(3)
      integer vide
      parameter (vide=-2**30)
      data mod3/2,3,1/
      x = c(1,ss1)-c(1,ss2)
      y = c(2,ss1)-c(2,ss2)
      do 10 i=1,nbac-1
        lst(1,i)=i+1
10    continue
      lst(1,nbac)=0
      ttlst = 1
20    continue
      ptlst  = ttlst
      pplst  = 0
30    continue
      if(ptlst.gt.0) then
        t1=lst(2,ptlst)
        a1=lst(3,ptlst)
        tt1 = nu(a1,t1)
        t2 = tt1/(2**3)
        a2 = tt1-t2*(2**3)
        i11 =   a1 -3
        i12 =   mod3(i11) 
        i13 =   mod3(i12)
        i21 =   a2 -3
        i22 =   mod3(i21)
        i23 =   mod3(i22)
        s1 = nu(i13,t1)
        s2 = nu(i11,t1)
        s3 = nu(i12,t1)
        s4 = nu(i23,t2)
        x41 = c(1,s4)-c(1,s1)
        y41 = c(2,s4)-c(2,s1)
        det2 = (c(1,s2)-c(1,s1))*y41-(c(2,s2)-c(2,s1))*x41
        det3 = (c(1,s3)-c(1,s1))*y41-(c(2,s3)-c(2,s1))*x41
        if(det2.gt.0.and.det3.lt.0) then
c         le quadrilataire est convexe on le retourne
c         update des sommets
c-------------------------
          nu(i12,t1) = s4
          nu(i22,t2) = s1
c         update du pointeur suivant
c-----------------------------------
          pslst=lst(1,ptlst)
          if(pslst.gt.0) then
            aas=lst(3,pslst)
            if(aas.eq.i22+3) then
              lst(2,pslst) = t1
              lst(3,pslst) = i11 + 3
            endif
          endif
c         update des aretes a1,a2
c-------------------------------
          tt1 = nu(i22+3,t2)
          nu(a1 ,t1) = tt1
          if(tt1.gt.0) then
            tt=tt1/(2**3)
            aa = tt1-(2**3)*tt
            nu(aa,tt)=  a1 +  (2**3) * t1
          elseif(tt1.ne.vide) then
            nu(2,-tt1)= a1 +  (2**3) * t1
          endif
          tt1 = nu(i12+3,t1)
          nu(a2 ,t2) = tt1
          if(tt1.gt.0) then
            tt=tt1/(2**3)
            aa=tt1-(2**3)*tt
            nu(aa,tt)= a2 +  (2**3) * t2
          elseif(tt1.ne.vide) then
            nu(2,-tt1)= a2 +  (2**3) * t2
          endif
          nu(i12+3,t1) =   i22+3 + (2**3)*t2
          nu(i22+3,t2) =   i12+3 + (2**3)*t1
          det1 = (c(1,s1)-c(1,ss1))*y-(c(2,s1)-c(2,ss1))*x
          det4 = (c(1,s4)-c(1,ss1))*y-(c(2,s4)-c(2,ss1))*x
          if(iop.ge.50) call mshdrw (c,nu,6,t1,iop)
          if(iop.ge.50) call mshdrw (c,nu,6,t2,iop)
          if(det1.lt.0.and.det4.gt.0) then
c           le sommets s4 est dans omega
            lst(2,ptlst) = t2
            lst(3,ptlst) = i22+3
          elseif(det1.gt.0.and.det4.lt.0) then
c           le sommets s1 est dans omega
            lst(2,ptlst) = t1
            lst(3,ptlst) = i12+3
          else
c           print *,'    on supprime l''arete dans  lst ',t1,a1,t2,a2
            if(pplst.eq.0) then
              ttlst = lst(1,ptlst)
              ptlst = ttlst
            else
              ptlst        = lst(1,ptlst)
              lst(1,pplst) = ptlst
            endif
            goto 30
          endif
        endif
        pplst = ptlst
        ptlst = lst(1,ptlst)
        goto 30
      endif
      if(ttlst.ne.0) goto 20
      nu(i12+3,t1) =  vide
      nu(i22+3,t2) =  vide
      t  = t2      
      ta = t1
      do 40 i=1,nbac
       call mshopt (c,nu,lst(2,i),4,nbs,iop,err)
       call mshopt (c,nu,lst(2,i),5,nbs,iop,err)
       call mshopt (c,nu,lst(2,i),6,nbs,iop,err)
40    continue
      end
