[interpolation] mesh2d gateway introduced
[scilab.git] / scilab / modules / interpolation / src / fortran / mesh2b.f
1 c Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
2 c Copyright (C) INRIA
3 c
4 c This file must be used under the terms of the CeCILL.
5 c This source file is licensed as described in the file COPYING, which
6 c you should have received as part of this distribution.  The terms
7 c are also available at
8 c http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
9
10
11       subroutine mesh2b(nbs,n6,n4,lfront,cr,c,nu,tri,front,nbt,err)
12 c-----------------------------------------------------------------------
13 c triangulation  2d a partir d'un ensemble de points et d'une frontiere
14 c de domaine decrite par ses composantes connexes
15 c-----------------------------------------------------------------------
16 c entrees :
17 c     cr     tableau des coordonnees des points c(1:2,1:nbs)
18 c     nbs    nombre de points
19 c      n6=6*(nbs+nbs-2)       n4=4*nbs-4
20 c     front  tableau definisant la frontiere par composantes connexes de la
21 c            frontiere (i1,i2) front(i1)=front(i2) front(i+1) est le
22 c            point frontiere suivant de front(i) pour i=i1,i2-1 telle que
23 c            la normale a la frontiere soit interne ( la  composante
24 c            connexe exterieure tourne dans le sens trigo )
25 c     lfront est la longueur du tableau front
26 c----------------------------------------------------------------------
27 c tableaux de travail
28 c     c(2,nbs) tableau d'entiers (copie de coordonnees)
29 c     tri(4*nbs-4) tableau d'entiers
30 c     nu (6*(2*nbs-2)) tableau d entiers contiendra le tableau
31 c             des connectivites ( au debut) a la sortie du sp
32 c
33 c sorties :
34 c     nbt         nombre de triangles generes
35 c     nu(1:3,nbt) sommets des triangles (tableau des connections)
36 c                        (!! nu a dimensionner en entree nu(1:6*(2*nbs-2))
37 c     err         si err = 0 alors pas de probleme  (ok)
38 c                        sinon nbt = 0 et pas de triangulation
39 c     c(1:2,nbs)  coordonnees des sommets (en entier) cr n est pas modifie
40 c---------------------------------------------------------------------------
41 c hecht-marrocco  inria-rocquencourt  (39 63 55 14)
42 c-------------------------------------------------------------------------
43       integer nbs,nbt,lfront,n6,n4,err
44       integer c(2,nbs),tri(n4),front(lfront),nu(n6)
45       double precision cr(2,nbs)
46       integer i,j,k,tete,i1,i2
47 c
48       err = 0
49       nbt = 0
50       do 1,i=1,nbs
51          c(1,i)=0
52          c(2,i)=0
53  1    continue
54       do 2,i=1,n6
55          nu(i)=0
56  2    continue
57 c preparation des donnees
58       call mshtri (cr,c,nbs,tri,tri(nbs+1),err)
59       if(err.ne.0) return
60 c maillage de l enveloppe convexe
61       call mshcxi (c,nu,tri,nbs,tete,err,n4)
62       if(err.ne.0) return
63 c
64       do 10 i=1,nbs
65       tri(i)=0
66 10    continue
67       i=tete
68 20    continue
69       j=nu(6*(i-1)+4)
70       tri(nu(6*(i-1)+1))=nu(6*(j-1)+1)
71       i=j
72       if(i.ne.tete) goto 20
73 c traitement frontiere
74       k=0
75       if(lfront.gt.0) then
76       call mshfrt(c,nu,nbs,front,lfront,tri,err,n4)
77       if(err.ne.0) return
78       do 30 i=1,nbs
79       tri(i)=0
80 30    continue
81       k  = 0
82       i2 = 0
83       do 40 i=1,lfront
84       i1=i2
85       i2=front(i)
86       if(i1.eq.k) then
87       k = -i2
88       elseif(i1.eq.-k) then
89       k=-k
90       tri(i1)=i2
91       else
92       tri(i1)=i2
93       endif
94 40    continue
95       endif
96 c     construction du tableau nu(1:3,1:nbt)
97       nbt=0
98       k = 0
99       do 200 j=1,6*(nbs+nbs-2),6
100       if(nu(j+5).ne.0) then
101       nbt=nbt + 1
102       do 190 i=0,2
103       k=k+1
104       nu(k)=nu(j+i)
105 190   continue
106       endif
107 200   continue
108       end
109 c**********************************************************************
110       integer function mshlcl(c,nu,tete,s,nbs)
111       integer nbs,c(2,nbs),nu(6,nbs+nbs-2),tete,s
112       integer x,y,pt,ppt,det
113       logical init
114 c
115       x=c(1,s)
116       y=c(2,s)
117       init=.true.
118       pt=tete
119 10    continue
120       ppt=pt
121       pt=nu(4,pt)
122       if(pt.ne.tete) then
123       det=x*c(2,nu(1,pt)) -y*c(1,nu(1,pt))
124       if(det.lt.0) then
125       init=.false.
126       goto 10
127       elseif(init.and.det.eq.0) then
128       goto 10
129       endif
130       endif
131       mshlcl=ppt
132       end
133 c**********************************************************************
134       subroutine mshtri (cr,c,nbs,tri,nu,err)
135       integer nbs,c(1:2,1:nbs),tri(1:nbs),nu(1:nbs),err
136       double precision cr(1:2,1:nbs)
137       integer iii,ic,xx,ip,i,j,jc,k,trik,tri3,det,ierr
138       double precision aa1,aa2,xmin,xmax,ymin,ymax
139       double precision precis
140       parameter (precis=2.**15-1.)
141 c
142       err = 0
143       ierr = 0
144       iii=1
145       xmin=cr(1,1)
146       ymin=cr(2,1)
147       xmax=cr(1,1)
148       ymax=cr(2,1)
149       do 10 ic=1,nbs
150        xmin=min(cr(1,ic),xmin)
151        ymin=min(cr(2,ic),ymin)
152        xmax=max(cr(1,ic),xmax)
153        ymax=max(cr(2,ic),ymax)
154        tri(ic)=ic
155        if(cr(1,ic).lt.cr(1,iii)) then
156         iii=ic
157        endif
158 10    continue
159       aa1 = precis/(xmax-xmin)
160       aa2 = precis/(ymax-ymin)
161       aa1 = min(aa1,aa2)
162       aa2 = aa1*(cr(2,iii)-ymin)
163       do 20 ic=1,nbs
164        c(1,ic) = nint(aa1*(cr(1,ic)-cr(1,iii)))
165        c(2,ic) = nint(aa1*(cr(2,ic)-ymin)-aa2)
166        nu(ic)= c(1,ic)**2 + c(2,ic)**2
167 20    continue
168 c----------------------------------------------------------
169       call mshtr1 (nu,tri,nbs)
170       ip = 1
171       xx=nu(ip)
172       do 30 jc=1,nbs
173       if(nu(jc).gt.xx)then
174       call mshtr1 (nu(ip),tri(ip),jc-ip)
175       do 25 i=ip,jc-2
176       if(nu(i).eq.nu(i+1)) then
177       ierr=ierr+1
178 c      print *,' error les points ',tri(i),tri(i+1),' sont egaux'
179       endif
180 25    continue
181       xx=nu(jc)
182       ip=jc
183       endif
184       ic=tri(jc)
185       nu(jc)=c(2,ic)
186 30    continue
187       call mshtr1 (nu(ip),tri(ip),nbs-ip)
188       do 35 i=ip,jc-2
189       if(nu(i).eq.nu(i+1)) then
190       ierr=ierr+1
191 c      print *,' error les points ',tri(i),tri(i+1),' sont egaux'
192       endif
193 35    continue
194       if(ierr.ne.0) then
195       err = 2
196 c      print *,' fatal error mshtri:il y a des points confondus'
197       return
198       endif
199       k=2
200 50    continue
201       if(k.le.nbs) then
202       k=k+1
203       det = c(1,tri(2))*c(2,tri(k)) - c(2,tri(2))*c(1,tri(k))
204       if(det.eq.0) goto 50
205       else
206 c      print *,'fatal error mshtri tous les points sont alignes'
207 c      print *,'tri =',(tri(k),k=1,nbs)
208       err = 3
209       return
210       endif
211       trik = tri(k)
212       do 60 j=k-1,3,-1
213       tri(j+1)=tri(j)
214 60    continue
215       tri(3)=trik
216       if(det.lt.0) then
217       tri3=tri(3)
218       tri(3)=tri(2)
219       tri(2)=tri3
220       endif
221       end
222 c**********************************************************************
223       subroutine mshtr1 (criter,record,n)
224       integer record(n)
225       integer criter(n)
226 c
227       integer i,l,r,j,n
228       integer rec
229       integer crit
230 c
231       if(n.le.1) return
232       l=n/2+1
233       r=n
234 2     if(l.le.1)goto 20
235       l=l-1
236       rec=record(l)
237       crit=criter(l)
238       goto 3
239 20    continue
240       rec=record(r)
241       crit=criter(r)
242       record(r)=record(1)
243       criter(r)=criter(1)
244       r=r-1
245       if(r.eq.1)goto 999
246 3     j=l
247 4     i=j
248       j=2*j
249       CRES=j-r
250       if (CRES .lt. 0) then
251          goto 5
252       elseif (CRES .eq. 0) then
253          goto 6
254       else
255          goto 8
256       endif
257 5     if(criter(j).lt.criter(j+1))j=j+1
258 6     if(crit.ge.criter(j))goto 8
259       record(i)=record(j)
260       criter(i)=criter(j)
261       goto 4
262 8     record(i)=rec
263       criter(i)=crit
264       goto 2
265 999   record(1)=rec
266       criter(1)=crit
267       return
268       end
269 c**********************************************************************
270       subroutine mshcvx(direct,c,nu,pfold,nbs,err)
271       integer nbs,c(2,nbs),nu(6,nbs+nbs-2),pfold,err
272       logical direct
273       integer pp,ps,i1,i2,i3,i4,i5,i6
274       integer pf,psf,ppf,s1,s2,s3,t,t4,t5,a4,a5,det,tt4,tt5
275       if(direct) then
276       pp=3
277       ps=4
278       i1=1
279       i2=3
280       i3=2
281       i4=6
282       i5=5
283       i6=4
284       else
285       pp=4
286       ps=3
287       i1=1
288       i2=2
289       i3=3
290       i4=4
291       i5=5
292       i6=6
293       endif
294 10    continue
295       ppf=pfold
296       pf =nu(ps,pfold)
297       psf=nu(ps,pf)
298       s1=nu(1,ppf)
299       s2=nu(1,pf)
300       s3=nu(1,psf)
301       det =   ( c(1,s2) - c(1,s1) ) * ( c(2,s3) - c(2,s1) )
302      &      - ( c(2,s2) - c(2,s1) ) * ( c(1,s3) - c(1,s1) )
303       if(((.not.direct).and.det.gt.0).or.(direct.and.det.lt.0)) then
304       if(direct) then
305       tt4 = nu(2,ppf)
306       tt5 = nu(2,pf)
307       else
308       tt4 = nu(2,pf)
309       tt5 = nu(2,psf)
310       endif
311       t4 = tt4/(2**3)
312       t5 = tt5/(2**3)
313       a4 = tt4 -8 * t4
314       a5 = tt5 -8 * t5
315       nu(ps,ppf) = psf
316       nu(pp,psf) = ppf
317       t = pf
318       if(direct) then
319       nu(2,ppf) = (2**3) * t + i6
320       else
321       nu(2,psf) = (2**3) * t + i6
322       endif
323       nu(i1,t ) = s1
324       nu(i2,t ) = s2
325       nu(i3,t ) = s3
326       nu(i4,t ) = (2**3) * t4 + a4
327       nu(i5,t ) = (2**3) * t5 + a5
328       if(direct) then
329       nu(i6,t ) = -ppf
330       else
331       nu(i6,t ) = -psf
332       endif
333       nu(a4,t4) = (2**3) * t + i4
334       nu(a5,t5) = (2**3) * t + i5
335       call mshopt (c,nu,t5,a5,nbs,err)
336       if(err.ne.0) return
337       goto 10
338       endif
339       end
340 c**********************************************************************
341       subroutine mshcxi (c,nu,tri,nbs,tete,err,n4)
342       integer nbs,c(2,nbs),nu(6,2*nbs-2),tri(n4),tete
343       integer mshlcl,err,n4
344       integer i,j,s,t,pf,ppf,psf,npf,pp,ps,taf,iaf,free,ttaf
345       parameter (pp=3,ps=4)
346       do 10 i=1,nbs+nbs-2
347       nu(1,i)=i+1
348       do 10 j=2,6
349       nu(j,i)=0
350 10    continue
351       nu(1,nbs+nbs-2)=0
352       free = 1
353       t=free
354       free = nu(1,free)
355       tete=free
356       pf  =free
357       do 20 i=1,3
358       nu(i  ,t) = tri(i)
359       nu(3+i,t) = -pf
360       ppf       = pf
361       free      = nu(1,pf)
362       pf        = free
363       if(i.eq.3) pf=tete
364       nu(1,ppf) = tri(i)
365       nu(2,ppf) = i + 3 + (2**3) * t
366       nu(ps,ppf) = pf
367       nu(pp,pf ) = ppf
368 20    continue
369       do 30 i=4,nbs
370       s=tri(i)
371       pf=mshlcl(c,nu,tete,s,nbs)
372       t=free
373       free = nu(1,free)
374       npf  = free
375       free = nu(1,free)
376       ppf  = nu(pp,pf)
377       psf  = nu(ps,pf)
378       ttaf  = nu(2,pf)
379       taf  = ttaf / (2**3)
380       iaf  = ttaf - (2**3) * taf
381       nu(1,t) = s
382       nu(2,t) = nu(1,psf)
383       nu(3,t) = nu(1,pf )
384       nu(4,t) = -npf
385       nu(5,t) = (2**3) * taf + iaf
386       nu(6,t) = -pf
387       nu(iaf,taf) = (2**3) * t + 5
388       nu(ps,npf) = psf
389       nu(ps,pf ) = npf
390       nu(pp,npf) = pf
391       nu(pp,psf) = npf
392       nu(1,npf)  = s
393       nu(2,npf)  = (2**3) * t + 4
394       nu(2,pf )  = (2**3) * t + 6
395       call mshopt (c,nu,t,5,nbs,err)
396       if(err.ne.0) return
397       call mshcvx (.true. ,c,nu,npf,nbs,err)
398       if(err.ne.0) return
399       call mshcvx (.false.,c,nu,npf,nbs,err)
400       if(err.ne.0) return
401 30    continue
402       end
403 c**********************************************************************
404       subroutine mshopt (c,nu,t,a,nbs,err)
405       integer nbs,c(2,nbs),nu(6,nbs+nbs-2),t,a,err
406       integer vide
407       parameter (vide=-2**30)
408       integer mxpile
409       parameter (mxpile=512)
410       integer pile(2,mxpile)
411       integer t1,t2,i,s1,s2,s3,s4,sin1,cos1,sin2,cos2,sgn
412       integer tt1,tt,i11,i12,i13,i21,i22,i23,a1,a2,aa,mod3(1:3)
413       double precision reel1,reel2
414       double precision reel8
415       data mod3/2,3,1/
416       i=1
417       pile(1,i) = t 
418       pile(2,i) = a
419 10    continue
420       if(i.gt.0) then
421       t1=pile(1,i)
422       a1=pile(2,i)
423       i=i-1
424       if(t1.le.0) goto 10
425       tt1 = nu(a1,t1)
426       if(tt1.le.0) goto 10
427       t2 = tt1/(2**3)
428       a2 = tt1-t2*(2**3)
429       i11 =   a1 -3
430       i12 =   mod3(i11) 
431       i13 =   mod3(i12)
432       i21 =   a2 -3
433       i22 =   mod3(i21)
434       i23 =   mod3(i22)
435       s1 = nu(i13,t1)
436       s2 = nu(i11,t1)
437       s3 = nu(i12,t1)
438       s4 = nu(i23,t2)
439         sin1 =   (c(2,s3)-c(2,s1)) * (c(1,s2)-c(1,s1))
440      &         - (c(1,s3)-c(1,s1)) * (c(2,s2)-c(2,s1))
441         cos1 =   (c(1,s3)-c(1,s1)) * (c(1,s3)-c(1,s2))
442      &         + (c(2,s3)-c(2,s1)) * (c(2,s3)-c(2,s2))
443         if(sin1.eq.0.and.cos1.eq.0) then
444 c          print *,'fatal error mshopt:'
445 c     &           ,'3 points confondus ',s1,s2,s3
446           err = 12
447           return
448         end if
449         sin2  =   (c(1,s4)-c(1,s1)) * (c(2,s2)-c(2,s1))
450      &          - (c(2,s4)-c(2,s1)) * (c(1,s2)-c(1,s1))
451         cos2  =   (c(1,s4)-c(1,s2)) * (c(1,s4)-c(1,s1))
452      &          + (c(2,s4)-c(2,s2)) * (c(2,s4)-c(2,s1))
453         reel1=float(cos2)*float(sin1)
454         reel2=float(cos1)*float(sin2)
455         if(abs(reel1)+abs(reel2).ge.2**30) then
456         reel8=dble(cos2)*dble(sin1)
457      &       +dble(cos1)*dble(sin2)
458         reel8=min(max(reel8,-1.d0),1.d0)
459         sgn=reel8
460         else
461         sgn = cos2*sin1 + cos1*sin2
462         endif
463         if(min(max(sgn,-1),+1)*sin1.ge.0) goto 10
464         nu(i12,t1) = s4
465         nu(i22,t2) = s1
466         tt1 = nu(i22+3,t2)
467         nu(a1 ,t1) = tt1
468         if(tt1.gt.0) then
469          tt=tt1/(2**3)
470          aa = tt1-(2**3)*tt
471          nu(aa,tt)=  a1 +  (2**3) * t1
472         elseif(tt1.ne.vide) then
473          nu(2,-tt1)= a1 +  (2**3) * t1
474         endif
475         tt1 = nu(i12+3,t1)
476         nu(a2 ,t2) = tt1
477         if(tt1.gt.0) then
478          tt=tt1/(2**3)
479          aa=tt1-(2**3)*tt
480          nu(aa,tt)= a2 +  (2**3) * t2
481         elseif(tt1.ne.vide) then
482          nu(2,-tt1)= a2 +  (2**3) * t2
483         endif
484         nu(i12+3,t1) =   i22+3 + (2**3)*t2
485         nu(i22+3,t2) =   i12+3 + (2**3)*t1
486         if(i+4.gt.mxpile) then
487 c          print *,' fatal error mshopt la pile est trop petite ',mxpile
488           err =13
489           return
490         endif
491 c
492         i=i+1
493         pile(1,i)=t1
494         pile(2,i)=a1
495         i=i+1
496         pile(1,i)=t2
497         pile(2,i)=a2
498         i=i+1
499         pile(1,i)=t1
500         pile(2,i)=i13+3
501         i=i+1
502         pile(1,i)=t2
503         pile(2,i)=i23+3
504         goto 10
505       endif
506       end
507 c**********************************************************************
508       subroutine mshfrt (c,nu,nbs,frt,lfrt,w,err,n4)
509       integer nbs,c(2,nbs),nu(6,nbs+nbs-2)
510       integer lfrt,frt(lfrt),err,w(n4)
511       integer i,ifrt,sinit,lnu,is,ie,tinter,nbac,nbaf,nbacpp
512       integer s0,s1,s2,ta,is1,it,s2t,s3t,det2,det3
513       integer p3(1:3)
514       integer vide
515       parameter (vide=-2**30)
516       logical fin
517       data p3/2,3,1/
518       if(lfrt.eq.0) return
519       tinter =0
520       ifrt=0
521       lnu = nbs+nbs-2
522 c       inite du tableau w
523       do 10 i=1,nbs
524        w(i)=-1
525 10    continue
526       nbaf = 0
527       s1 = 0
528       sinit= 0
529       fin =.true.
530       do 20 i=1,lfrt
531       s0 = s1
532       s1 = frt(i)
533       if(s1.le.0.or.s1.gt.nbs) then
534       err=5
535 c      print *,' fatal error mshfrt '
536 c      print *,' le tableau des la frontiere est mauvais en ',i,s1
537       return
538       endif
539       if(s0.eq.sinit) then
540       if(fin) then
541       sinit=s1
542       else
543       nbaf = nbaf + 1
544       if(w(s0).ne.-1) then
545 c      print *,'fatal error mshfrt : la frontiere est croisee '
546 c     &       ,' en ',s0                    
547       err=6
548       endif
549       w(s0)=i
550       endif
551       fin=.not.fin
552       else
553       nbaf = nbaf + 1
554       if(w(s0).ne.-1) then
555 c      print *,'fatal error mshfrt : la frontiere est croisee '
556 c     &       ,' en ',s0                    
557       err=6
558       endif
559       w(s0)=i
560       endif
561 20    continue
562       if(sinit.ne.s1) then
563 c      print *,'warning mshfrt:la frontiere n''est pas fermee'
564 c     &         ,' on la ferme avec l''arete ',s1,sinit
565       if(w(s1).ne.-1) then
566 c      print *,'fatal error mshfrt : la frontiere est croisee '
567 c     &           ,' en ',s1
568       err=6
569       endif
570       w(s1)=sinit
571       nbaf = nbaf + 1
572       endif
573       nbac = 0
574       nbacpp = 1
575 30    continue
576       if(err.ne.0) return
577       if(nbac.lt.nbaf) then
578       if(nbacpp.eq.0) then
579       err = 7
580 c      print *,' fatal error mshfrt :l''algorithme boucle :'
581 c     &           ,nbaf,nbac
582 c      print *,' la frontiere est certainement mal orientee '
583       return
584       endif
585 c     on s'occupe des aretes du maillage et frontiere de omega
586 c---------------------------------------------------------------------
587       nbacpp = 0
588       do 60 ie=1,lnu
589       if(nu(5,ie).ne.0) then
590       do 50 is=1,3
591       s1  =nu(    is ,ie)
592       s2t =nu( p3(is),ie)
593       if(w(s1).gt.0) then
594       s2   = frt(w(s1))
595       if(s2.eq.s2t) then
596       tinter = ie
597       nbacpp = nbacpp + 1 
598       w(s1) = 0
599       if(nu(is+3,ie).gt.0) then
600       ta = nu(is+3,ie) /(2**3)
601       i  = nu(is+3,ie)-(2**3) * ta
602       nu(i,ta)=vide
603       endif
604       nu(is+3,ie)=vide
605       else
606       it   = ie
607       is1  = is
608       s3t  = nu(p3(p3(is)),it)
609       det2 =  (c(1,s2t)-c(1,s1))*(c(2,s2)-c(2,s1))
610      &      - (c(2,s2t)-c(2,s1))*(c(1,s2)-c(1,s1))
611       det3 =  (c(1,s3t)-c(1,s1))*(c(2,s2)-c(2,s1))
612      &      - (c(2,s3t)-c(2,s1))*(c(1,s2)-c(1,s1))
613       if(det2.ge.0.and.det3.le.0) then
614       if(det2.eq.0) then
615       if(w(s2t).eq.-1) then
616       print *,' fatal error mshfrt: le point ',s2t
617      &       ,' qui ne doit pas etre frontiere , l''est'
618       err = 10
619       endif
620       goto 50
621       endif
622       if(det3.eq.0) then
623       if(w(s3t).eq.-1) then
624       print *,' fatal error mshfrt: le point ',s3t
625      &      ,' qui ne doit pas etre frontiere , l''est'
626       err = 10
627       endif
628       goto 50
629       endif
630       call mshfr1 (c,nu,nbs,it,is1,s2,err)
631       if(err.ne.0) return
632       tinter=it
633       w(s1) = 0
634       nbacpp = nbacpp + 1
635       endif
636       endif
637       endif
638 50    continue
639       endif
640 60    continue
641       nbac = nbac + nbacpp
642       goto 30
643       endif
644       i=2
645       w(1)=tinter
646       w(2)=3
647       nu(1,tinter) = -nu(1,tinter)
648 70    continue
649       if(i.gt.0) then
650       w(i)=w(i)+1
651       if(w(i).le.6) then
652       ta=nu(w(i),w(i-1))
653       if(ta.gt.0) then
654       ta = ta / (2**3)
655       if(nu(1,ta).gt.0) then
656       w(i+1)=ta
657       w(i+2)=3
658       i=i+2
659       nu(1,ta)=-nu(1,ta)
660       endif
661       endif
662       else
663       i=i-2
664       endif
665       goto  70
666       endif
667       do 90 ie=1,lnu
668       if(nu(1,ie).lt.0) then
669       nu(1,ie)=-nu(1,ie)
670       else
671       do 80 i=1,6
672       nu(i,ie)=0
673 80    continue
674       endif
675 90    continue
676       end
677 c**********************************************************************
678       subroutine mshfr1 (c,nu,nbs,it1,is1,s2,err)
679       integer nbs,c(2,nbs),nu(6,nbs+nbs-2),is1,s2,err,it1
680       integer lstmx
681       parameter (lstmx=256)
682       integer lst(3,lstmx)
683       integer s1,s3,x,y,det,nbac,s2t,s3t,t,ta
684       integer l1,l2,l3,la,p3(1:5)
685       logical direct
686       data p3 /2,3,1,2,3/
687       direct = .true.
688       t = it1
689       s1 = nu(is1,t)
690       x = c(1,s2)-c(1,s1)
691       y = c(2,s2)-c(2,s1)
692       nbac = 0
693       l1 = is1
694       l2 = p3(l1)
695       l3 = p3(l2)
696       s2t = nu(l2,t)
697       s3t = nu(l3,t)
698       la = l2 + 3
699 20    continue
700       nbac = nbac + 1
701       if(nbac.gt.lstmx) then
702         print *,' fatal error mshfr1 : lst trop petit ',nbac,lstmx
703         err =8
704         return
705       endif
706       lst(2,nbac) = t
707       lst(3,nbac) = la
708       ta = nu(la,t)
709       if(ta.le.0) then
710         print *,' fatal error mshfr1:la frontiere est croisee en ',t
711         err =9
712         return
713       endif
714       t  = ta/8
715       la = ta-8*t 
716       s3 = nu(p3(la-2),t)
717       if(s3.ne.s2) then
718         det = x*(c(2,s3)-c(2,s1))-y*(c(1,s3)-c(1,s1))
719         if(det.gt.0) then
720           la = 3+p3(la-3)
721         elseif(det.lt.0) then
722           la = 3+p3(la-2)
723         else
724           print *,' fatal error mshfr1: le point ',s3
725      &           ,' qui ne doit pas etre frontiere , l''est'
726           err = 10
727           return
728         endif
729         goto 20
730       endif
731       call mshfr2 (c,nu,nbs,lst,nbac,it1,s1,s2)
732       return
733       end
734 c**********************************************************************
735       subroutine mshfr2 (c,nu,nbs,lst,nbac,t,ss1,ss2)
736       integer nbs,nbac,c(2,nbs),nu(6,nbs+nbs-2),lst(3,nbac)
737       integer t,ss1,ss2
738       integer ptlst,ttlst,pslst,pplst,s1,s2,s3,s4,x41,y41,x,y
739       integer i,t1,a1,tt1,t2,a2,tt,i11,i12,i13,i21,i22,i23,aas,aa
740       integer det1,det4,det2,det3
741       integer mod3(3)
742       integer vide
743       parameter (vide=-2**30)
744       data mod3/2,3,1/
745       x = c(1,ss1)-c(1,ss2)
746       y = c(2,ss1)-c(2,ss2)
747       do 10 i=1,nbac-1
748       lst(1,i)=i+1
749 10    continue
750       lst(1,nbac)=0
751       ttlst = 1
752 20    continue
753       ptlst  = ttlst
754       pplst  = 0
755 30    continue
756       if(ptlst.gt.0) then
757       t1=lst(2,ptlst)
758       a1=lst(3,ptlst)
759       tt1 = nu(a1,t1)
760       t2 = tt1/(2**3)
761       a2 = tt1-t2*(2**3)
762       i11 =   a1 -3
763       i12 =   mod3(i11) 
764       i13 =   mod3(i12)
765       i21 =   a2 -3
766       i22 =   mod3(i21)
767       i23 =   mod3(i22)
768       s1 = nu(i13,t1)
769       s2 = nu(i11,t1)
770       s3 = nu(i12,t1)
771       s4 = nu(i23,t2)
772       x41 = c(1,s4)-c(1,s1)
773       y41 = c(2,s4)-c(2,s1)
774       det2 = (c(1,s2)-c(1,s1))*y41-(c(2,s2)-c(2,s1))*x41
775       det3 = (c(1,s3)-c(1,s1))*y41-(c(2,s3)-c(2,s1))*x41
776       if(det2.gt.0.and.det3.lt.0) then
777       nu(i12,t1) = s4
778       nu(i22,t2) = s1
779       pslst=lst(1,ptlst)
780       if(pslst.gt.0) then
781       aas=lst(3,pslst)
782       if(aas.eq.i22+3) then
783       lst(2,pslst) = t1
784       lst(3,pslst) = i11 + 3
785       endif
786       endif
787       tt1 = nu(i22+3,t2)
788       nu(a1 ,t1) = tt1
789       if(tt1.gt.0) then
790       tt=tt1/(2**3)
791       aa = tt1-(2**3)*tt
792       nu(aa,tt)=  a1 +  (2**3) * t1
793       elseif(tt1.ne.vide) then
794       nu(2,-tt1)= a1 +  (2**3) * t1
795       endif
796       tt1 = nu(i12+3,t1)
797       nu(a2 ,t2) = tt1
798       if(tt1.gt.0) then
799       tt=tt1/(2**3)
800       aa=tt1-(2**3)*tt
801       nu(aa,tt)= a2 +  (2**3) * t2
802       elseif(tt1.ne.vide) then
803       nu(2,-tt1)= a2 +  (2**3) * t2
804       endif
805       nu(i12+3,t1) =   i22+3 + (2**3)*t2
806       nu(i22+3,t2) =   i12+3 + (2**3)*t1
807       det1 = (c(1,s1)-c(1,ss1))*y-(c(2,s1)-c(2,ss1))*x
808       det4 = (c(1,s4)-c(1,ss1))*y-(c(2,s4)-c(2,ss1))*x
809       if(det1.lt.0.and.det4.gt.0) then
810       lst(2,ptlst) = t2
811       lst(3,ptlst) = i22+3
812       elseif(det1.gt.0.and.det4.lt.0) then
813       lst(2,ptlst) = t1
814       lst(3,ptlst) = i12+3
815       else
816       if(pplst.eq.0) then
817       ttlst = lst(1,ptlst)
818       ptlst = ttlst
819       else
820       ptlst        = lst(1,ptlst)
821       lst(1,pplst) = ptlst
822       endif
823       goto 30
824       endif
825       endif
826       pplst = ptlst
827       ptlst = lst(1,ptlst)
828       goto 30
829       endif
830       if(ttlst.ne.0) goto 20
831       nu(i12+3,t1) =  vide
832       nu(i22+3,t2) =  vide
833       t = t2      
834       do 40 i=1,nbac
835       call mshopt (c,nu,lst(2,i),4,nbs,ierr)
836       call mshopt (c,nu,lst(2,i),5,nbs,ierr)
837       call mshopt (c,nu,lst(2,i),6,nbs,ierr)
838 40    continue
839       end