Salome HOME
f66104e34c2d2b10052f1993048587f31a7a459e
[modules/smesh.git] / src / MEFISTO2 / trte.f
1       subroutine qutr2d( p1, p2, p3, qualite )
2 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3 c but :     calculer la qualite d'un triangle de r**2
4 c -----     2 coordonnees des 3 sommets en double precision
5 c
6 c entrees :
7 c ---------
8 c p1,p2,p3 : les 3 coordonnees des 3 sommets du triangle
9 c            sens direct pour une surface et qualite >0
10 c sorties :
11 c ---------
12 c qualite: valeur de la qualite du triangle entre 0 et 1 (equilateral)
13 c          1 etant la qualite optimale
14 c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
15 c auteur : alain perronnet analyse numerique upmc paris     janvier 1995
16 c2345x7..............................................................012
17       parameter  ( d2uxr3 = 3.4641016151377544d0 )
18 c                  d2uxr3 = 2 * sqrt(3)
19       double precision  p1(2), p2(2), p3(2), qualite, a, b, c, p
20 c
21 c     la longueur des 3 cotes
22       a = sqrt( (p2(1)-p1(1))**2 + (p2(2)-p1(2))**2 )
23       b = sqrt( (p3(1)-p2(1))**2 + (p3(2)-p2(2))**2 )
24       c = sqrt( (p1(1)-p3(1))**2 + (p1(2)-p3(2))**2 )
25 c
26 c     demi perimetre
27       p = (a+b+c) * 0.5d0
28 c
29       if ( (a*b*c) .ne. 0d0 ) then
30 c        critere : 2 racine(3) * rayon_inscrit / plus longue arete
31          qualite = d2uxr3 * sqrt( abs( (p-a) / p * (p-b) * (p-c) ) )
32      %          / max(a,b,c)
33       else
34          qualite = 0d0
35       endif
36 c
37 c
38 c     autres criteres possibles:
39 c     critere : 2 * rayon_inscrit / rayon_circonscrit
40 c     qualite = 8d0 * (p-a) * (p-b) * (p-c) / (a * b * c)
41 c
42 c     critere : 3*sqrt(3.) * ray_inscrit / demi perimetre
43 c     qualite = 3*sqrt(3.) * sqrt ((p-a)*(p-b)*(p-c) / p**3)
44 c
45 c     critere : 2*sqrt(3.) * ray_inscrit / max( des aretes )
46 c     qualite = 2*sqrt(3.) * sqrt( (p-a)*(p-b)*(p-c) / p ) / max(a,b,c)
47       end
48
49
50       double precision function surtd2( p1 , p2 , p3 )
51 c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
52 c but : calcul de la surface d'un triangle defini par 3 points de R**2
53 c -----
54 c parametres d entree :
55 c ---------------------
56 c p1 p2 p3 : les 3 fois 2 coordonnees des sommets du triangle
57 c
58 c parametre resultat :
59 c --------------------
60 c surtd2 : surface du triangle
61 c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
62 c auteur : alain perronnet analyse numerique upmc paris     fevrier 1992
63 c2345x7..............................................................012
64       double precision  p1(2), p2(2), p3(2)
65 c
66 c     la surface du triangle
67       surtd2 = ( ( p2(1)-p1(1) ) * ( p3(2)-p1(2) )
68      %         - ( p2(2)-p1(2) ) * ( p3(1)-p1(1) ) ) * 0.5d0
69       end
70
71       integer function nopre3( i )
72 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
73 c but :   numero precedent i dans le sens circulaire  1 2 3 1 ...
74 c -----
75 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
76 c auteur : alain perronnet  analyse numerique paris upmc    fevrier 1992
77 c2345x7..............................................................012
78       if( i .eq. 1 ) then
79          nopre3 = 3
80       else
81          nopre3 = i - 1
82       endif
83       end
84
85       integer function nosui3( i )
86 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
87 c but :   numero suivant i dans le sens circulaire  1 2 3 1 ...
88 c -----
89 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
90 c auteur : alain perronnet  analyse numerique paris upmc    fevrier 1992
91 c2345x7..............................................................012
92       if( i .eq. 3 ) then
93          nosui3 = 1
94       else
95          nosui3 = i + 1
96       endif
97       end
98
99       subroutine provec( v1 , v2 , v3 )
100 c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
101 c but :    v3 vecteur = produit vectoriel de 2 vecteurs de r ** 3
102 c -----
103 c entrees:
104 c --------
105 c v1, v2 : les 2 vecteurs de 3 composantes
106 c
107 c sortie :
108 c --------
109 c v3     : vecteur = v1  produit vectoriel v2
110 cc++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
111 c auteur : perronnet alain upmc analyse numerique paris        mars 1987
112 c2345x7..............................................................012
113       double precision    v1(3), v2(3), v3(3)
114 c
115       v3( 1 ) = v1( 2 ) * v2( 3 ) - v1( 3 ) * v2( 2 )
116       v3( 2 ) = v1( 3 ) * v2( 1 ) - v1( 1 ) * v2( 3 )
117       v3( 3 ) = v1( 1 ) * v2( 2 ) - v1( 2 ) * v2( 1 )
118 c
119       return
120       end
121
122       subroutine norme1( n, v, ierr )
123 c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
124 c but :   normalisation euclidienne a 1 d un vecteur v de n composantes
125 c -----
126 c entrees :
127 c ---------
128 c n       : nombre de composantes du vecteur
129 c
130 c modifie :
131 c ---------
132 c v       : le vecteur a normaliser a 1
133 c
134 c sortie  :
135 c ---------
136 c ierr    : 1 si la norme de v est egale a 0
137 c           0 si pas d'erreur
138 c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
139 c auteur : alain perronnet analyse numerique paris             mars 1987
140 c ......................................................................
141       double precision  v( n ), s, sqrt
142 c
143       s = 0.0d0
144       do 10 i=1,n
145          s = s + v( i ) * v( i )
146    10 continue
147 c
148 c     test de nullite de la norme du vecteur
149 c     --------------------------------------
150       if( s .le. 0.0d0 ) then
151 c        norme nulle du vecteur non normalisable a 1
152          ierr = 1
153          return
154       endif
155 c
156       s = 1.0d0 / sqrt( s )
157       do 20 i=1,n
158          v( i ) = v ( i ) * s
159    20 continue
160 c
161       ierr = 0
162       end
163
164
165       subroutine insoar( mxsomm, mosoar, mxsoar, n1soar, nosoar )
166 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
167 c but :    initialiser le tableau nosoar pour le hachage des aretes
168 c -----
169 c
170 c entrees:
171 c --------
172 c mxsomm : plus grand numero de sommet d'une arete au cours du calcul
173 c mosoar : nombre maximal d'entiers par arete du tableau nosoar
174 c mxsoar : nombre maximal d'aretes stockables dans le tableau nosoar
175 c          avec mxsoar>=3*mxsomm
176 c
177 c sorties:
178 c --------
179 c n1soar : numero de la premiere arete vide dans le tableau nosoar
180 c          une arete i de nosoar est vide  <=>  nosoar(1,i)=0
181 c          chainage des aretes vides amont et aval
182 c          l'arete vide qui precede=nosoar(4,i)
183 c          l'arete vide qui suit   =nosoar(5,i)
184 c nosoar : numero des 2 sommets, no ligne, 2 triangles de l'arete,
185 c          chainage momentan'e d'aretes, chainage du hachage des aretes
186 c          hachage des aretes = min( nosoar(1), nosoar(2) )
187 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
188 c auteur : alain perronnet  analyse numerique paris upmc       mars 1997
189 c2345x7..............................................................012
190       integer   nosoar(mosoar,mxsoar)
191 c
192 c     initialisation des aretes 1 a mxsomm
193       do 10 i=1,mxsomm
194 c
195 c        sommet 1 = 0 <=> temoin d'arete vide pour le hachage
196          nosoar( 1, i ) = 0
197 c
198 c        arete sur aucune ligne
199          nosoar( 3, i ) = 0
200 c
201 c        la position de l'arete interne ou frontaliere est inconnue
202          nosoar( 6, i ) = -2
203 c
204 c        fin de chainage du hachage pas d'arete suivante
205          nosoar( mosoar, i ) = 0
206 c
207  10   continue
208 c
209 c     la premiere arete vide chainee est la mxsomm+1 du tableau
210 c     car ces aretes ne sont pas atteignables par le hachage direct
211       n1soar = mxsomm + 1
212 c
213 c     initialisation des aretes vides et des chainages
214       do 20 i = n1soar, mxsoar
215 c
216 c        sommet 1 = 0 <=> temoin d'arete vide pour le hachage
217          nosoar( 1, i ) = 0
218 c
219 c        arete sur aucune ligne
220          nosoar( 3, i ) = 0
221 c
222 c        chainage sur l'arete vide qui precede
223 c        (si arete occupee cela deviendra le no du triangle 1 de l'arete)
224          nosoar( 4, i ) = i-1
225 c
226 c        chainage sur l'arete vide qui suit
227 c        (si arete occupee cela deviendra le no du triangle 2 de l'arete)
228          nosoar( 5, i ) = i+1
229 c
230 c        chainages des aretes frontalieres ou internes ou ...
231          nosoar( 6, i ) = -2
232 c
233 c        fin de chainage du hachage
234          nosoar( mosoar, i ) = 0
235 c
236  20   continue
237 c
238 c     la premiere arete vide n'a pas de precedent
239       nosoar( 4, n1soar ) = 0
240 c
241 c     la derniere arete vide est mxsoar sans arete vide suivante
242       nosoar( 5, mxsoar ) = 0
243       end
244
245
246       subroutine azeroi ( l , ntab )
247 c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
248 c but : initialisation a zero d un tableau ntab de l variables entieres
249 c -----
250 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
251 c auteur : alain perronnet analyse numerique upmc paris septembre 1988
252 c23456---------------------------------------------------------------012
253       integer ntab(l)
254       do 1 i = 1 , l
255          ntab( i ) = 0
256     1 continue
257       end
258
259
260       subroutine fasoar( ns1,    ns2,    nt1,    nt2,    nolign,
261      %                   mosoar, mxsoar, n1soar, nosoar, noarst,
262      %                   noar,   ierr )
263 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
264 c but :    former l'arete de sommet ns1-ns2 dans le hachage du tableau
265 c -----    nosoar des aretes de la triangulation
266 c
267 c entrees:
268 c --------
269 c ns1 ns2: numero pxyd des 2 sommets de l'arete
270 c nt1    : numero du triangle auquel appartient l'arete
271 c          nt1=-1 si numero inconnu
272 c nt2    : numero de l'eventuel second triangle de l'arete si connu
273 c          nt2=-1 si numero inconnu
274 c nolign : numero de la ligne de l'arete dans ladefi(wulftr-1+nolign)
275 c          =0 si l'arete n'est une arete de ligne
276 c          ce numero est ajoute seulement si l'arete est creee
277 c mosoar : nombre maximal d'entiers par arete du tableau nosoar
278 c mxsoar : nombre maximal d'aretes stockables dans le tableau nosoar
279 c
280 c modifies:
281 c ---------
282 c n1soar : numero de la premiere arete vide dans le tableau nosoar
283 c          une arete i de nosoar est vide  <=>  nosoar(1,i)=0
284 c          chainage des aretes vides amont et aval
285 c          l'arete vide qui precede=nosoar(4,i)
286 c          l'arete vide qui suit   =nosoar(5,i)
287 c nosoar : numero des 2 sommets, no ligne, 2 triangles de l'arete,
288 c          chainage momentan'e d'aretes, chainage du hachage des aretes
289 c          hachage des aretes = min( nosoar(1), nosoar(2) )
290 c noarst : noarst(np) numero d'une arete du sommet np
291 c
292 c ierr   : si < 0  en entree pas d'affichage en cas d'erreur du type
293 c         "arete appartenant a plus de 2 triangles et a creer!"
294 c          si >=0  en entree       affichage de ce type d'erreur
295 c
296 c sorties:
297 c --------
298 c noar   : >0 numero de l'arete retrouvee ou ajoutee
299 c ierr   : =0 si pas d'erreur
300 c          =1 si le tableau nosoar est sature
301 c          =2 si arete a creer et appartenant a 2 triangles distincts
302 c             des triangles nt1 et nt2
303 c          =3 si arete appartenant a 2 triangles distincts
304 c             differents des triangles nt1 et nt2
305 c          =4 si arete appartenant a 2 triangles distincts
306 c             dont le second n'est pas le triangle nt2
307 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
308 c auteur : alain perronnet  analyse numerique paris upmc       mars 1997
309 c2345x7..............................................................012
310       common / unites / lecteu, imprim, nunite(30)
311       integer           nosoar(mosoar,mxsoar), noarst(*)
312       integer           nu2sar(2)
313 c
314 c     ajout eventuel de l'arete s1 s2 dans nosoar
315       nu2sar(1) = ns1
316       nu2sar(2) = ns2
317 c
318 c     hachage de l'arete de sommets nu2sar
319       call hasoar( mosoar, mxsoar, n1soar, nosoar, nu2sar, noar )
320 c     en sortie: noar>0 => no arete retrouvee
321 c                    <0 => no arete ajoutee
322 c                    =0 => saturation du tableau nosoar
323 c
324       if( noar .eq. 0 ) then
325 c
326 c        saturation du tableau nosoar
327          write(imprim,*) 'fasoar: tableau nosoar sature'
328          ierr = 1
329          return
330 c
331       else if( noar .lt. 0 ) then
332 c
333 c        l'arete a ete ajoutee. initialisation des autres informations
334          noar = -noar
335 c        le numero de la ligne de l'arete
336          nosoar(3,noar) = nolign
337 c        le triangle 1 de l'arete => le triangle nt1
338          nosoar(4,noar) = nt1
339 c        le triangle 2 de l'arete => le triangle nt2
340          nosoar(5,noar) = nt2
341 c
342 c        le sommet appartient a l'arete noar
343          noarst( nu2sar(1) ) = noar
344          noarst( nu2sar(2) ) = noar
345 c
346       else
347 c
348 c        l'arete a ete retrouvee.
349 c        si elle appartient a 2 triangles differents de nt1 et nt2
350 c        alors il y a une erreur
351          if( nosoar(4,noar) .gt. 0 .and.
352      %       nosoar(5,noar) .gt. 0 ) then
353              if( nosoar(4,noar) .ne. nt1 .and.
354      %           nosoar(4,noar) .ne. nt2 .or.
355      %           nosoar(5,noar) .ne. nt1 .and.
356      %           nosoar(5,noar) .ne. nt2 ) then
357 c                arete appartenant a plus de 2 triangles => erreur
358                  if( ierr .ge. 0 ) then
359                     write(imprim,*) 'erreur fasoar: arete ',noar,
360      %              ' dans 2 triangles et a creer!'
361                  endif
362                  ierr = 2
363                  return
364              endif
365          endif
366 c
367 c        mise a jour du numero des triangles de l'arete noar
368 c        le triangle 2 de l'arete => le triangle nt1
369          if( nosoar(4,noar) .lt. 0 ) then
370 c            pas de triangle connu pour cette arete
371              n = 4
372          else
373 c            deja un triangle connu. ce nouveau est le second
374              if( nosoar(5,noar) .gt. 0  .and.  nt1 .gt. 0 .and.
375      %          nosoar(5,noar) .ne. nt1 ) then
376 c               arete appartenant a plus de 2 triangles => erreur
377                 write(imprim,*) 'erreur fasoar: arete ',noar,
378      %          ' dans plus de 2 triangles'
379                 ierr = 3
380                 return
381              endif
382              n = 5
383          endif
384          nosoar(n,noar) = nt1
385 c
386 c        cas de l'arete frontaliere retrouvee comme diagonale d'un quadrangle
387          if( nt2 .gt. 0 ) then
388 c           l'arete appartient a 2 triangles
389             if( nosoar(5,noar) .gt. 0  .and.
390      %          nosoar(5,noar) .ne. nt2 ) then
391 c               arete appartenant a plus de 2 triangles => erreur
392                 write(imprim,*) 'erreur fasoar: arete ',noar,
393      %         ' dans plus de 2 triangles'
394                 ierr = 4
395                 return
396             endif
397             nosoar(5,noar) = nt2
398          endif
399 c
400       endif
401 c
402 c     pas d'erreur
403       ierr = 0
404       end
405
406       subroutine fq1inv( x, y, s, xc, yc, ierr )
407 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
408 c but :   calcul des 2 coordonnees (xc,yc) dans le carre (0,1)
409 c -----   image par f:carre unite-->quadrangle appartenant a q1**2
410 c         par une resolution directe due a nicolas thenault
411 c
412 c entrees:
413 c --------
414 c x,y   : coordonnees du point image dans le quadrangle de sommets s
415 c s     : les 2 coordonnees des 4 sommets du quadrangle
416 c
417 c sorties:
418 c --------
419 c xc,yc : coordonnees dans le carre dont l'image par f vaut (x,y)
420 c ierr  : 0 si calcul sans erreur, 1 si quadrangle degenere
421 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
422 c auteurs: thenault tulenew  analyse numerique paris        janvier 1998
423 c modifs : perronnet alain   analyse numerique paris        janvier 1998
424 c234567..............................................................012
425       real             s(1:2,1:4), dist(2)
426       double precision a,b,c,d,alpha,beta,gamma,delta,x0,y0,t(2),u,v,w
427 c
428       a = s(1,1)
429       b = s(1,2) - s(1,1)
430       c = s(1,4) - s(1,1)
431       d = s(1,1) - s(1,2) + s(1,3) - s(1,4)
432 c
433       alpha = s(2,1)
434       beta  = s(2,2) - s(2,1)
435       gamma = s(2,4) - s(2,1)
436       delta = s(2,1) - s(2,2) + s(2,3) - s(2,4)
437 c
438       u = beta  * c - b * gamma
439       if( u .eq. 0 ) then
440 c        quadrangle degenere
441          ierr = 1
442          return
443       endif
444       v = delta * c - d * gamma
445       w = b * delta - beta * d
446 c
447       x0 = c * (y-alpha) - gamma * (x-a)
448       y0 = b * (y-alpha) - beta  * (x-a)
449 c
450       a = v  * w
451       b = u  * u - w * x0 - v * y0
452       c = x0 * y0
453 c
454       if( a .ne. 0 ) then
455 c
456          delta = sqrt( b*b-4*a*c )
457          if( b .ge. 0.0 ) then
458             t(2) = -b - delta
459          else
460             t(2) = -b + delta
461          endif
462 c        la racine de plus grande valeur absolue
463 c       (elle donne le plus souvent le point exterieur au carre unite
464 c        donc a tester en second pour reduire les calculs)
465          t(2) = t(2) / ( 2 * a )
466 c        calcul de la seconde racine a partir de la somme => plus stable
467          t(1) = - b/a - t(2)
468 c
469          do 10 i=1,2
470 c
471 c           la solution i donne t elle un point interne au carre unite?
472             xc = ( x0 - v * t(i) ) / u
473             yc = ( w * t(i) - y0 ) / u
474             if( 0.0 .le. xc .and. xc .le. 1.0 ) then
475                if( 0.0 .le. yc .and. yc .le. 1.0 ) goto 9000
476             endif
477 c
478 c           le point (xc,yc) n'est pas dans le carre unite
479 c           cela peut etre du aux erreurs d'arrondi
480 c           => choix par le minimum de la distance aux bords du carre
481             dist(i) = max( 0.0, -xc, xc-1.0, -yc, yc-1.0 )
482 c
483  10      continue
484 c
485          if( dist(1) .gt. dist(2) ) then
486 c           f(xc,yc) pour la racine 2 est plus proche de x,y
487 c           xc yc sont deja calcules
488             goto 9000
489          endif
490 c
491       else if ( b .ne. 0 ) then
492          t(1) = - c / b
493       else
494          t(1) = 0
495       endif
496 c
497 c     les 2 coordonnees du point dans le carre unite
498       xc = ( x0 - v * t(1) ) / u
499       yc = ( w * t(1) - y0 ) / u
500 c
501  9000 ierr = 0
502       return
503       end
504
505
506       subroutine ptdatr( point, pxyd, nosotr, nsigne )
507 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
508 c but :    le point est il dans le triangle de sommets nosotr
509 c -----
510 c
511 c entrees:
512 c --------
513 c point  : les 2 coordonnees du point
514 c pxyd   : les 2 coordonnees et distance souhaitee des points du maillage
515 c nosotr : le numero des 3 sommets du triangle
516 c
517 c sorties:
518 c --------
519 c nsigne : >0 si le point est dans le triangle ou sur une des 3 aretes
520 c          =0 si le triangle est degenere ou indirect ou ne contient pas le poin
521 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
522 c auteur : alain perronnet  analyse numerique paris upmc       mars 1997
523 c....................................................................012
524       integer           nosotr(3)
525       double precision  point(2), pxyd(3,*)
526       double precision  xp,yp, x1,x2,x3, y1,y2,y3, d,dd, cb1,cb2,cb3
527 c
528       xp = point( 1 )
529       yp = point( 2 )
530 c
531       n1 = nosotr( 1 )
532       x1 = pxyd( 1 , n1 )
533       y1 = pxyd( 2 , n1 )
534 c
535       n2 = nosotr( 2 )
536       x2 = pxyd( 1 , n2 )
537       y2 = pxyd( 2 , n2 )
538 c
539       n3 = nosotr( 3 )
540       x3 = pxyd( 1 , n3 )
541       y3 = pxyd( 2 , n3 )
542 c
543 c     2 fois la surface du triangle = determinant de la matrice
544 c     de calcul des coordonnees barycentriques du point p
545       d  = ( x2 - x1 ) * ( y3 - y1 ) - ( x3 - x1 ) * ( y2 - y1 )
546 c
547       if( d .gt. 0 ) then
548 c
549 c        triangle non degenere
550 c        =====================
551 c        calcul des 3 coordonnees barycentriques du
552 c        point xp yp dans le triangle
553          cb1 = ( ( x2-xp ) * ( y3-yp ) - ( x3-xp ) * ( y2-yp ) ) / d
554          cb2 = ( ( x3-xp ) * ( y1-yp ) - ( x1-xp ) * ( y3-yp ) ) / d
555          cb3 = 1d0 - cb1 -cb2
556 ccc         cb3 = ( ( x1-xp ) * ( y2-yp ) - ( x2-xp ) * ( y1-yp ) ) / d
557 c
558 ccc         if( cb1 .ge. -0.00005d0 .and. cb1 .le. 1.00005d0 .and.
559          if( cb1 .ge. 0d0 .and. cb1 .le. 1d0 .and.
560      %       cb2 .ge. 0d0 .and. cb2 .le. 1d0 .and.
561      %       cb3 .ge. 0d0 .and. cb3 .le. 1d0 ) then
562 c
563 c           le triangle nosotr contient le point
564             nsigne = 1
565          else
566             nsigne = 0
567          endif
568 c
569       else
570 c
571 c        triangle degenere
572 c        =================
573 c        le point est il du meme cote que le sommet oppose de chaque arete?
574          nsigne = 0
575          do 10 i=1,3
576 c           le sinus de l'angle p1 p2-p1 point
577             x1  = pxyd(1,n1)
578             y1  = pxyd(2,n1)
579             d   = ( pxyd(1,n2) - x1 ) * ( point(2) - y1 )
580      %          - ( pxyd(2,n2) - y1 ) * ( point(1) - x1 )
581             dd  = ( pxyd(1,n2) - x1 ) * ( pxyd(2,n3) - y1 )
582      %          - ( pxyd(2,n2) - y1 ) * ( pxyd(1,n3) - x1 )
583             cb1 = ( pxyd(1,n2) - x1 ) ** 2
584      %          + ( pxyd(2,n2) - y1 ) ** 2
585             cb2 = ( point(1) - x1 ) ** 2
586      %          + ( point(2) - y1 ) ** 2
587             cb3 = ( pxyd(1,n3) - x1 ) ** 2
588      %          + ( pxyd(2,n3) - y1 ) ** 2
589             if( abs( dd ) .le. 1e-4 * sqrt( cb1 * cb3 ) ) then
590 c              le point 3 est sur l'arete 1-2
591 c              le point doit y etre aussi
592                if( abs( d ) .le. 1e-4 * sqrt( cb1 * cb2 ) ) then
593 c                 point sur l'arete
594                   nsigne = nsigne + 1
595                endif
596             else
597 c              le point 3 n'est pas sur l'arete . test des signes
598                if( d * dd .ge. 0 ) then
599                   nsigne = nsigne + 1
600                endif
601             endif
602 c           permutation circulaire des 3 sommets et aretes
603             n  = n1
604             n1 = n2
605             n2 = n3
606             n3 = n
607  10      continue
608          if( nsigne .ne. 3 ) nsigne = 0
609       endif
610       end
611
612       integer function nosstr( p, pxyd, nt, letree )
613 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
614 c but :    calculer le numero 0 a 3 du sous-triangle te contenant
615 c -----    le point p
616 c
617 c entrees:
618 c --------
619 c p      : point de r**2 contenu dans le te nt de letree
620 c pxyd   : x y distance des points
621 c nt     : numero letree du te de te voisin a calculer
622 c letree : arbre-4 des triangles equilateraux (te) fond de la triangulation
623 c      letree(0,0)  no du 1-er te vide dans letree
624 c      letree(0,1) : maximum du 1-er indice de letree (ici 8)
625 c      letree(0,2) : maximum declare du 2-eme indice de letree (ici mxtree)
626 c      letree(0:8,1) : racine de l'arbre  (triangle sans sur triangle)
627 c      si letree(0,.)>0 alors
628 c         letree(0:3,j) : no (>0) letree des 4 sous-triangles du triangle j
629 c      sinon
630 c         letree(0:3,j) :-no pxyd des 1 \85a 4 points internes au triangle j
631 c                         0  si pas de point
632 c                       ( j est alors une feuille de l'arbre )
633 c      letree(4,j) : no letree du sur-triangle du triangle j
634 c      letree(5,j) : 0 1 2 3 no du sous-triangle j pour son sur-triangle
635 c      letree(6:8,j) : no pxyd des 3 sommets du triangle j
636 c
637 c sorties :
638 c ---------
639 c nosstr : 0 si le sous-triangle central contient p
640 c          i =1,2,3 numero du sous-triangle contenant p
641 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
642 c auteur : alain perronnet  analyse numerique paris upmc    fevrier 1992
643 c2345x7..............................................................012
644       integer           letree(0:8,0:*)
645       double precision  pxyd(3,*), p(2),
646      %                  x1, y1, x21, y21, x31, y31, d, xe, ye
647 c
648 c     le numero des 3 sommets du triangle
649       ns1 = letree( 6, nt )
650       ns2 = letree( 7, nt )
651       ns3 = letree( 8, nt )
652 c
653 c     les coordonnees entre 0 et 1 du point p
654       x1  = pxyd(1,ns1)
655       y1  = pxyd(2,ns1)
656 c
657       x21 = pxyd(1,ns2) - x1
658       y21 = pxyd(2,ns2) - y1
659 c
660       x31 = pxyd(1,ns3) - x1
661       y31 = pxyd(2,ns3) - y1
662 c
663       d   = 1.0 / ( x21 * y31 - x31 * y21 )
664 c
665       xe  = ( ( p(1) - x1 ) * y31 - ( p(2) - y1 ) * x31 ) * d
666       ye  = ( ( p(2) - y1 ) * x21 - ( p(1) - x1 ) * y21 ) * d
667 c
668       if( xe .gt. 0.5d0 ) then
669 c        sous-triangle droit
670          nosstr = 2
671       else if( ye .gt. 0.5d0 ) then
672 c        sous-triangle haut
673          nosstr = 3
674       else if( xe+ye .lt. 0.5d0 ) then
675 c        sous-triangle gauche
676          nosstr = 1
677       else
678 c        sous-triangle central
679          nosstr = 0
680       endif
681       end
682
683
684       integer function notrpt( p, pxyd, notrde, letree )
685 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
686 c but :    calculer le numero letree du sous-triangle feuille contenant
687 c -----    le point p a partir du te notrde de letree
688 c
689 c entrees:
690 c --------
691 c p      : point de r**2 contenu dans le te nt de letree
692 c pxyd   : x y distance des points
693 c notrde : numero letree du triangle depart de recherche (1=>racine)
694 c letree : arbre-4 des triangles equilateraux (te) fond de la triangulation
695 c      letree(0,0)  no du 1-er te vide dans letree
696 c      letree(0,1) : maximum du 1-er indice de letree (ici 8)
697 c      letree(0,2) : maximum declare du 2-eme indice de letree (ici mxtree)
698 c      letree(0:8,1) : racine de l'arbre  (triangle sans sur triangle)
699 c      si letree(0,.)>0 alors
700 c         letree(0:3,j) : no (>0) letree des 4 sous-triangles du triangle j
701 c      sinon
702 c         letree(0:3,j) :-no pxyd des 1 \85 4 points internes au triangle j
703 c                         0  si pas de point
704 c                        ( j est alors une feuille de l'arbre )
705 c      letree(4,j) : no letree du sur-triangle du triangle j
706 c      letree(5,j) : 0 1 2 3 no du sous-triangle j pour son sur-triangle
707 c      letree(6:8,j) : no pxyd des 3 sommets du triangle j
708 c
709 c sorties :
710 c ---------
711 c notrpt : numero letree du triangle contenant le point p
712 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
713 c auteur : alain perronnet  analyse numerique paris upmc    fevrier 1992
714 c2345x7..............................................................012
715       integer           letree(0:8,0:*)
716       double precision  pxyd(1:3,*), p(2)
717 c
718 c     la racine depart de la recherche
719       notrpt = notrde
720 c
721 c     tant que la feuille n'est pas atteinte descendre l'arbre
722  10   if( letree(0,notrpt) .gt. 0 ) then
723 c
724 c        recherche du sous-triangle contenant p
725          nsot = nosstr( p, pxyd, notrpt, letree )
726 c
727 c        le numero letree du sous-triangle
728          notrpt = letree( nsot, notrpt )
729          goto 10
730 c
731       endif
732       end
733
734
735       subroutine teajpt( ns,   nbsomm, mxsomm, pxyd, letree,
736      &                   ntrp, ierr )
737 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
738 c but :    ajout du point ns de pxyd dans letree
739 c -----
740 c
741 c entrees:
742 c --------
743 c ns     : numero du point a ajouter dans letree
744 c mxsomm : nombre maximal de points declarables dans pxyd
745 c pxyd   : tableau des coordonnees des points
746 c          par point : x  y  distance_souhaitee
747 c
748 c modifies :
749 c ----------
750 c nbsomm : nombre actuel de points dans pxyd
751 c
752 c letree : arbre-4 des triangles equilateraux (te) fond de la triangulation
753 c      letree(0,0) : no du 1-er te vide dans letree
754 c      letree(0,1) : maximum du 1-er indice de letree (ici 8)
755 c      letree(0,2) : maximum declare du 2-eme indice de letree (ici mxtree)
756 c      letree(0:8,1) : racine de l'arbre  (triangle sans sur triangle)
757 c      si letree(0,.)>0 alors
758 c         letree(0:3,j) : no (>0) letree des 4 sous-triangles du triangle j
759 c      sinon
760 c         letree(0:3,j) :-no pxyd des 1 \85a 4 points internes au triangle j
761 c                         0  si pas de point
762 c                        ( j est alors une feuille de l'arbre )
763 c      letree(4,j) : no letree du sur-triangle du triangle j
764 c      letree(5,j) : 0 1 2 3 no du sous-triangle j pour son sur-triangle
765 c      letree(6:8,j) : no pxyd des 3 sommets du triangle j
766 c
767 c sorties :
768 c ---------
769 c ntrp    : numero letree du triangle te ou a ete ajoute le point
770 c ierr    : 0 si pas d'erreur,  51 saturation letree, 52 saturation pxyd
771 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
772 c auteur : alain perronnet  analyse numerique paris upmc    fevrier 1992
773 c2345x7..............................................................012
774       integer           letree(0:8,0:*)
775       double precision  pxyd(3,mxsomm)
776 c
777 c     depart de la racine
778       ntrp = 1
779 c
780 c     recherche du triangle contenant le point pxyd(ns)
781  1    ntrp = notrpt( pxyd(1,ns), pxyd, ntrp, letree )
782 c
783 c     existe t il un point libre
784       do 10 i=0,3
785          if( letree(i,ntrp) .eq. 0 ) then
786 c           la place i est libre
787             letree(i,ntrp) = -ns
788             return
789          endif
790  10   continue
791 c
792 c     pas de place libre => 4 sous-triangles sont crees
793 c                           a partir des 3 milieux des aretes
794       call te4ste( nbsomm, mxsomm, pxyd, ntrp, letree, ierr )
795       if( ierr .ne. 0 ) return
796 c
797 c     ajout du point ns
798       goto 1
799       end
800
801       subroutine n1trva( nt, lar, letree, notrva, lhpile )
802 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
803 c but :    calculer le numero letree du triangle voisin du te nt
804 c -----    par l'arete lar (1 a 3 ) de nt
805 c          attention : notrva n'est pas forcement minimal
806 c
807 c entrees:
808 c --------
809 c nt     : numero letree du te de te voisin a calculer
810 c lar    : numero 1 a 3 de l'arete du triangle nt
811 c letree : arbre-4 des triangles equilateraux (te) fond de la triangulation
812 c   letree(0,0)  no du 1-er te vide dans letree
813 c   letree(0,1) : maximum du 1-er indice de letree (ici 8)
814 c   letree(0,2) : maximum declare du 2-eme indice de letree (ici mxtree)
815 c   letree(0:8,1) : racine de l'arbre  (triangle sans sur-triangle)
816 c   si letree(0,.)>0 alors
817 c      letree(0:3,j) : no (>0) letree des 4 sous-triangles du triangle j
818 c   sinon
819 c      letree(0:3,j) :-no pxyd des 1 a 4 points internes au triangle j
820 c                      0  si pas de point
821 c                     ( j est alors une feuille de l'arbre )
822 c   letree(4,j) : no letree du sur-triangle du triangle j
823 c   letree(5,j) : 0 1 2 3 no du sous-triangle j pour son sur-triangle
824 c   letree(6:8,j) : no pxyd des 3 sommets du triangle j
825 c
826 c sorties :
827 c ---------
828 c notrva  : >0 numero letree du te voisin par l'arete lar
829 c           =0 si pas de te voisin (racine , ... )
830 c lhpile  : =0 si nt et notrva ont meme taille
831 c           >0 nt est 4**lhpile fois plus petit que notrva
832 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
833 c auteur : alain perronnet  analyse numerique paris upmc    fevrier 1992
834 c2345x7..............................................................012
835       integer   letree(0:8,0:*)
836       integer   lapile(1:64)
837 c
838 c     initialisation de la pile
839 c     le triangle est empile
840       lapile(1) = nt
841       lhpile = 1
842 c
843 c     tant qu'il existe un sur-triangle
844  10   ntr  = lapile( lhpile )
845       if( ntr .eq. 1 ) then
846 c        racine atteinte => pas de triangle voisin
847          notrva = 0
848          lhpile = lhpile - 1
849          return
850       endif
851 c
852 c     le type du triangle ntr
853       nty  = letree( 5, ntr )
854 c     l'eventuel sur-triangle
855       nsut = letree( 4, ntr )
856 c
857       if( nty .eq. 0 ) then
858 c
859 c        triangle de type 0 => triangle voisin de type precedent(lar)
860 c                              dans le sur-triangle de ntr
861 c                              ce triangle remplace ntr dans lapile
862          lapile( lhpile ) = letree( nopre3(lar), nsut )
863          goto 20
864       endif
865 c
866 c     triangle ntr de type nty>0
867       if( nosui3(nty) .eq. lar ) then
868 c
869 c        le triangle voisin par lar est le triangle 0
870          lapile( lhpile ) = letree( 0, nsut )
871          goto 20
872       endif
873 c
874 c     triangle sans voisin direct => passage par le sur-triangle
875       if( nsut .eq. 0 ) then
876 c
877 c        ntr est la racine => pas de triangle voisin par cette arete
878          notrva = 0
879          return
880       else
881 c
882 c        le sur-triangle est empile
883          lhpile = lhpile + 1
884          lapile(lhpile) = nsut
885          goto 10
886       endif
887 c
888 c     descente aux sous-triangles selon la meme arete
889  20   notrva = lapile( lhpile )
890 c
891  30   lhpile = lhpile - 1
892       if( letree(0,notrva) .le. 0 ) then
893 c        le triangle est une feuille de l'arbre 0 sous-triangle
894 c        lhpile = nombre de differences de niveaux dans l'arbre
895          return
896       else
897 c        le triangle a 4 sous-triangles
898          if( lhpile .gt. 0 ) then
899 c
900 c           bas de pile non atteint
901             nty  = letree( 5, lapile(lhpile) )
902             if( nty .eq. lar ) then
903 c              l'oppose est suivant(nty) de notrva
904                notrva = letree( nosui3(nty) , notrva )
905             else
906 c              l'oppose est precedent(nty) de notrva
907                notrva = letree( nopre3(nty) , notrva )
908             endif
909             goto 30
910          endif
911       endif
912 c
913 c     meme niveau dans l'arbre lhpile = 0
914       end
915
916
917       subroutine cenced( xy1, xy2, xy3, cetria, ierr )
918 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
919 c but : calcul des coordonnees du centre du cercle circonscrit
920 c ----- du triangle defini par ses 3 sommets de coordonnees
921 c       xy1 xy2 xy3 ainsi que le carre du rayon de ce cercle
922 c
923 c entrees :
924 c ---------
925 c xy1 xy2 xy3 : les 2 coordonnees des 3 sommets du triangle
926 c ierr   : <0  => pas d'affichage si triangle degenere
927 c          >=0 =>       affichage si triangle degenere
928 c
929 c sortie :
930 c --------
931 c cetria : cetria(1)=abcisse  du centre
932 c          cetria(2)=ordonnee du centre
933 c          cetria(3)=carre du rayon   1d28 si triangle degenere
934 c ierr   : 0 si triangle non degenere
935 c          1 si triangle degenere
936 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
937 c auteur : perronnet alain upmc analyse numerique paris        juin 1995
938 c2345x7..............................................................012
939       parameter        (epsurf=1d-7)
940       common / unites / lecteu,imprim,nunite(30)
941       double precision  x1,y1,x21,y21,x31,y31,
942      %                  aire2,xc,yc,rot,
943      %                  xy1(2),xy2(2),xy3(2),cetria(3)
944 c
945 c     le calcul de 2 fois l'aire du triangle
946 c     attention l'ordre des 3 sommets est direct ou non
947       x1  = xy1(1)
948       x21 = xy2(1) - x1
949       x31 = xy3(1) - x1
950 c
951       y1  = xy1(2)
952       y21 = xy2(2) - y1
953       y31 = xy3(2) - y1
954 c
955       aire2  = x21 * y31 - x31 * y21
956 c
957 c     recherche d'un test relatif peu couteux
958 c     pour reperer la degenerescence du triangle
959       if( abs(aire2) .le.
960      %    epsurf*(abs(x21)+abs(x31))*(abs(y21)+abs(y31)) ) then
961 c        triangle de qualite trop faible
962          if( ierr .ge. 0 ) then
963 c            nblgrc(nrerr) = 1
964 c            kerr(1) = 'erreur cenced: triangle degenere'
965 c            call lereur
966             write(imprim,*) 'erreur cenced: triangle degenere'
967             write(imprim,10000)  xy1,xy2,xy3,aire2
968          endif
969 10000 format( 3(' x=',g24.16,' y=',g24.16/),' aire*2=',g24.16)
970          cetria(1) = 0d0
971          cetria(2) = 0d0
972          cetria(3) = 1d28
973          ierr = 1
974          return
975       endif
976 c
977 c     les 2 coordonnees du centre intersection des 2 mediatrices
978 c     x = (x1+x2)/2 + lambda * (y2-y1)
979 c     y = (y1+y2)/2 - lambda * (x2-x1)
980 c     x = (x1+x3)/2 + rot    * (y3-y1)
981 c     y = (y1+y3)/2 - rot    * (x3-x1)
982 c     ==========================================================
983       rot = ((xy2(1)-xy3(1))*x21 + (xy2(2)-xy3(2))*y21) / (2 * aire2)
984 c
985       xc = ( x1 + xy3(1) ) * 0.5d0 + rot * y31
986       yc = ( y1 + xy3(2) ) * 0.5d0 - rot * x31
987 c
988       cetria(1) = xc
989       cetria(2) = yc
990 c
991 c     le carre du rayon
992       cetria(3) = (x1-xc) ** 2 + (y1-yc) ** 2
993 c
994 c     pas d'erreur rencontree
995       ierr = 0
996       end
997
998
999       double precision function angled( p1, p2, p3 )
1000 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1001 c but :   calculer l'angle (p1p2,p1p3) en radians
1002 c -----
1003 c
1004 c entrees :
1005 c ---------
1006 c p1,p2,p3 : les 2 coordonnees des 3 sommets de l'angle
1007 c               sens direct pour une surface >0
1008 c sorties :
1009 c ---------
1010 c angled :  angle (p1p2,p1p3) en radians entre [0 et 2pi]
1011 c           0 si p1=p2 ou p1=p3
1012 c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1013 c auteur : alain perronnet analyse numerique upmc paris     fevrier 1992
1014 c2345x7..............................................................012
1015       double precision  p1(2),p2(2),p3(2),x21,y21,x31,y31,a1,a2,d,c
1016 c
1017 c     les cotes
1018       x21 = p2(1) - p1(1)
1019       y21 = p2(2) - p1(2)
1020       x31 = p3(1) - p1(1)
1021       y31 = p3(2) - p1(2)
1022 c
1023 c     longueur des cotes
1024       a1 = x21 * x21 + y21 * y21
1025       a2 = x31 * x31 + y31 * y31
1026       d  = sqrt( a1 * a2 )
1027       if( d .eq. 0 ) then
1028          angled = 0
1029          return
1030       endif
1031 c
1032 c     cosinus de l'angle
1033       c  = ( x21 * x31 + y21 * y31 ) / d
1034       if( c .le. -1.d0 ) then
1035 c        tilt sur apollo si acos( -1 -eps )
1036          angled = atan( 1.d0 ) * 4.d0
1037          return
1038       else if( c .ge. 1.d0 ) then
1039 c        tilt sur apollo si acos( 1 + eps )
1040          angled = 0
1041          return
1042       endif
1043 c
1044       angled = acos( c )
1045       if( x21 * y31 - x31 * y21 .lt. 0 ) then
1046 c        demi plan inferieur
1047          angled = 8.d0 * atan( 1.d0 ) - angled
1048       endif
1049       end
1050
1051
1052       subroutine teajte( mxsomm, nbsomm, pxyd,   comxmi,
1053      %                   aretmx, mxtree, letree,
1054      %                   ierr )
1055 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1056 c but :    initialisation des tableaux letree
1057 c -----    ajout des sommets 1 a nbsomm (valeur en entree) dans letree
1058 c
1059 c entrees:
1060 c --------
1061 c mxsomm : nombre maximal de sommets permis pour la triangulation
1062 c mxtree : nombre maximal de triangles equilateraux (te) declarables
1063 c aretmx : longueur maximale des aretes des triangles equilateraux
1064 c
1065 c entrees et sorties :
1066 c --------------------
1067 c nbsomm : nombre de sommets apres identification
1068 c pxyd   : tableau des coordonnees 2d des points
1069 c          par point : x  y  distance_souhaitee
1070 c          tableau reel(3,mxsomm)
1071 c
1072 c sorties:
1073 c --------
1074 c comxmi : coordonnees minimales et maximales des points frontaliers
1075 c letree : arbre-4 des triangles equilateraux (te) fond de la triangulation
1076 c          letree(0,0) : no du 1-er te vide dans letree
1077 c          letree(0,1) : maximum du 1-er indice de letree (ici 8)
1078 c          letree(0,2) : maximum declare du 2-eme indice de letree (ici mxtree)
1079 c          letree(0:8,1) : racine de l'arbre  (triangle sans sur triangle)
1080 c          si letree(0,.)>0 alors
1081 c             letree(0:3,j) : no (>0) letree des 4 sous-triangles du triangle j
1082 c          sinon
1083 c             letree(0:3,j) :-no pxyd des 1 a 4 points internes au triangle j
1084 c                             0  si pas de point
1085 c                             ( j est alors une feuille de l'arbre )
1086 c          letree(4,j) : no letree du sur-triangle du triangle j
1087 c          letree(5,j) : 0 1 2 3 no du sous-triangle j pour son sur-triangle
1088 c          letree(6:8,j) : no pxyd des 3 sommets du triangle j
1089 c
1090 c ierr   :  0 si pas d'erreur
1091 c          51 saturation letree
1092 c          52 saturation pxyd
1093 c           7 tous les points sont alignes
1094 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1095 c auteur : alain perronnet  analyse numerique paris upmc    juillet 1994
1096 c....................................................................012
1097       integer           letree(0:8,0:mxtree)
1098       double precision  pxyd(3,mxsomm)
1099       double precision  comxmi(3,2)
1100       double precision  a(2),s,aretmx,rac3
1101 c
1102 c     protection du nombre de sommets avant d'ajouter ceux de tetree
1103       nbsofr = nbsomm
1104       do 1 i = 1, nbsomm 
1105          comxmi(1,1) = min( comxmi(1,1), pxyd(1,i) )
1106          comxmi(1,2) = max( comxmi(1,2), pxyd(1,i) )
1107          comxmi(2,1) = min( comxmi(2,1), pxyd(2,i) )
1108          comxmi(2,2) = max( comxmi(2,2), pxyd(2,i) )
1109  1    continue
1110 c
1111 c     creation de l'arbre tee
1112 c     =======================
1113 c     la premiere colonne vide de letree
1114       letree(0,0) = 2
1115 c     chainage des te vides
1116       do 4 i = 2 , mxtree
1117          letree(0,i) = i+1
1118  4    continue
1119       letree(0,mxtree) = 0
1120 c     les maxima des 2 indices de letree
1121       letree(1,0) = 8
1122       letree(2,0) = mxtree
1123 c
1124 c     la racine
1125 c     aucun point interne au triangle equilateral (te) 1
1126       letree(0,1) = 0
1127       letree(1,1) = 0
1128       letree(2,1) = 0
1129       letree(3,1) = 0
1130 c     pas de sur-triangle
1131       letree(4,1) = 0
1132       letree(5,1) = 0
1133 c     le numero pxyd des 3 sommets du te 1
1134       letree(6,1) = nbsomm + 1
1135       letree(7,1) = nbsomm + 2
1136       letree(8,1) = nbsomm + 3
1137 c
1138 c     calcul de la largeur et hauteur du rectangle englobant
1139 c     ======================================================
1140       a(1) = comxmi(1,2) - comxmi(1,1)
1141       a(2) = comxmi(2,2) - comxmi(2,1)
1142 c     la longueur de la diagonale
1143       s = sqrt( a(1)**2 + a(2)**2 )
1144       do 60 k=1,2
1145          if( a(k) .lt. 1e-4 * s ) then
1146 c            nblgrc(nrerr) = 1
1147             write(imprim,*) 'tous les points sont alignes'
1148 c            call lereur
1149             ierr = 7
1150             return
1151          endif
1152  60   continue
1153 c
1154 c     le maximum des ecarts
1155       s = s + s
1156 c
1157 c     le triangle equilateral englobant
1158 c     =================================
1159 c     ecart du rectangle au triangle equilateral
1160       rac3 = sqrt( 3.0d0 )
1161       arete = a(1) + 2 * aretmx + 2 * ( a(2) + aretmx ) / rac3
1162 c
1163 c     le point nbsomm + 1 en bas a gauche
1164       nbsomm = nbsomm + 1
1165       pxyd(1,nbsomm) = (comxmi(1,1)+comxmi(1,2))*0.5d0 - arete*0.5d0
1166       pxyd(2,nbsomm) =  comxmi(2,1) - aretmx
1167       pxyd(3,nbsomm) = s
1168 c
1169 c     le point nbsomm + 2 en bas a droite
1170       nbsomm = nbsomm + 1
1171       pxyd(1,nbsomm) = pxyd(1,nbsomm-1) + arete
1172       pxyd(2,nbsomm) = pxyd(2,nbsomm-1)
1173       pxyd(3,nbsomm) = s
1174 c
1175 c     le point nbsomm + 3 sommet au dessus
1176       nbsomm = nbsomm + 1
1177       pxyd(1,nbsomm) = pxyd(1,nbsomm-2) + arete * 0.5d0
1178       pxyd(2,nbsomm) = pxyd(2,nbsomm-2) + arete * 0.5d0 * rac3
1179       pxyd(3,nbsomm) = s
1180 c
1181 c     ajout des sommets des lignes pour former letree
1182 c     ===============================================
1183       do 150 i=1,nbsofr
1184 c        ajout du point i de pxyd a letree
1185          call teajpt(  i, nbsomm, mxsomm, pxyd, letree,
1186      &                nt, ierr )
1187          if( ierr .ne. 0 ) return
1188  150  continue
1189 c
1190       return
1191       end
1192
1193
1194       subroutine tetaid( nutysu, dx, dy, longai, ierr )
1195 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1196 c but :     calculer la longueur de l'arete ideale en dx,dy
1197 c -----
1198 c
1199 c entrees:
1200 c --------
1201 c nutysu : numero de traitement de areteideale() selon le type de surface
1202 c          0 pas d'emploi de la fonction areteideale() => aretmx active
1203 c          1 il existe une fonction areteideale(xyz,xyzdir)
1204 c          ... autres options a definir ...
1205 c dx, dy : abscisse et ordonnee dans le plan du point (reel2!)
1206 c
1207 c sorties:
1208 c --------
1209 c longai : longueur de l'areteideale(xyz,xyzdir) autour du point xyz
1210 c ierr   : 0 si pas d'erreur, <>0 sinon
1211 c          1 calcul incorrect de areteideale(xyz,xyzdir)
1212 c          2 longueur calculee nulle
1213 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1214 c auteur : alain perronnet  analyse numerique paris upmc       mars 1997
1215 c2345x7..............................................................012
1216       common / unites / lecteu, imprim, nunite(30)
1217 c
1218       double precision  areteideale
1219       double precision  dx, dy, longai
1220       double precision  xyz(3), xyzd(3), d0
1221 c
1222       ierr = 0
1223       if( nutysu .gt. 0 ) then
1224          d0 = longai
1225 c        le point ou se calcule la longueur
1226          xyz(1) = dx
1227          xyz(2) = dy
1228 c        z pour le calcul de la longueur (inactif ici!)
1229          xyz(3) = 0d0
1230 c        la direction pour le calcul de la longueur (inactif ici!)
1231          xyzd(1) = 0d0
1232          xyzd(2) = 0d0
1233          xyzd(3) = 0d0
1234
1235          longai = areteideale(xyz,xyzd)
1236          if( longai .lt. 0d0 ) then
1237             write(imprim,10000) xyz
1238 10000       format('attention: longueur de areteideale(',
1239      %              g14.6,',',g14.6,',',g14.6,')<=0! => rendue >0' )
1240             longai = -longai
1241          endif
1242          if( longai .eq. 0d0 ) then
1243             write(imprim,10001) xyz
1244 10001       format('erreur: longueur de areteideale(',
1245      %              g14.6,',',g14.6,',',g14.6,')=0!' )
1246             ierr = 2
1247             longai = d0
1248          endif
1249       endif
1250       end
1251
1252
1253       subroutine tehote( nutysu,
1254      %                   nbarpi, mxsomm, nbsomm, pxyd,
1255      %                   comxmi, aretmx,
1256      %                   letree, mxqueu, laqueu,
1257      %                   ierr )
1258 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1259 c but :     homogeneisation de l'arbre des te a un saut de taille au plus
1260 c -----     prise en compte des distances souhaitees autour des sommets initiaux
1261 c
1262 c entrees:
1263 c --------
1264 c nutysu : numero de traitement de areteideale() selon le type de surface
1265 c          0 pas d'emploi de la fonction areteideale() => aretmx active
1266 c          1 il existe une fonction areteideale()
1267 c            dont seules les 2 premieres composantes de uv sont actives
1268 c          autres options a definir...
1269 c nbarpi : nombre de sommets de la frontiere + nombre de points internes
1270 c          imposes par l'utilisateur
1271 c mxsomm : nombre maximal de sommets permis pour la triangulation  et te
1272 c mxqueu : nombre d'entiers utilisables dans laqueu
1273 c comxmi : minimum et maximum des coordonnees de l'objet
1274 c aretmx : longueur maximale des aretes des triangles equilateraux
1275 c permtr : perimetre de la ligne enveloppe dans le plan
1276 c          avant mise a l'echelle a 2**20
1277 c
1278 c modifies :
1279 c ----------
1280 c nbsomm : nombre de sommets apres identification
1281 c pxyd   : tableau des coordonnees 2d des points
1282 c          par point : x  y  distance_souhaitee
1283 c letree : arbre-4 des triangles equilateraux (te) fond de la triangulation
1284 c          letree(0,0) : no du 1-er te vide dans letree
1285 c          letree(1,0) : maximum du 1-er indice de letree (ici 8)
1286 c          letree(2,0) : maximum declare du 2-eme indice de letree (ici mxtree)
1287 c          letree(0:8,1) : racine de l'arbre  (triangle sans sur triangle)
1288 c          si letree(0,.)>0 alors
1289 c             letree(0:3,j) : no (>0) letree des 4 sous-triangles du triangle j
1290 c          sinon
1291 c             letree(0:3,j) :-no pxyd des 1 a 4 points internes au triangle j
1292 c                             0  si pas de point
1293 c                             ( j est alors une feuille de l'arbre )
1294 c          letree(4,j) : no letree du sur-triangle du triangle j
1295 c          letree(5,j) : 0 1 2 3 no du sous-triangle j pour son sur-triangle
1296 c          letree(6:8,j) : no pxyd des 3 sommets du triangle j
1297 c
1298 c auxiliaire :
1299 c ------------
1300 c laqueu : mxqueu entiers servant de queue pour le parcours de letree
1301 c
1302 c sorties:
1303 c --------
1304 c ierr   :  0 si pas d'erreur
1305 c          51 si saturation letree dans te4ste
1306 c          52 si saturation pxyd   dans te4ste
1307 c          >0 si autre erreur
1308 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1309 c auteur : alain perronnet  analyse numerique paris upmc      avril 1997
1310 c2345x7..............................................................012
1311       double precision  ampli
1312       parameter        (ampli=1.34d0)
1313       common / unites / lecteu, imprim, intera, nunite(29)
1314 c
1315       double precision  pxyd(3,mxsomm), d2, aretm2
1316       double precision  comxmi(3,2),aretmx,a,s,xrmin,xrmax,yrmin,yrmax
1317       double precision  dmin, dmax
1318       integer           letree(0:8,0:*)
1319 c
1320       integer           laqueu(1:mxqueu),lequeu
1321 c                       lequeu : entree dans la queue
1322 c                       lhqueu : longueur de la queue
1323 c                       gestion circulaire
1324 c
1325       integer           nuste(3)
1326       equivalence      (nuste(1),ns1),(nuste(2),ns2),(nuste(3),ns3)
1327 c
1328 c     existence ou non de la fonction 'taille_ideale' des aretes
1329 c     autour du point.  ici la carte est supposee isotrope
1330 c     ==========================================================
1331 c     attention: si la fonction taille_ideale existe
1332 c                alors pxyd(3,*) est la taille_ideale dans l'espace initial
1333 c                sinon pxyd(3,*) est la distance calculee dans le plan par
1334 c                propagation a partir des tailles des aretes de la frontiere
1335 c
1336       if( nutysu .gt. 0 ) then
1337 c
1338 c        la fonction taille_ideale(x,y,z) existe
1339 c        ---------------------------------------
1340 c        initialisation de la distance souhaitee autour des points 1 a nbsomm
1341          do 1 i=1,nbsomm
1342 c           calcul de pxyzd(3,i)
1343             call tetaid( nutysu, pxyd(1,i), pxyd(2,i),
1344      %                   pxyd(3,i), ierr )
1345             if( ierr .ne. 0 ) goto 9999
1346  1       continue
1347 c
1348       else
1349 c
1350 c        la fonction taille_ideale(x,y,z) n'existe pas
1351 c        ---------------------------------------------
1352 c        prise en compte des distances souhaitees dans le plan
1353 c        autour des points frontaliers et des points internes imposes
1354 c        toutes les autres distances souhaitees ont ete mis a aretmx
1355 c        lors de l'execution du sp teqini
1356          do 3 i=1,nbarpi
1357 c           le sommet i n'est pas un sommet de letree => sommet frontalier
1358 c           recherche du sous-triangle minimal feuille contenant le point i
1359             nte = 1
1360  2          nte = notrpt( pxyd(1,i), pxyd, nte, letree )
1361 c           la distance au sommet le plus eloigne est elle inferieure
1362 c           a la distance souhaitee?
1363             ns1 = letree(6,nte)
1364             ns2 = letree(7,nte)
1365             ns3 = letree(8,nte)
1366             d2  = max( ( pxyd(1,i)-pxyd(1,ns1) )**2 +
1367      %                 ( pxyd(2,i)-pxyd(2,ns1) )**2
1368      %               , ( pxyd(1,i)-pxyd(1,ns2) )**2 +
1369      %                 ( pxyd(2,i)-pxyd(2,ns2) )**2
1370      %               , ( pxyd(1,i)-pxyd(1,ns3) )**2 +
1371      %                 ( pxyd(2,i)-pxyd(2,ns3) )**2 )
1372             if( d2 .gt. pxyd(3,i)**2 ) then
1373 c              le triangle nte trop grand doit etre subdivise en 4 sous-triangle
1374                call te4ste( nbsomm, mxsomm, pxyd, nte, letree,
1375      &                      ierr )
1376                if( ierr .ne. 0 ) return
1377                goto 2
1378             endif
1379  3       continue
1380       endif
1381 c
1382 c     le sous-triangle central de la racine est decoupe systematiquement
1383 c     ==================================================================
1384       nte = 2
1385       if( letree(0,2) .le. 0 ) then
1386 c        le sous-triangle central de la racine n'est pas subdivise
1387 c        il est donc decoupe en 4 soustriangles
1388          nbsom0 = nbsomm
1389          call te4ste( nbsomm, mxsomm, pxyd, nte, letree,
1390      %                ierr )
1391          if( ierr .ne. 0 ) return
1392          do 4 i=nbsom0+1,nbsomm
1393 c           mise a jour de taille_ideale des nouveaux sommets de te
1394             call tetaid( nutysu, pxyd(1,i), pxyd(2,i), pxyd(3,i), ierr )
1395             if( ierr .ne. 0 ) goto 9999
1396  4       continue
1397       endif
1398 c
1399 c     le carre de la longueur de l'arete de triangles equilateraux
1400 c     souhaitee pour le fond de la triangulation
1401       aretm2 = (aretmx*ampli) ** 2
1402 c
1403 c     tout te contenu dans le rectangle englobant doit avoir un
1404 c     cote < aretmx et etre de meme taille que les te voisins
1405 c     s'il contient un point; sinon un seul saut de taille est permis
1406 c     ===============================================================
1407 c     le rectangle englobant pour selectionner les te "internes"
1408 c     le numero des 3 sommets du te englobant racine de l'arbre des te
1409       ns1 = letree(6,1)
1410       ns2 = letree(7,1)
1411       ns3 = letree(8,1)
1412       a   = aretmx * 0.01d0
1413 c     abscisse du milieu de l'arete gauche du te 1
1414       s      = ( pxyd(1,ns1) + pxyd(1,ns3) ) / 2
1415       xrmin  = min( s, comxmi(1,1) - aretmx ) - a
1416 c     abscisse du milieu de l'arete droite du te 1
1417       s      = ( pxyd(1,ns2) + pxyd(1,ns3) ) / 2
1418       xrmax  = max( s, comxmi(1,2) + aretmx ) + a
1419       yrmin  = comxmi(2,1) - aretmx
1420 c     ordonnee de la droite passant par les milieus des 2 aretes
1421 c     droite gauche du te 1
1422       s      = ( pxyd(2,ns1) + pxyd(2,ns3) ) / 2
1423       yrmax  = max( s, comxmi(2,2) + aretmx ) + a
1424 c
1425 c     cas particulier de 3 ou 4 ou peu d'aretes frontalieres
1426       if( nbarpi .le. 8 ) then
1427 c        tout le triangle englobant (racine) est a prendre en compte
1428          xrmin = pxyd(1,ns1) - a
1429          xrmax = pxyd(1,ns2) + a
1430          yrmin = pxyd(2,ns1) - a
1431          yrmax = pxyd(2,ns3) + a
1432       endif
1433 c
1434       nbs0   = nbsomm
1435       nbiter = -1
1436 c
1437 c     initialisation de la queue
1438   5   nbiter = nbiter + 1
1439       lequeu = 1
1440       lhqueu = 0
1441 c     la racine de letree initialise la queue
1442       laqueu(1) = 1
1443 c
1444 c     tant que la longueur de la queue est >=0 traiter le debut de queue
1445  10   if( lhqueu .ge. 0 ) then
1446 c
1447 c        le triangle te a traiter
1448          i   = lequeu - lhqueu
1449          if( i .le. 0 ) i = mxqueu + i
1450          nte = laqueu( i )
1451 c        la longueur de la queue est reduite
1452          lhqueu = lhqueu - 1
1453 c
1454 c        nte est il un sous-triangle feuille minimal ?
1455  15      if( letree(0,nte) .gt. 0 ) then
1456 c
1457 c           non les 4 sous-triangles sont mis dans la queue
1458             if( lhqueu + 4 .ge. mxqueu ) then
1459                write(imprim,*) 'tehote: saturation de la queue'
1460                ierr = 7
1461                return
1462             endif
1463             do 20 i=3,0,-1
1464 c              ajout du sous-triangle i
1465                lhqueu = lhqueu + 1
1466                lequeu = lequeu + 1
1467                if( lequeu .gt. mxqueu ) lequeu = lequeu - mxqueu
1468                laqueu( lequeu ) = letree( i, nte )
1469  20         continue
1470             goto 10
1471 c
1472          endif
1473 c
1474 c        ici nte est un triangle minimal non subdivise
1475 c        ---------------------------------------------
1476 c        le te est il dans le cadre englobant de l'objet ?
1477          ns1 = letree(6,nte)
1478          ns2 = letree(7,nte)
1479          ns3 = letree(8,nte)
1480          if( pxyd(1,ns1) .gt. pxyd(1,ns2) ) then
1481             dmin = pxyd(1,ns2)
1482             dmax = pxyd(1,ns1)
1483          else
1484             dmin = pxyd(1,ns1)
1485             dmax = pxyd(1,ns2)
1486          endif
1487          if( (xrmin .le. dmin .and. dmin .le. xrmax) .or.
1488      %       (xrmin .le. dmax .and. dmax .le. xrmax) ) then
1489             if( pxyd(2,ns1) .gt. pxyd(2,ns3) ) then
1490                dmin = pxyd(2,ns3)
1491                dmax = pxyd(2,ns1)
1492             else
1493                dmin = pxyd(2,ns1)
1494                dmax = pxyd(2,ns3)
1495             endif
1496             if( (yrmin .le. dmin .and. dmin .le. yrmax) .or.
1497      %          (yrmin .le. dmax .and. dmax .le. yrmax) ) then
1498 c
1499 c              nte est un te feuille et interne au rectangle englobant
1500 c              =======================================================
1501 c              le carre de la longueur de l'arete du te de numero nte
1502                d2 = (pxyd(1,ns1)-pxyd(1,ns2)) ** 2 +
1503      %              (pxyd(2,ns1)-pxyd(2,ns2)) ** 2
1504 c
1505                if( nutysu .eq. 0 ) then
1506 c
1507 c                 il n'existe pas de fonction 'taille_ideale'
1508 c                 -------------------------------------------
1509 c                 si la taille effective de l'arete du te est superieure a aretmx
1510 c                 alors le te est decoupe
1511                   if( d2 .gt. aretm2 ) then
1512 c                    le triangle nte trop grand doit etre subdivise
1513 c                    en 4 sous-triangles
1514                      call te4ste( nbsomm,mxsomm, pxyd,
1515      %                            nte, letree, ierr )
1516                      if( ierr .ne. 0 ) return
1517                      goto 15
1518                   endif
1519 c
1520                else
1521 c
1522 c                 il existe ici une fonction 'taille_ideale'
1523 c                 ------------------------------------------
1524 c                 si la taille effective de l'arete du te est superieure au mini
1525 c                 des 3 tailles_ideales aux sommets  alors le te est decoupe
1526                   do 28 i=1,3
1527                      if( d2 .gt. (pxyd(3,nuste(i))*ampli)**2 ) then
1528 c                       le triangle nte trop grand doit etre subdivise
1529 c                       en 4 sous-triangles
1530                         nbsom0 = nbsomm
1531                         call te4ste( nbsomm, mxsomm, pxyd,
1532      &                               nte, letree, ierr )
1533                         if( ierr .ne. 0 ) return
1534                         do 27 j=nbsom0+1,nbsomm
1535 c                          mise a jour de taille_ideale des nouveaux sommets de
1536                            call tetaid( nutysu, pxyd(1,j), pxyd(2,j),
1537      %                                  pxyd(3,j), ierr )
1538                            if( ierr .ne. 0 ) goto 9999
1539  27                     continue
1540                         goto 15
1541                      endif
1542  28               continue
1543                endif
1544 c
1545 c              recherche du nombre de niveaux entre nte et les te voisins par se
1546 c              si la difference de subdivisions excede 1 alors le plus grand des
1547 c              =================================================================
1548  29            do 30 i=1,3
1549 c
1550 c                 noteva triangle voisin de nte par l'arete i
1551                   call n1trva( nte, i, letree, noteva, niveau )
1552                   if( noteva .le. 0 ) goto 30
1553 c                 il existe un te voisin
1554                   if( niveau .gt. 0 ) goto 30
1555 c                 nte a un te voisin plus petit ou egal
1556                   if( letree(0,noteva) .le. 0 ) goto 30
1557 c                 nte a un te voisin noteva subdivise au moins une fois
1558 c
1559                   if( nbiter .gt. 0 ) then
1560 c                    les 2 sous triangles voisins sont-ils subdivises?
1561                      ns2 = letree(i,noteva)
1562                      if( letree(0,ns2) .le. 0 ) then
1563 c                       ns2 n'est pas subdivise
1564                         ns2 = letree(nosui3(i),noteva)
1565                         if( letree(0,ns2) .le. 0 ) then
1566 c                          les 2 sous-triangles ne sont pas subdivises
1567                            goto 30
1568                         endif
1569                      endif
1570                   endif
1571 c
1572 c                 saut>1 => le triangle nte doit etre subdivise en 4 sous-triang
1573 c                 --------------------------------------------------------------
1574                   nbsom0 = nbsomm
1575                   call te4ste( nbsomm,mxsomm, pxyd, nte, letree,
1576      &                         ierr )
1577                   if( ierr .ne. 0 ) return
1578                   if( nutysu .gt. 0 ) then
1579                      do 32 j=nbsom0+1,nbsomm
1580 c                       mise a jour de taille_ideale des nouveaux sommets de te
1581                         call tetaid( nutysu, pxyd(1,j), pxyd(2,j),
1582      %                               pxyd(3,j), ierr )
1583                         if( ierr .ne. 0 ) goto 9999
1584  32                  continue
1585                   endif
1586                   goto 15
1587 c
1588  30            continue
1589             endif
1590          endif
1591          goto 10
1592       endif
1593       if( nbs0 .lt. nbsomm ) then
1594          nbs0 = nbsomm
1595          goto 5
1596       endif
1597       return
1598 c
1599 c     pb dans le calcul de la fonction taille_ideale
1600
1601  9999 write(imprim,*) 'pb dans le calcul de taille_ideale'
1602 c      nblgrc(nrerr) = 1
1603 c      kerr(1) = 'pb dans le calcul de taille_ideale'
1604 c      call lereur
1605       return
1606       end
1607
1608
1609       subroutine tetrte( comxmi, aretmx, nbarpi, mxsomm, pxyd,
1610      %                   mxqueu, laqueu, letree,
1611      %                   mosoar, mxsoar, n1soar, nosoar,
1612      %                   moartr, mxartr, n1artr, noartr, noarst,
1613      %                   ierr  )
1614 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1615 c but :    trianguler les triangles equilateraux feuilles et
1616 c -----    les points de la frontiere et les points internes imposes
1617 c
1618 c attention: la triangulation finale n'est pas de type delaunay!
1619 c
1620 c entrees:
1621 c --------
1622 c comxmi : minimum et maximum des coordonnees de l'objet
1623 c aretmx : longueur maximale des aretes des triangles equilateraux
1624 c nbarpi : nombre de sommets de la frontiere + nombre de points internes
1625 c          imposes par l'utilisateur
1626 c mxsomm : nombre maximal de sommets declarables dans pxyd
1627 c pxyd   : tableau des coordonnees 2d des points
1628 c          par point : x  y  distance_souhaitee
1629 c
1630 c mxqueu : nombre d'entiers utilisables dans laqueu
1631 c mosoar : nombre maximal d'entiers par arete du tableau nosoar
1632 c mxsoar : nombre maximal d'aretes stockables dans le tableau nosoar
1633 c moartr : nombre maximal d'entiers par arete du tableau noartr
1634 c mxartr : nombre maximal de triangles stockables dans le tableau noartr
1635 c letree : arbre-4 des triangles equilateraux (te) fond de la triangulation
1636 c          letree(0,0) : no du 1-er te vide dans letree
1637 c          letree(0,1) : maximum du 1-er indice de letree (ici 8)
1638 c          letree(0,2) : maximum declare du 2-eme indice de letree (ici mxtree)
1639 c          letree(0:8,1) : racine de l'arbre  (triangle sans sur triangle)
1640 c          si letree(0,.)>0 alors
1641 c             letree(0:3,j) : no (>0) letree des 4 sous-triangles du triangle j
1642 c          sinon
1643 c             letree(0:3,j) :-no pxyd des 1 a 4 points internes au triangle j
1644 c                             0  si pas de point
1645 c                             ( j est alors une feuille de l'arbre )
1646 c          letree(4,j) : no letree du sur-triangle du triangle j
1647 c          letree(5,j) : 0 1 2 3 no du sous-triangle j pour son sur-triangle
1648 c          letree(6:8,j) : no pxyd des 3 sommets du triangle j
1649 c
1650 c modifies:
1651 c ---------
1652 c n1soar : numero de la premiere arete vide dans le tableau nosoar
1653 c          une arete i de nosoar est vide  <=>  nosoar(1,i)=0
1654 c nosoar : numero des 2 sommets , no ligne, 2 triangles de l'arete,
1655 c          chainage des aretes frontalieres, chainage du hachage des aretes
1656 c          hachage des aretes = nosoar(1)+nosoar(2)*2
1657 c noarst : noarst(i) numero d'une arete de sommet i
1658 c
1659 c auxiliaire :
1660 c ------------
1661 c laqueu : mxqueu entiers servant de queue pour le parcours de letree
1662 c
1663 c sorties:
1664 c --------
1665 c n1artr : numero du premier triangle vide dans le tableau noartr
1666 c          le chainage des triangles vides se fait sur noartr(2,.)
1667 c noartr : les 3 aretes des triangles +-arete1, +-arete2, +-arete3
1668 c          arete1 = 0 si triangle vide => arete2 = triangle vide suivant
1669 c ierr   : =0 si pas d'erreur
1670 c          =1 si le tableau nosoar est sature
1671 c          =2 si le tableau noartr est sature
1672 c          =3 si aucun des triangles ne contient l'un des points internes d'un t
1673 c          =5 si saturation de la queue de parcours de l'arbre des te
1674 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1675 c auteur : alain perronnet  analyse numerique paris upmc       mars 1997
1676 c2345x7..............................................................012
1677       common / unites / lecteu, imprim, intera, nunite(29)
1678 c
1679       double precision  pxyd(3,mxsomm)
1680       double precision  comxmi(3,2),aretmx,a,s,xrmin,xrmax,yrmin,yrmax
1681       double precision  dmin, dmax
1682 c
1683       integer           nosoar(mosoar,mxsoar),
1684      %                  noartr(moartr,mxartr),
1685      %                  noarst(mxsomm)
1686 c
1687       integer           letree(0:8,0:*)
1688       integer           laqueu(1:mxqueu)
1689 c                       lequeu:entree dans la queue en gestion circulaire
1690 c                       lhqueu:longueur de la queue en gestion circulaire
1691 c
1692       integer           milieu(3), nutr(1:13)
1693 c
1694 c     le rectangle englobant pour selectionner les te "internes"
1695 c     le numero des 3 sommets du te englobant racine de l'arbre des te
1696       ns1 = letree(6,1)
1697       ns2 = letree(7,1)
1698       ns3 = letree(8,1)
1699       a   = aretmx * 0.01d0
1700 c     abscisse du milieu de l'arete gauche du te 1
1701       s      = ( pxyd(1,ns1) + pxyd(1,ns3) ) / 2
1702       xrmin  = min( s, comxmi(1,1) - aretmx ) - a
1703 c     abscisse du milieu de l'arete droite du te 1
1704       s      = ( pxyd(1,ns2) + pxyd(1,ns3) ) / 2
1705       xrmax  = max( s, comxmi(1,2) + aretmx ) + a
1706       yrmin  = comxmi(2,1) - aretmx
1707 c     ordonnee de la droite passant par les milieus des 2 aretes
1708 c     droite gauche du te 1
1709       s      = ( pxyd(2,ns1) + pxyd(2,ns3) ) / 2
1710       yrmax  = max( s, comxmi(2,2) + aretmx ) + a
1711 c
1712 c     cas particulier de 3 ou 4 ou peu d'aretes frontalieres
1713       if( nbarpi .le. 8 ) then
1714 c        tout le triangle englobant (racine) est a prendre en compte
1715          xrmin = pxyd(1,ns1) - a
1716          xrmax = pxyd(1,ns2) + a
1717          yrmin = pxyd(2,ns1) - a
1718          yrmax = pxyd(2,ns3) + a
1719       endif
1720 c
1721 c     initialisation du tableau noartr
1722       do 5 i=1,mxartr
1723 c        le numero de l'arete est inconnu
1724          noartr(1,i) = 0
1725 c        le chainage sur le triangle vide suivant
1726          noartr(2,i) = i+1
1727  5    continue
1728       noartr(2,mxartr) = 0
1729       n1artr = 1
1730 c
1731 c     parcours des te jusqu'a trianguler toutes les feuilles (triangles eq)
1732 c     =====================================================================
1733 c     initialisation de la queue sur les te
1734       ierr   = 0
1735       lequeu = 1
1736       lhqueu = 0
1737 c     la racine de letree initialise la queue
1738       laqueu(1) = 1
1739 c
1740 c     tant que la longueur de la queue est >=0 traiter le debut de queue
1741  10   if( lhqueu .ge. 0 ) then
1742 c
1743 c        le triangle te a traiter
1744          i   = lequeu - lhqueu
1745          if( i .le. 0 ) i = mxqueu + i
1746          nte = laqueu( i )
1747 c        la longueur est reduite
1748          lhqueu = lhqueu - 1
1749 c
1750 c        nte est il un sous-triangle feuille (minimal) ?
1751  15      if( letree(0,nte) .gt. 0 ) then
1752 c           non les 4 sous-triangles sont mis dans la queue
1753             if( lhqueu + 4 .ge. mxqueu ) then
1754                write(imprim,*) 'tetrte: saturation de la queue'
1755                ierr = 5
1756                return
1757             endif
1758             do 20 i=3,0,-1
1759 c              ajout du sous-triangle i
1760                lhqueu = lhqueu + 1
1761                lequeu = lequeu + 1
1762                if( lequeu .gt. mxqueu ) lequeu = lequeu - mxqueu
1763                laqueu( lequeu ) = letree( i, nte )
1764  20         continue
1765             goto 10
1766          endif
1767 c
1768 c        ici nte est un triangle minimal non subdivise
1769 c        ---------------------------------------------
1770 c        le te est il dans le cadre englobant de l'objet ?
1771          ns1 = letree(6,nte)
1772          ns2 = letree(7,nte)
1773          ns3 = letree(8,nte)
1774          if( pxyd(1,ns1) .gt. pxyd(1,ns2) ) then
1775             dmin = pxyd(1,ns2)
1776             dmax = pxyd(1,ns1)
1777          else
1778             dmin = pxyd(1,ns1)
1779             dmax = pxyd(1,ns2)
1780          endif
1781          if( (xrmin .le. dmin .and. dmin .le. xrmax) .or.
1782      %       (xrmin .le. dmax .and. dmax .le. xrmax) ) then
1783             if( pxyd(2,ns1) .gt. pxyd(2,ns3) ) then
1784                dmin = pxyd(2,ns3)
1785                dmax = pxyd(2,ns1)
1786             else
1787                dmin = pxyd(2,ns1)
1788                dmax = pxyd(2,ns3)
1789             endif
1790             if( (yrmin .le. dmin .and. dmin .le. yrmax) .or.
1791      %          (yrmin .le. dmax .and. dmax .le. yrmax) ) then
1792 c
1793 c              te minimal et interne au rectangle englobant
1794 c              --------------------------------------------
1795 c              recherche du nombre de niveaux entre nte et les te voisins
1796 c              par ses aretes
1797                nbmili = 0
1798                do 30 i=1,3
1799 c
1800 c                 a priori pas de milieu de l'arete i du te nte
1801                   milieu(i) = 0
1802 c
1803 c                 recherche de noteva te voisin de nte par l'arete i
1804                   call n1trva( nte, i, letree, noteva, niveau )
1805 c                 noteva  : >0 numero letree du te voisin par l'arete i
1806 c                           =0 si pas de te voisin (racine , ... )
1807 c                 niveau  : =0 si nte et noteva ont meme taille
1808 c                           >0 nte est 4**niveau fois plus petit que noteva
1809                   if( noteva .gt. 0 ) then
1810 c                    il existe un te voisin
1811                      if( letree(0,noteva) .gt. 0 ) then
1812 c                       noteva est plus petit que nte
1813 c                       => recherche du numero du milieu du cote=sommet du te no
1814 c                       le sous-te 0 du te noteva
1815                         nsot = letree(0,noteva)
1816 c                       le numero dans pxyd du milieu de l'arete i de nte
1817                         milieu( i ) = letree( 5+nopre3(i), nsot )
1818                         nbmili = nbmili + 1
1819                      endif
1820                   endif
1821 c
1822  30            continue
1823 c
1824 c              triangulation du te nte en fonction du nombre de ses milieux
1825                goto( 50, 100, 200, 300 ) , nbmili + 1
1826 c
1827 c              0 milieu => 1 triangle = le te nte
1828 c              ----------------------------------
1829  50            call f0trte( letree(0,nte),  pxyd,
1830      %                      mosoar, mxsoar, n1soar, nosoar,
1831      %                      moartr, mxartr, n1artr, noartr,
1832      %                      noarst,
1833      %                      nbtr,   nutr,   ierr )
1834                if( ierr .ne. 0 ) return
1835                goto 10
1836 c
1837 c              1 milieu => 2 triangles = 2 demi te
1838 c              -----------------------------------
1839  100           call f1trte( letree(0,nte),  pxyd,   milieu,
1840      %                      mosoar, mxsoar, n1soar, nosoar,
1841      %                      moartr, mxartr, n1artr, noartr,
1842      %                      noarst,
1843      %                      nbtr,   nutr,   ierr )
1844                if( ierr .ne. 0 ) return
1845                goto 10
1846 c
1847 c              2 milieux => 3 triangles
1848 c              -----------------------------------
1849  200           call f2trte( letree(0,nte),  pxyd,   milieu,
1850      %                      mosoar, mxsoar, n1soar, nosoar,
1851      %                      moartr, mxartr, n1artr, noartr,
1852      %                      noarst,
1853      %                      nbtr,   nutr,   ierr )
1854                if( ierr .ne. 0 ) return
1855                goto 10
1856 c
1857 c              3 milieux => 4 triangles = 4 quart te
1858 c              -------------------------------------
1859  300           call f3trte( letree(0,nte),  pxyd,   milieu,
1860      %                      mosoar, mxsoar, n1soar, nosoar,
1861      %                      moartr, mxartr, n1artr, noartr,
1862      %                      noarst,
1863      %                      nbtr,   nutr,   ierr )
1864                if( ierr .ne. 0 ) return
1865                goto 10
1866             endif
1867          endif
1868          goto 10
1869       endif
1870       end
1871
1872
1873       subroutine aisoar( mosoar, mxsoar, nosoar, na1 )
1874 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1875 c but :    chainer en colonne lchain les aretes non vides et
1876 c -----    non frontalieres du tableau nosoar
1877 c
1878 c entrees:
1879 c --------
1880 c mosoar : nombre maximal d'entiers par arete dans le tableau nosoar
1881 c mxsoar : nombre maximal d'aretes frontalieres declarables
1882 c
1883 c modifies :
1884 c ----------
1885 c nosoar : numero des 2 sommets , no ligne, 2 triangles, chainages en +
1886 c          nosoar(lchain,i)=arete interne suivante
1887 c
1888 c sortie :
1889 c --------
1890 c na1    : numero dans nosoar de la premiere arete interne
1891 c          les suivantes sont nosoar(lchain,na1), ...
1892 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1893 c auteur : alain perronnet  analyse numerique paris upmc       mars 1997
1894 c....................................................................012
1895       parameter (lchain=6)
1896       integer    nosoar(mosoar,mxsoar)
1897 c
1898 c     formation du chainage des aretes internes a echanger eventuellement
1899 c     recherche de la premiere arete non vide et non frontaliere
1900       do 10 na1=1,mxsoar
1901          if( nosoar(1,na1) .gt. 0 .and. nosoar(3,na1) .le. 0 ) goto 15
1902  10   continue
1903 c
1904 c     protection de la premiere arete non vide et non frontaliere
1905  15   na0 = na1
1906       do 20 na=na1+1,mxsoar
1907          if( nosoar(1,na) .gt. 0 .and. nosoar(3,na) .le. 0 ) then
1908 c           arete interne => elle est chainee a partir de la precedente
1909             nosoar(lchain,na0) = na
1910             na0 = na
1911          endif
1912  20   continue
1913 c
1914 c     la derniere arete interne n'a pas de suivante
1915       nosoar(lchain,na0) = 0
1916       end
1917
1918
1919       subroutine tedela( pxyd,   noarst,
1920      %                   mosoar, mxsoar, n1soar, nosoar, n1ardv,
1921      %                   moartr, mxartr, n1artr, noartr, modifs )
1922 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1923 c but :    pour toutes les aretes chainees dans nosoar(lchain,*)
1924 c -----    du tableau nosoar
1925 c          echanger la diagonale des 2 triangles si le sommet oppose
1926 c          a un triangle ayant en commun une arete appartient au cercle
1927 c          circonscrit de l'autre (violation boule vide delaunay)
1928 c
1929 c entrees:
1930 c --------
1931 c pxyd   : tableau des x  y  distance_souhaitee de chaque sommet
1932 c
1933 c modifies :
1934 c ----------
1935 c noarst : noarst(i) numero d'une arete de sommet i
1936 c mosoar : nombre maximal d'entiers par arete dans le tableau nosoar
1937 c mxsoar : nombre maximal d'aretes frontalieres declarables
1938 c n1soar : numero de la premiere arete vide dans le tableau nosoar
1939 c nosoar : numero des 2 sommets , no ligne, 2 triangles, chainages en +
1940 c n1ardv : numero dans nosoar de la premiere arete du chainage
1941 c          des aretes a rendre delaunay
1942 c
1943 c moartr : nombre d'entiers par triangle dans le tableau noartr
1944 c mxartr : nombre maximal de triangles declarables dans noartr
1945 c n1artr : numero du premier triangle vide dans le tableau noartr
1946 c          le chainage des triangles vides se fait sur noartr(2,.)
1947 c noartr : les 3 aretes des triangles +-arete1, +-arete2, +-arete3
1948 c          arete1 = 0 si triangle vide => arete2 = triangle vide suivant
1949 c modifs : nombre d'echanges de diagonales pour maximiser la qualite
1950 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1951 c auteur : alain perronnet  analyse numerique paris upmc       mars 1997
1952 c....................................................................012
1953       parameter        (lchain=6)
1954       common / unites / lecteu, imprim, nunite(30)
1955       double precision  pxyd(3,*), surtd2, s123, s142, s143, s234,
1956      %                  s12, s34, a12, cetria(3), r0
1957       integer           nosoar(mosoar,mxsoar),
1958      %                  noartr(moartr,mxartr),
1959      %                  noarst(*)
1960 c
1961 c     le nombre d'echanges de diagonales pour minimiser l'aire
1962       modifs = 0
1963       r0     = 0
1964 c
1965 c     la premiere arete du chainage des aretes a rendre delaunay
1966       na0 = n1ardv
1967 c
1968 c     tant que la pile des aretes a echanger eventuellement est non vide
1969 c     ==================================================================
1970  20   if( na0 .gt. 0 ) then
1971 c
1972 c        l'arete a traiter
1973          na  = na0
1974 c        la prochaine arete a traiter
1975          na0 = nosoar(lchain,na0)
1976 c
1977 c        l'arete est marquee traitee avec le numero -1
1978          nosoar(lchain,na) = -1
1979 c
1980 c        l'arete est elle active?
1981          if( nosoar(1,na) .eq. 0 ) goto 20
1982 c
1983 c        si arete frontaliere pas d'echange possible
1984          if( nosoar(3,na) .gt. 0 ) goto 20
1985 c
1986 c        existe-t-il 2 triangles ayant cette arete commune?
1987          if( nosoar(4,na) .le. 0 .or. nosoar(5,na) .le. 0 ) goto 20
1988 c
1989 c        aucun des 2 triangles est-il desactive?
1990          if( noartr(1,nosoar(4,na)) .eq. 0 .or.
1991      %       noartr(1,nosoar(5,na)) .eq. 0 ) goto 20
1992 c
1993 c        l'arete appartient a deux triangles actifs
1994 c        le numero des 4 sommets du quadrangle des 2 triangles
1995          call mt4sqa( na, moartr, noartr, mosoar, nosoar,
1996      %                ns1, ns2, ns3, ns4 )
1997          if( ns4 .eq. 0 ) goto 20
1998 c
1999 c        carre de la longueur de l'arete ns1 ns2
2000          a12 = (pxyd(1,ns2)-pxyd(1,ns1))**2+(pxyd(2,ns2)-pxyd(2,ns1))**2
2001 c
2002 c        comparaison de la somme des aires des 2 triangles
2003 c        -------------------------------------------------
2004 c        calcul des surfaces des triangles 123 et 142 de cette arete
2005          s123=surtd2( pxyd(1,ns1), pxyd(1,ns2), pxyd(1,ns3) )
2006          s142=surtd2( pxyd(1,ns1), pxyd(1,ns4), pxyd(1,ns2) )
2007          s12 = abs( s123 ) + abs( s142 )
2008          if( s12 .le. 0.001*a12 ) goto 20
2009 c
2010 c        calcul des surfaces des triangles 143 et 234 de cette arete
2011          s143=surtd2( pxyd(1,ns1), pxyd(1,ns4), pxyd(1,ns3) )
2012          s234=surtd2( pxyd(1,ns2), pxyd(1,ns3), pxyd(1,ns4) )
2013          s34 = abs( s234 ) + abs( s143 )
2014 c
2015          if( abs(s34-s12) .gt. 1d-15*s34 ) goto 20
2016 c
2017 c        quadrangle convexe : le critere de delaunay intervient
2018 c        ------------------   ---------------------------------
2019 c        calcul du centre et rayon de la boule circonscrite a 123
2020 c        pas d'affichage si le triangle est degenere
2021          ierr = -1
2022          call cenced( pxyd(1,ns1), pxyd(1,ns2), pxyd(1,ns3), cetria,
2023      %                ierr )
2024          if( ierr .gt. 0 ) then
2025 c           ierr=1 si triangle degenere  => abandon
2026             goto 20
2027          endif
2028 c
2029          if( (cetria(1)-pxyd(1,ns4))**2+(cetria(2)-pxyd(2,ns4))**2
2030      %       .lt. cetria(3) ) then
2031 c
2032 c           protection contre une boucle infinie sur le meme cercle
2033             if( r0 .eq. cetria(3) ) goto 20
2034 c
2035 c           oui: ns4 est dans le cercle circonscrit a ns1 ns2 ns3
2036 c           => ns3 est aussi dans le cercle circonscrit de ns1 ns2 ns4
2037 c
2038 cccc           les 2 triangles d'arete na sont effaces
2039 ccc            do 25 j=4,5
2040 ccc               nt = nosoar(j,na)
2041 cccc              trace du triangle nt
2042 ccc               call mttrtr( pxyd, nt, moartr, noartr, mosoar, nosoar,
2043 ccc     %                      ncnoir, ncjaun )
2044 ccc 25         continue
2045 c
2046 c           echange de la diagonale 12 par 34 des 2 triangles
2047             call te2t2t( na,     mosoar, n1soar, nosoar, noarst,
2048      %                   moartr, noartr, na34 )
2049             if( na34 .eq. 0 ) goto 20
2050             r0 = cetria(3)
2051 c
2052 c           l'arete na34 est marquee traitee
2053             nosoar(lchain,na34) = -1
2054             modifs = modifs + 1
2055 c
2056 c           les aretes internes peripheriques des 2 triangles sont enchainees
2057             do 60 j=4,5
2058                nt = nosoar(j,na34)
2059 cccc              trace du triangle nt
2060 ccc               call mttrtr( pxyd, nt, moartr, noartr, mosoar, nosoar,
2061 ccc     %                      ncoran, ncgric )
2062                do 50 i=1,3
2063                   n = abs( noartr(i,nt) )
2064                   if( n .ne. na34 ) then
2065                      if( nosoar(3,n)      .eq.  0  .and.
2066      %                   nosoar(lchain,n) .eq. -1 ) then
2067 c                        cette arete marquee est chainee pour etre traitee
2068                          nosoar(lchain,n) = na0
2069                          na0 = n
2070                      endif
2071                   endif
2072  50            continue
2073  60         continue
2074             goto 20
2075          endif
2076 c
2077 c        retour en haut de la pile des aretes a traiter
2078          goto 20
2079       endif
2080       end
2081
2082
2083       subroutine terefr( nbarpi, pxyd,
2084      %                   mosoar, mxsoar, n1soar, nosoar,
2085      %                   moartr, n1artr, noartr, noarst,
2086      %                   mxarcf, n1arcf, noarcf, larmin, notrcf,
2087      %                   nbarpe, ierr )
2088 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2089 c but :   recherche des aretes de la frontiere non dans la triangulation
2090 c -----   triangulation frontale pour les reobtenir
2091 c
2092 c         attention: le chainage lchain de nosoar devient celui des cf
2093 c
2094 c entrees:
2095 c --------
2096 c          le tableau nosoar
2097 c nbarpi : numero du dernier point interne impose par l'utilisateur
2098 c pxyd   : tableau des coordonnees 2d des points
2099 c          par point : x  y  distance_souhaitee
2100 c mosoar : nombre maximal d'entiers par arete et
2101 c          indice dans nosoar de l'arete suivante dans le hachage
2102 c mxsoar : nombre maximal d'aretes stockables dans le tableau nosoar
2103 c          attention: mxsoar>3*mxsomm obligatoire!
2104 c moartr : nombre maximal d'entiers par arete du tableau noartr
2105 c mxarcf : nombre de variables des tableaux n1arcf, noarcf, larmin, notrcf
2106 c
2107 c modifies:
2108 c ---------
2109 c n1soar : no de l'eventuelle premiere arete libre dans le tableau nosoar
2110 c          chainage des vides suivant en 3 et precedant en 2 de nosoar
2111 c nosoar : numero des 2 sommets , no ligne, 2 triangles de l'arete,
2112 c          chainage des aretes frontalieres, chainage du hachage des aretes
2113 c          hachage des aretes = nosoar(1)+nosoar(2)*2
2114 c          avec mxsoar>=3*mxsomm
2115 c          une arete i de nosoar est vide <=> nosoar(1,i)=0 et
2116 c          nosoar(2,arete vide)=l'arete vide qui precede
2117 c          nosoar(3,arete vide)=l'arete vide qui suit
2118 c n1artr : numero du premier triangle vide dans le tableau noartr
2119 c          le chainage des triangles vides se fait sur noartr(2,.)
2120 c noartr : les 3 aretes des triangles +-arete1, +-arete2, +-arete3
2121 c          arete1 = 0 si triangle vide => arete2 = triangle vide suivant
2122 c noarst : noarst(i) numero d'une arete de sommet i
2123 c
2124 c
2125 c auxiliaires :
2126 c -------------
2127 c n1arcf : tableau (0:mxarcf) auxiliaire d'entiers
2128 c noarcf : tableau (3,mxarcf) auxiliaire d'entiers
2129 c larmin : tableau (mxarcf)   auxiliaire d'entiers
2130 c notrcf : tableau (mxarcf)   auxiliaire d'entiers
2131 c
2132 c sortie :
2133 c --------
2134 c nbarpe : nombre d'aretes perdues puis retrouvees
2135 c ierr   : =0 si pas d'erreur
2136 c          >0 si une erreur est survenue
2137 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2138 c auteur : alain perronnet  analyse numerique paris upmc       mars 1997
2139 c....................................................................012
2140       parameter        (lchain=6)
2141       common / unites / lecteu,imprim,intera,nunite(29)
2142       double precision  pxyd(3,*)
2143       integer           nosoar(mosoar,mxsoar),
2144      %                  noartr(moartr,*),
2145      %                  noarst(*),
2146      %                  n1arcf(0:mxarcf),
2147      %                  noarcf(3,mxarcf),
2148      %                  larmin(mxarcf),
2149      %                  notrcf(mxarcf)
2150 c
2151 c     le nombre d'aretes de la frontiere non arete de la triangulation
2152       nbarpe = 0
2153 c
2154 c     initialisation du chainage des aretes des cf => 0 arete de cf
2155       do 10 narete=1,mxsoar
2156          nosoar( lchain, narete) = -1
2157  10   continue
2158 c
2159 c     boucle sur l'ensemble des aretes actuelles
2160 c     ==========================================
2161       do 30 narete=1,mxsoar
2162 c
2163          if( nosoar(3,narete) .gt. 0 ) then
2164 c           arete appartenant a une ligne => frontaliere
2165 c
2166             if(nosoar(4,narete) .le. 0 .or. nosoar(5,narete) .le. 0)then
2167 c              l'arete narete frontaliere n'appartient pas a 2 triangles
2168 c              => elle est perdue
2169                nbarpe = nbarpe + 1
2170 c
2171 c              le numero des 2 sommets de l'arete frontaliere perdue
2172                ns1 = nosoar( 1, narete )
2173                ns2 = nosoar( 2, narete )
2174 c               write(imprim,10000) ns1,(pxyd(j,ns1),j=1,2),
2175 c     %                             ns2,(pxyd(j,ns2),j=1,2)
2176 10000          format(' arete perdue a forcer',
2177      %               (t24,'sommet=',i6,' x=',g13.5,' y=',g13.5))
2178 c
2179 c              traitement de cette arete perdue ns1-ns2
2180                call tefoar( narete, nbarpi, pxyd,
2181      %                      mosoar, mxsoar, n1soar, nosoar,
2182      %                      moartr, n1artr, noartr, noarst,
2183      %                      mxarcf, n1arcf, noarcf, larmin, notrcf,
2184      %                      ierr )
2185                if( ierr .ne. 0 ) return
2186 c
2187 c              fin du traitement de cette arete perdue et retrouvee
2188             endif
2189          endif
2190 c
2191  30   continue
2192       end
2193
2194
2195       subroutine tesuex( nblftr, nulftr,
2196      %                   ndtri0, nbsomm, pxyd, nslign,
2197      %                   mosoar, mxsoar, nosoar,
2198      %                   moartr, mxartr, n1artr, noartr, noarst,
2199      %                   nbtria, letrsu, ierr  )
2200 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2201 c but :    supprimer du tableau noartr les triangles externes au domaine
2202 c -----    en annulant le numero de leur 1-ere arete dans noartr
2203 c          et en les chainant comme triangles vides
2204 c
2205 c entrees:
2206 c --------
2207 c nblftr : nombre de  lignes fermees definissant la surface
2208 c nulftr : numero des lignes fermees definissant la surface
2209 c ndtri0 : plus grand numero dans noartr d'un triangle
2210 c pxyd   : tableau des coordonnees 2d des points
2211 c          par point : x  y  distance_souhaitee
2212 c nslign : tableau du numero de sommet dans sa ligne pour chaque
2213 c          sommet frontalier
2214 c          numero du point dans le lexique point si interne impose
2215 c          0 si le point est interne non impose par l'utilisateur
2216 c         -1 si le sommet est externe au domaine
2217 c mosoar : nombre maximal d'entiers par arete et
2218 c          indice dans nosoar de l'arete suivante dans le hachage
2219 c mxsoar : nombre maximal d'aretes stockables dans le tableau nosoar
2220 c          attention: mxsoar>3*mxsomm obligatoire!
2221 c nosoar : numero des 2 sommets , no ligne, 2 triangles de l'arete,
2222 c          chainage des aretes frontalieres, chainage du hachage des aretes
2223 c          hachage des aretes = nosoar(1)+nosoar(2)*2
2224 c          avec mxsoar>=3*mxsomm
2225 c          une arete i de nosoar est vide <=> nosoar(1,i)=0 et
2226 c          nosoar(2,arete vide)=l'arete vide qui precede
2227 c          nosoar(3,arete vide)=l'arete vide qui suit
2228 c moartr : nombre maximal d'entiers par arete du tableau noartr
2229 c mxartr : nombre maximal de triangles declarables
2230 c n1artr : numero du premier triangle vide dans le tableau noartr
2231 c          le chainage des triangles vides se fait sur noartr(2,.)
2232 c noartr : les 3 aretes des triangles +-arete1, +-arete2, +-arete3
2233 c          arete1 = 0 si triangle vide => arete2 = triangle vide suivant
2234 c noarst : noarst(i) numero nosoar d'une arete de sommet i
2235 c
2236 c sorties:
2237 c --------
2238 c nbtria : nombre de triangles internes au domaine
2239 c letrsu : letrsu(nt)=numero du triangle interne, 0 sinon
2240 c noarst : noarst(i) numero nosoar d'une arete du sommet i (modifi'e)
2241 c ierr   : 0 si pas d'erreur, >0 sinon
2242 cc++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2243 c auteur : alain perronnet  analyse numerique paris upmc        mai 1999
2244 c2345x7..............................................................012
2245       double precision  pxyd(3,*)
2246       integer           nulftr(nblftr),nslign(nbsomm),
2247      %                  nosoar(mosoar,mxsoar),
2248      %                  noartr(moartr,mxartr),
2249      %                  noarst(*)
2250       integer           letrsu(1:ndtri0)
2251       double precision  dmin
2252 c
2253 c     les triangles sont a priori non marques
2254       do 5 nt=1,ndtri0
2255          letrsu(nt) = 0
2256  5    continue
2257 c
2258 c     les aretes sont marquees non chainees
2259       do 10 noar1=1,mxsoar
2260          nosoar(6,noar1) = -2
2261  10   continue
2262 c
2263 c     recherche du sommet de la triangulation de plus petite abscisse
2264 c     ===============================================================
2265       ntmin = 0
2266       dmin  = 1d38
2267       do 20 i=1,nbsomm
2268          if( pxyd(1,i) .lt. dmin ) then
2269 c           le nouveau minimum
2270             noar1 = noarst(i)
2271             if( noar1 .gt. 0 ) then
2272 c              le sommet appartient a une arete de triangle
2273                if( nosoar(4,noar1) .gt. 0 ) then
2274 c                 le nouveau minimum
2275                   dmin  = pxyd(1,i)
2276                   ntmin = i
2277                endif
2278             endif
2279          endif
2280  20   continue
2281 c
2282 c     une arete de sommet ntmin
2283       noar1 = noarst( ntmin )
2284 c     un triangle d'arete noar1
2285       ntmin = nosoar( 4, noar1 )
2286       if( ntmin .le. 0 ) then
2287 c         nblgrc(nrerr) = 1
2288 c         kerr(1) = 'pas de triangle d''abscisse minimale'
2289 c         call lereur
2290          write(imprim,*) 'pas de triangle d''abscisse minimale'
2291          ierr = 2
2292          goto 9990
2293       endif
2294 c
2295 c     chainage des 3 aretes du triangle ntmin
2296 c     =======================================
2297 c     la premiere arete du chainage des aretes traitees
2298       noar1 = abs( noartr(1,ntmin) )
2299       na0   = abs( noartr(2,ntmin) )
2300 c     elle est chainee sur la seconde arete du triangle ntmin
2301       nosoar(6,noar1) = na0
2302 c     les 2 autres aretes du triangle ntmin sont chainees
2303       na1 = abs( noartr(3,ntmin) )
2304 c     la seconde est chainee sur la troisieme arete
2305       nosoar(6,na0) = na1
2306 c     la troisieme n'a pas de suivante
2307       nosoar(6,na1) = 0
2308 c
2309 c     le triangle ntmin est a l'exterieur du domaine
2310 c     tous les triangles externes sont marques -123 456 789
2311 c     les triangles de l'autre cote d'une arete sur une ligne
2312 c     sont marques: no de la ligne de l'arete * signe oppose
2313 c     =======================================================
2314       ligne0 = 0
2315       ligne  = -123 456 789
2316 c
2317  40   if( noar1 .ne. 0 ) then
2318 c
2319 c        l'arete noar1 du tableau nosoar est a traiter
2320 c        ---------------------------------------------
2321          noar = noar1
2322 c        l'arete suivante devient la premiere a traiter ensuite
2323          noar1 = nosoar(6,noar1)
2324 c        l'arete noar est traitee
2325          nosoar(6,noar) = -3
2326 c
2327          do 60 i=4,5
2328 c
2329 c           l'un des 2 triangles de l'arete
2330             nt = nosoar(i,noar)
2331             if( nt .gt. 0 ) then
2332 c
2333 c              triangle deja traite pour une ligne anterieure?
2334                if(     letrsu(nt)  .ne. 0      .and.
2335      %             abs(letrsu(nt)) .ne. ligne ) goto 60
2336 c
2337 cccc              trace du triangle nt en couleur ligne0
2338 ccc               call mttrtr( pxyd,   nt, moartr, noartr, mosoar, nosoar,
2339 ccc     %                      ligne0, ncnoir )
2340 c
2341 c              le triangle est marque avec la valeur de ligne
2342                letrsu(nt) = ligne
2343 c
2344 c              chainage eventuel des autres aretes de ce triangle
2345 c              si ce n'est pas encore fait
2346                do 50 j=1,3
2347 c
2348 c                 le numero na de l'arete j du triangle nt dans nosoar
2349                   na = abs( noartr(j,nt) )
2350                   if( nosoar(6,na) .ne. -2 ) goto 50
2351 c
2352 c                 le numero de 1 a nblftr dans nulftr de la ligne de l'arete
2353                   nl = nosoar(3,na)
2354 c
2355 c                 si l'arete est sur une ligne fermee differente de celle envelo
2356 c                 et non marquee alors examen du triangle oppose
2357                   if( nl .gt. 0 ) then
2358 c
2359                      if( nl .eq. ligne0 ) goto 50
2360 c
2361 c                    arete frontaliere de ligne non traitee
2362 c                    => passage de l'autre cote de la ligne
2363 c                    le triangle de l'autre cote de la ligne est recherche
2364                      if( nt .eq. abs( nosoar(4,na) ) ) then
2365                         nt2 = 5
2366                      else
2367                         nt2 = 4
2368                      endif
2369                      nt2 = abs( nosoar(nt2,na) )
2370                      if( nt2 .gt. 0 ) then
2371 c
2372 c                       le triangle nt2 de l'autre cote est marque avec le
2373 c                       avec le signe oppose de celui de ligne
2374                         if( ligne .ge. 0 ) then
2375                            lsigne = -1
2376                         else
2377                            lsigne =  1
2378                         endif
2379                         letrsu(nt2) = lsigne * nl
2380 c
2381 c                       temoin de ligne a traiter ensuite dans nulftr
2382                         nulftr(nl) = -abs( nulftr(nl) )
2383 c
2384 cccc                       trace du triangle nt2 en jaune borde de magenta
2385 ccc                        call mttrtr( pxyd,nt2,
2386 ccc     %                               moartr,noartr,mosoar,nosoar,
2387 ccc     %                               ncjaun, ncmage )
2388 c
2389 c                       l'arete est traitee
2390                         nosoar(6,na) = -3
2391 c
2392                      endif
2393 c
2394 c                    l'arete est traitee
2395                      goto 50
2396 c
2397                   endif
2398 c
2399 c                 arete non traitee => elle est chainee
2400                   nosoar(6,na) = noar1
2401                   noar1 = na
2402 c
2403  50            continue
2404 c
2405             endif
2406  60      continue
2407 c
2408          goto 40
2409       endif
2410 c     les triangles de la ligne fermee ont tous ete marques
2411 c     plus d'arete chainee
2412 c
2413 c     recherche d'une nouvelle ligne fermee a traiter
2414 c     ===============================================
2415  65   do 70 nl=1,nblftr
2416          if( nulftr(nl) .lt. 0 ) goto 80
2417  70   continue
2418 c     plus de ligne fermee a traiter
2419       goto 110
2420 c
2421 c     tous les triangles de cette composante connexe
2422 c     entre ligne et ligne0 vont etre marques
2423 c     ==============================================
2424 c     remise en etat du numero de ligne
2425 c     nl est le numero de la ligne dans nulftr a traiter
2426  80   nulftr(nl) = -nulftr(nl)
2427       do 90 nt2=1,ndtri0
2428          if( abs(letrsu(nt2)) .eq. nl ) goto 92
2429  90   continue
2430 c
2431 c     recherche de l'arete j du triangle nt2 avec ce numero de ligne nl
2432  92   do 95 j=1,3
2433 c
2434 c        le numero de l'arete j du triangle dans nosoar
2435          noar1 = 0
2436          na0   = abs( noartr(j,nt2) )
2437          if( nl .eq. nosoar(3,na0) ) then
2438 c
2439 c           na0 est l'arete de ligne nl
2440 c           l'arete suivante du triangle nt2
2441             i   = mod(j,3) + 1
2442 c           le numero dans nosoar de l'arete i de nt2
2443             na1 = abs( noartr(i,nt2) )
2444             if( nosoar(6,na1) .eq. -2 ) then
2445 c              arete non traitee => elle est la premiere du chainage
2446                noar1 = na1
2447 c              pas de suivante dans ce chainage
2448                nosoar(6,na1) = 0
2449             else
2450                na1 = 0
2451             endif
2452 c
2453 c           l'eventuelle seconde arete suivante
2454             i  = mod(i,3) + 1
2455             na = abs( noartr(i,nt2) )
2456             if( nosoar(6,na) .eq. -2 ) then
2457                if( na1 .eq. 0 ) then
2458 c                 1 arete non traitee et seule a chainer
2459                   noar1 = na
2460                   nosoar(6,na) = 0
2461                else
2462 c                 2 aretes a chainer
2463                   noar1 = na
2464                   nosoar(6,na) = na1
2465                endif
2466             endif
2467 c
2468             if( noar1 .gt. 0 ) then
2469 c
2470 c              il existe au moins une arete a visiter pour ligne
2471 c              marquage des triangles internes a la ligne nl
2472                ligne  = letrsu(nt2)
2473                ligne0 = nl
2474                goto 40
2475 c
2476             else
2477 c
2478 c              nt2 est le seul triangle de la ligne fermee
2479                goto 65
2480 c
2481             endif
2482          endif
2483  95   continue
2484 c
2485 c     reperage des sommets internes ou externes dans nslign
2486 c     nslign(sommet externe au domaine)=-1
2487 c     nslign(sommet interne au domaine)= 0
2488 c     =====================================================
2489  110  do 170 ns1=1,nbsomm
2490 c        tout sommet non sur la frontiere ou interne impose
2491 c        est suppose externe
2492          if( nslign(ns1) .eq. 0 ) nslign(ns1) = -1
2493  170  continue
2494 c
2495 c     les triangles externes sont marques vides dans le tableau noartr
2496 c     ================================================================
2497       nbtria = 0
2498       do 200 nt=1,ndtri0
2499 c
2500          if( letrsu(nt) .le. 0 ) then
2501 c
2502 c           triangle nt externe
2503             if( noartr(1,nt) .ne. 0 ) then
2504 c              la premiere arete est annulee
2505                noartr(1,nt) = 0
2506 c              le triangle nt est considere comme etant vide
2507                noartr(2,nt) = n1artr
2508                n1artr = nt
2509             endif
2510 c
2511          else
2512 c
2513 c           triangle nt interne
2514             nbtria = nbtria + 1
2515             letrsu(nt) = nbtria
2516 c
2517 c           marquage des 3 sommets du triangle nt
2518             do 190 i=1,3
2519 c              le numero nosoar de l'arete i du triangle nt
2520                noar = abs( noartr(i,nt) )
2521 c              le numero des 2 sommets
2522                ns1 = nosoar(1,noar)
2523                ns2 = nosoar(2,noar)
2524 c              mise a jour du numero d'une arete des 2 sommets de l'arete
2525                noarst( ns1 ) = noar
2526                noarst( ns2 ) = noar
2527 c              ns1 et ns2 sont des sommets de la triangulation du domaine
2528                if( nslign(ns1) .lt. 0 ) nslign(ns1)=0
2529                if( nslign(ns2) .lt. 0 ) nslign(ns2)=0
2530  190        continue
2531 c
2532          endif
2533 c
2534  200  continue
2535 c     ici tout sommet externe ns verifie nslign(ns)=-1
2536 c
2537 c     les triangles externes sont mis a zero dans nosoar
2538 c     ==================================================
2539       do 300 noar=1,mxsoar
2540 c
2541          if( nosoar(1,noar) .gt. 0 ) then
2542 c
2543 c           le second triangle de l'arete noar
2544             nt = nosoar(5,noar)
2545             if( nt .gt. 0 ) then
2546 c              si le triangle nt est externe
2547 c              alors il est supprime pour l'arete noar
2548                if( letrsu(nt) .le. 0 ) nosoar(5,noar)=0
2549             endif
2550 c
2551 c           le premier triangle de l'arete noar
2552             nt = nosoar(4,noar)
2553             if( nt .gt. 0 ) then
2554                if( letrsu(nt) .le. 0 ) then
2555 c                 si le triangle nt est externe
2556 c                 alors il est supprime pour l'arete noar
2557 c                 et l'eventuel triangle oppose prend sa place
2558 c                 en position 4 de nosoar
2559                   if( nosoar(5,noar) .gt. 0 ) then
2560                      nosoar(4,noar)=nosoar(5,noar)
2561                      nosoar(5,noar)=0
2562                   else
2563                      nosoar(4,noar)=0
2564                   endif
2565                endif
2566             endif
2567          endif
2568 c
2569  300  continue
2570 c
2571 c     remise en etat pour eviter les modifications de ladefi
2572  9990 do 9991 nl=1,nblftr
2573          if( nulftr(nl) .lt. 0 ) nulftr(nl)=-nulftr(nl)
2574  9991 continue
2575       return
2576       end
2577
2578
2579
2580       subroutine trp1st( ns,     noarst, mosoar, nosoar, moartr, noartr,
2581      %                   mxpile, lhpile, lapile )
2582 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2583 c but :   recherche des triangles de noartr partageant le sommet ns
2584 c -----
2585 c         limite: un camembert de centre ns entame 2 fois
2586 c                 ne donne que l'une des parties
2587 c
2588 c entrees:
2589 c --------
2590 c ns     : numero du sommet
2591 c noarst : noarst(i) numero d'une arete de sommet i
2592 c mosoar : nombre maximal d'entiers par arete et
2593 c          indice dans nosoar de l'arete suivante dans le hachage
2594 c nosoar : numero des 2 sommets , no ligne, 2 triangles de l'arete,
2595 c          chainage des aretes frontalieres, chainage du hachage des aretes
2596 c moartr : nombre maximal d'entiers par arete du tableau noartr
2597 c noartr : les 3 aretes des triangles +-arete1, +-arete2, +-arete3
2598 c mxpile : nombre maximal de triangles empilables
2599 c
2600 c sorties :
2601 c --------
2602 c lhpile : >0 nombre de triangles empiles
2603 c          =0       si impossible de tourner autour du point
2604 c          =-lhpile si apres butee sur la frontiere il y a a nouveau
2605 c          butee sur la frontiere . a ce stade on ne peut dire si tous
2606 c          les triangles ayant ce sommet ont ete recenses
2607 c          ce cas arrive seulement si le sommet est sur la frontiere
2608 c lapile : numero dans noartr des triangles de sommet ns
2609 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2610 c auteur : alain perronnet  analyse numerique paris upmc       mars 1997
2611 c....................................................................012
2612       common / unites / lecteu, imprim, nunite(30)
2613       integer           noartr(moartr,*),
2614      %                  nosoar(mosoar,*),
2615      %                  noarst(*)
2616       integer           lapile(1:mxpile)
2617       integer           nosotr(3)
2618 c
2619 c     la premiere arete de sommet ns
2620       nar = noarst( ns )
2621       if( nar .le. 0 ) then
2622          write(imprim,*) 'trp1st: sommet',ns,' sans arete'
2623          goto 9999
2624       endif
2625 c
2626 c     l'arete nar est elle active?
2627       if( nosoar(1,nar) .le. 0 ) then
2628 ccc         write(imprim,*) 'trp1st: arete vide',nar,
2629 ccc     %                  ' st1:', nosoar(1,nar),' st2:',nosoar(2,nar)
2630          goto 9999
2631       endif
2632 c
2633 c     le premier triangle de sommet ns
2634       nt0 = abs( nosoar(4,nar) )
2635       if( nt0 .le. 0 ) then
2636          write(imprim,*) 'trp1st: sommet',ns,' dans aucun triangle'
2637          goto 9999
2638       endif
2639 c
2640 c     le triangle est il interne?
2641       if( noartr(1,nt0) .eq. 0 ) goto 9999
2642 c
2643 c     le numero des 3 sommets du triangle nt0 dans le sens direct
2644       call nusotr( nt0, mosoar, nosoar, moartr, noartr, nosotr )
2645 c
2646 c     reperage du sommet ns dans le triangle nt0
2647       do 5 nar=1,3
2648          if( nosotr(nar) .eq. ns ) goto 10
2649  5    continue
2650       nta = nt0
2651       goto 9995
2652 c
2653 c     ns retrouve : le triangle nt0 est empile
2654  10   lhpile = 1
2655       lapile(1) = nt0
2656       nta = nt0
2657 c
2658 c     recherche dans le sens des aiguilles d'une montre
2659 c     (sens indirect) du triangle nt1 de l'autre cote de l'arete
2660 c     nar du triangle et en tournant autour du sommet ns
2661 c     ==========================================================
2662       noar = abs( noartr(nar,nt0) )
2663 c     le triangle nt1 oppose du triangle nt0 par l'arete noar
2664       if( nosoar(4,noar) .eq. nt0 ) then
2665          nt1 = nosoar(5,noar)
2666       else
2667          nt1 = nosoar(4,noar)
2668       endif
2669 c
2670 c     la boucle sur les triangles nt1 de sommet ns dans le sens indirect
2671 c     ==================================================================
2672       if( nt1 .gt. 0 ) then
2673 c
2674          if( noartr(1,nt1) .eq. 0 ) goto 30
2675 c
2676 c        le triangle nt1 n'a pas ete detruit. il est actif
2677 c        le triangle oppose par l'arete noar existe
2678 c        le numero des 3 sommets du triangle nt1 dans le sens direct
2679  15      call nusotr( nt1, mosoar, nosoar, moartr, noartr, nosotr )
2680 c
2681 c        reperage du sommet ns dans nt1
2682          do 20 nar=1,3
2683             if( nosotr(nar) .eq. ns ) goto 25
2684  20      continue
2685          nta = nt1
2686          goto 9995
2687 c
2688 c        nt1 est empile
2689  25      if( lhpile .ge. mxpile ) goto 9990
2690          lhpile = lhpile + 1
2691          lapile(lhpile) = nt1
2692 c
2693 c        le triangle nt1 de l'autre cote de l'arete de sommet ns
2694 c        sauvegarde du precedent triangle dans nta
2695          nta  = nt1
2696          noar = abs( noartr(nar,nt1) )
2697          if( nosoar(4,noar) .eq. nt1 ) then
2698             nt1 = nosoar(5,noar)
2699          else
2700             nt1 = nosoar(4,noar)
2701          endif
2702          if( nt1 .le. 0   ) goto 30
2703 c        le triangle suivant est a l'exterieur
2704          if( nt1 .ne. nt0 ) goto 15
2705 c
2706 c        recherche terminee par arrivee sur nt0
2707 c        les triangles forment un "cercle" de "centre" ns
2708          return
2709 c
2710       endif
2711 c
2712 c     pas de triangle voisin a nt1
2713 c     ============================
2714 c     le parcours passe par 1 des triangles exterieurs
2715 c     le parcours est inverse par l'arete de gauche
2716 c     le triangle nta est le premier triangle empile
2717  30   lhpile = 1
2718       lapile(lhpile) = nta
2719 c
2720 c     le numero des 3 sommets du triangle nta dans le sens direct
2721       call nusotr( nta, mosoar, nosoar, moartr, noartr, nosotr )
2722       do 32 nar=1,3
2723          if( nosotr(nar) .eq. ns ) goto 33
2724  32   continue
2725       goto 9995
2726 c
2727 c     l'arete qui precede (rotation / ns dans le sens direct)
2728  33   if( nar .eq. 1 ) then
2729          nar = 3
2730       else
2731          nar = nar - 1
2732       endif
2733 c
2734 c     le triangle voisin de nta dans le sens direct
2735       noar = abs( noartr(nar,nta) )
2736       if( nosoar(4,noar) .eq. nta ) then
2737          nt1 = nosoar(5,noar)
2738       else
2739          nt1 = nosoar(4,noar)
2740       endif
2741       if( nt1 .le. 0 ) then
2742 c        un seul triangle contient ns
2743          goto 70
2744       endif
2745 c
2746 c     boucle sur les triangles de sommet ns dans le sens direct
2747 c     ==========================================================
2748  40   if( noartr(1,nt1) .eq. 0 ) goto 70
2749 c
2750 c     le triangle nt1 n'a pas ete detruit. il est actif
2751 c     le numero des 3 sommets du triangle nt1 dans le sens direct
2752       call nusotr( nt1, mosoar, nosoar, moartr, noartr, nosotr )
2753 c
2754 c     reperage du sommet ns dans nt1
2755       do 50 nar=1,3
2756          if( nosotr(nar) .eq. ns ) goto 60
2757  50   continue
2758       nta = nt1
2759       goto 9995
2760 c
2761 c     nt1 est empile
2762  60   if( lhpile .ge. mxpile ) goto 9990
2763       lhpile = lhpile + 1
2764       lapile(lhpile) = nt1
2765 c
2766 c     l'arete qui precede dans le sens direct
2767       if( nar .eq. 1 ) then
2768          nar = 3
2769       else
2770          nar = nar - 1
2771       endif
2772 c
2773 c     l'arete de sommet ns dans nosoar
2774       noar = abs( noartr(nar,nt1) )
2775 c
2776 c     le triangle voisin de nta dans le sens direct
2777       nta  = nt1
2778       if( nosoar(4,noar) .eq. nt1 ) then
2779          nt1 = nosoar(5,noar)
2780       else
2781          nt1 = nosoar(4,noar)
2782       endif
2783       nta = nt1
2784       if( nt1 .gt. 0 ) goto 40
2785 c
2786 c     butee sur le trou => fin des triangles de sommet ns
2787 c     ----------------------------------------------------
2788  70   lhpile = -lhpile
2789 c     impossible ici de trouver les autres triangles de sommet ns
2790 c     les triangles de sommet ns ne forment pas une boule de centre ns
2791       return
2792 c
2793 c     saturation de la pile des triangles
2794 c     -----------------------------------
2795  9990 write(imprim,*) 'trp1st:saturation pile des triangles autour ',
2796      %'sommet',ns
2797       goto 9999
2798 c
2799 c     erreur triangle ne contenant pas le sommet ns
2800 c     ----------------------------------------------
2801  9995 write(imprim,*) 'trp1st:triangle ',nta,' st=',
2802      %   (nosotr(nar),nar=1,3),' sans le sommet' ,ns
2803 c
2804  9999 lhpile = 0
2805       return
2806       end
2807
2808
2809
2810       subroutine nusotr( nt, mosoar, nosoar, moartr, noartr, nosotr )
2811 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2812 c but :    calcul du numero des 3 sommets du triangle nt de noartr
2813 c -----    dans le sens direct (aire>0 si non degenere)
2814 c
2815 c entrees:
2816 c --------
2817 c nt     : numero du triangle dans le tableau noartr
2818 c mosoar : nombre maximal d'entiers par arete
2819 c nosoar : numero des 2 sommets , no ligne, 2 triangles, chainages en +
2820 c          sommet 1 = 0 si arete vide => sommet 2 = arete vide suivante
2821 c moartr : nombre maximal d'entiers par arete du tableau noartr
2822 c noartr : les 3 aretes des triangles +-arete1, +-arete2, +-arete3
2823 c          arete1=0 si triangle vide => arete2=triangle vide suivant
2824 c
2825 c sorties:
2826 c --------
2827 c nosotr : numero (dans le tableau pxyd) des 3 sommets du triangle
2828 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2829 c auteur : alain perronnet  analyse numerique paris upmc       mars 1997
2830 c2345x7..............................................................012
2831       integer     nosoar(mosoar,*), noartr(moartr,*), nosotr(3)
2832 c
2833 c     les 2 sommets de l'arete 1 du triangle nt dans le sens direct
2834       na = noartr( 1, nt )
2835       if( na .gt. 0 ) then
2836          nosotr(1) = 1
2837          nosotr(2) = 2
2838       else
2839          nosotr(1) = 2
2840          nosotr(2) = 1
2841          na = -na
2842       endif
2843       nosotr(1) = nosoar( nosotr(1), na )
2844       nosotr(2) = nosoar( nosotr(2), na )
2845 c
2846 c     l'arete suivante
2847       na = abs( noartr(2,nt) )
2848 c
2849 c     le sommet nosotr(3 du triangle 123
2850       nosotr(3) = nosoar( 1, na )
2851       if( nosotr(3) .eq. nosotr(1) .or. nosotr(3) .eq. nosotr(2) ) then
2852          nosotr(3) = nosoar(2,na)
2853       endif
2854       end
2855
2856
2857       subroutine tesusp( nbarpi, pxyd,   noarst,
2858      %                   mosoar, mxsoar, n1soar, nosoar,
2859      %                   moartr, mxartr, n1artr, noartr,
2860      %                   mxarcf, n1arcf, noarcf, larmin, notrcf, liarcf,
2861      %                   nbstsu, ierr )
2862 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2863 c but :   supprimer de la triangulation les sommets de te trop proches
2864 c -----   soit d'un sommet frontalier ou point interne impose
2865 c         soit d'une arete frontaliere
2866 c
2867 c         attention: le chainage lchain de nosoar devient celui des cf
2868 c
2869 c entrees:
2870 c --------
2871 c nbarpi : numero du dernier point interne impose par l'utilisateur
2872 c pxyd   : tableau des coordonnees 2d des points
2873 c          par point : x  y  distance_souhaitee
2874 c mosoar : nombre maximal d'entiers par arete et
2875 c          indice dans nosoar de l'arete suivante dans le hachage
2876 c mxsoar : nombre maximal d'aretes stockables dans le tableau nosoar
2877 c          attention: mxsoar>3*mxsomm obligatoire!
2878 c moartr : nombre maximal d'entiers par arete du tableau noartr
2879 c mxarcf : nombre de variables des tableaux n1arcf, noarcf, larmin, notrcf
2880 c
2881 c modifies:
2882 c ---------
2883 c noarst : noarst(i) numero d'une arete de sommet i
2884 c n1soar : no de l'eventuelle premiere arete libre dans le tableau nosoar
2885 c          chainage des vides suivant en 3 et precedant en 2 de nosoar
2886 c nosoar : numero des 2 sommets , no ligne, 2 triangles de l'arete,
2887 c          chainage des aretes frontalieres, chainage du hachage des aretes
2888 c          hachage des aretes = nosoar(1)+nosoar(2)*2
2889 c          avec mxsoar>=3*mxsomm
2890 c          une arete i de nosoar est vide <=> nosoar(1,i)=0 et
2891 c          nosoar(2,arete vide)=l'arete vide qui precede
2892 c          nosoar(3,arete vide)=l'arete vide qui suit
2893 c n1artr : numero du premier triangle vide dans le tableau noartr
2894 c          le chainage des triangles vides se fait sur noartr(2,.)
2895 c noartr : les 3 aretes des triangles +-arete1, +-arete2, +-arete3
2896 c          arete1 = 0 si triangle vide => arete2 = triangle vide suivant
2897 c
2898 c
2899 c auxiliaires :
2900 c -------------
2901 c n1arcf : tableau (0:mxarcf) auxiliaire d'entiers
2902 c noarcf : tableau (3,mxarcf) auxiliaire d'entiers
2903 c larmin : tableau ( mxarcf ) auxiliaire d'entiers
2904 c notrcf : tableau ( mxarcf ) auxiliaire d'entiers
2905 c liarcf : tableau ( mxarcf ) auxiliaire d'entiers
2906 c
2907 c sortie :
2908 c --------
2909 c nbstsu : nombre de sommets de te supprimes
2910 c ierr   : =0 si pas d'erreur
2911 c          >0 si une erreur est survenue
2912 c          11 algorithme defaillant
2913 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2914 c auteur : alain perronnet  analyse numerique paris upmc       mars 1997
2915 c....................................................................012
2916 c      parameter       ( quamal=0.3 ) => ok
2917 c      parameter       ( quamal=0.4 ) => pb pour le test ocean
2918 c      parameter       ( quamal=0.5 ) => pb pour le test ocean
2919 c
2920       parameter       ( quamal=0.333, lchain=6 )
2921       common / unites / lecteu,imprim,intera,nunite(29)
2922       double precision  pxyd(3,*), qualit
2923       integer           nosoar(mosoar,mxsoar),
2924      %                  noartr(moartr,*),
2925      %                  noarst(*),
2926      %                  n1arcf(0:mxarcf),
2927      %                  noarcf(3,mxarcf),
2928      %                  larmin(mxarcf),
2929      %                  notrcf(mxarcf),
2930      %                  liarcf(mxarcf)
2931 c
2932       integer           nosotr(3)
2933       equivalence      (nosotr(1),ns1), (nosotr(2),ns2),
2934      %                 (nosotr(3),ns3)
2935 c
2936 c     le nombre de sommets de te supprimes
2937       nbstsu = 0
2938 c
2939 c     initialisation du chainage des aretes des cf => 0 arete de cf
2940       do 10 narete=1,mxsoar
2941          nosoar( lchain, narete ) = -1
2942  10   continue
2943 c
2944 c     boucle sur l'ensemble des sommets frontaliers ou points internes
2945 c     ================================================================
2946       do 100 ns = 1, nbarpi
2947 c
2948 cccc        le nombre de sommets supprimes pour ce sommet ns
2949 ccc         nbsuns = 0
2950 c
2951 c        la qualite minimale au dessous de laquelle le point proche
2952 c        interne est supprime
2953          quaopt = quamal
2954 c
2955 c        une arete de sommet ns
2956  15      narete = noarst( ns )
2957          if( narete .le. 0 ) then
2958 c           erreur: le point appartient a aucune arete
2959             write(imprim,*) 'sommet ',ns,' dans aucune arete'
2960             pause
2961             ierr = 11
2962             return
2963          endif
2964 c
2965 c        recherche des triangles de sommet ns
2966 c        ils doivent former un contour ferme de type etoile
2967          call trp1st( ns, noarst, mosoar, nosoar, moartr, noartr,
2968      %                mxarcf, nbtrcf, notrcf )
2969          if( nbtrcf .le. 0 ) then
2970 c           erreur: impossible de trouver tous les triangles de sommet ns
2971 c           seule une partie est a priori retrouvee
2972             nbtrcf = -nbtrcf
2973          endif
2974 c
2975 c        boucle sur les triangles de l'etoile du sommet ns
2976          quamin = 2.0
2977          do 20 i=1,nbtrcf
2978 c
2979 c           le numero des 3 sommets du triangle nt
2980             nt = notrcf(i)
2981             call nusotr( nt, mosoar, nosoar, moartr, noartr,
2982      %                   nosotr )
2983 c           nosotr(1:3) est en equivalence avec ns1, ns2, ns3
2984 c
2985 c           la qualite du triangle ns1 ns2 ns3
2986             call qutr2d( pxyd(1,ns1), pxyd(1,ns2), pxyd(1,ns3), qualit )
2987             if( qualit .lt. quamin ) then
2988                quamin = qualit
2989                ntqmin = nt
2990             endif
2991  20      continue
2992 c
2993 c        bilan sur la qualite des triangles de sommet ns
2994          if( quamin .lt. quaopt ) then
2995 c
2996 c           recherche du sommet de ntqmin le plus proche et non frontalier
2997 c           ==============================================================
2998 c           le numero des 3 sommets du triangle nt
2999             call nusotr( ntqmin, mosoar, nosoar, moartr, noartr,
3000      %                   nosotr )
3001             nste   = 0
3002             quamin = 1e28
3003             do 30 j=1,3
3004                if( nosotr(j) .ne. ns .and. nosotr(j) .gt. nbarpi ) then
3005                   d = (pxyd(1,nosotr(j))-pxyd(1,ns))**2
3006      %              + (pxyd(2,nosotr(j))-pxyd(2,ns))**2
3007                   if( d .lt. quamin ) then
3008                      quamin = d
3009                      nste   = j
3010                   endif
3011                endif
3012  30         continue
3013 c
3014             if( nste .gt. 0 ) then
3015 c
3016 c              nste est le sommet le plus proche de ns de ce
3017 c              triangle de mauvaise qualite et sommet non encore traite
3018                nste = nosotr( nste )
3019 c
3020 c              nste est un sommet de triangle equilateral
3021 c              => le sommet nste va etre supprime
3022 c              ==========================================
3023                call te1stm( nste,   pxyd,   noarst,
3024      %                      mosoar, mxsoar, n1soar, nosoar,
3025      %                      moartr, mxartr, n1artr, noartr,
3026      %                      mxarcf, n1arcf, noarcf,
3027      %                      larmin, notrcf, liarcf, ierr )
3028                if( ierr .eq. 0 ) then
3029 c                 un sommet de te supprime de plus
3030                   nbstsu = nbstsu + 1
3031                else if( ierr .lt. 0 ) then
3032 c                 le sommet nste est externe donc non supprime
3033 c                 ou bien le sommet nste est le centre d'un cf dont toutes
3034 c                 les aretes simples sont frontalieres
3035 c                 dans les 2 cas le sommet n'est pas supprime
3036                   ierr = 0
3037                   goto 100
3038                else
3039 c                 erreur motivant un arret de la triangulation
3040                   return
3041                endif
3042 c
3043 c              boucle jusqu'a obtenir une qualite suffisante
3044 c              si triangulation tres irreguliere =>
3045 c              destruction de beaucoup de points internes
3046 c              les 2 variables suivantes brident ces destructions massives
3047 ccc               nbsuns = nbsuns + 1
3048                quaopt = quaopt * 0.8
3049 ccc               if( nbsuns .le. 5 ) goto 15
3050                goto 15
3051             endif
3052          endif
3053 c
3054  100  continue
3055       end
3056
3057
3058       subroutine teamqa( nutysu,
3059      %                   noarst, mosoar, mxsoar, n1soar, nosoar,
3060      %                   moartr, mxartr, n1artr, noartr,
3061      %                   mxtrcf, notrcf, nostbo,
3062      %                   n1arcf, noarcf, larmin,
3063      %                   comxmi, nbarpi, nbsomm, mxsomm, pxyd, nslign,
3064      %                   ierr )
3065 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3066 c but:    si la taille de l'arete moyenne est >ampli*taille souhaitee
3067 c ----    alors ajout d'un sommet barycentre du plus grand triangle
3068 c               de sommet ns
3069 c         si la taille de l'arete moyenne est <ampli/2*taille souhaitee
3070 c         alors suppression du sommet ns
3071 c         sinon le sommet ns devient le barycentre pondere de ses voisins
3072 c
3073 c         remarque: ampli est defini dans $mefisto/mail/tehote.f
3074 c         et doit avoir la meme valeur pour eviter trop de modifications
3075 c
3076 c entrees:
3077 c --------
3078 c nutysu : numero de traitement de areteideale() selon le type de surface
3079 c          0 pas d'emploi de la fonction areteideale() => aretmx active
3080 c          1 il existe une fonction areteideale()
3081 c            dont seules les 2 premieres composantes de uv sont actives
3082 c          autres options a definir...
3083 c noarst : noarst(i) numero d'une arete de sommet i
3084 c mosoar : nombre maximal d'entiers par arete et
3085 c          indice dans nosoar de l'arete suivante dans le hachage
3086 c mxsoar : nombre maximal d'aretes frontalieres declarables
3087 c n1soar : numero de la premiere arete vide dans le tableau nosoar
3088 c nosoar : numero des 2 sommets , no ligne, 2 triangles de l'arete,
3089 c          chainage des aretes frontalieres, chainage du hachage des aretes
3090 c moartr : nombre maximal d'entiers par arete du tableau noartr
3091 c mxartr : nombre maximal de triangles declarables dans noartr
3092 c n1artr : numero du premier triangle vide dans le tableau noartr
3093 c          le chainage des triangles vides se fait sur noartr(2,.)
3094 c noartr : les 3 aretes des triangles +-arete1, +-arete2, +-arete3
3095 c mxtrcf : nombre maximal de triangles empilables
3096 c nbarpi : numero du dernier sommet frontalier ou interne impose
3097 c nslign : tableau du numero de sommet dans sa ligne pour chaque
3098 c          sommet frontalier
3099 c          numero du point dans le lexique point si interne impose
3100 c          0 si le point est interne non impose par l'utilisateur
3101 c         -1 si le sommet est externe au domaine
3102 c comxmi : min et max des coordonneees des sommets du maillage
3103 c
3104 c modifies :
3105 c ----------
3106 c nbsomm : nombre actuel de sommets de la triangulation
3107 c          (certains sommets internes ont ete desactives ou ajoutes)
3108 c pxyd   : tableau des coordonnees 2d des points
3109 c
3110 c auxiliaires:
3111 c ------------
3112 c notrcf : tableau ( mxtrcf ) auxiliaire d'entiers
3113 c          numero dans noartr des triangles de sommet ns
3114 c nostbo : tableau ( mxtrcf ) auxiliaire d'entiers
3115 c          numero dans pxyd des sommets des aretes simples de la boule
3116 c n1arcf : tableau (0:mxtrcf) auxiliaire d'entiers
3117 c noarcf : tableau (3,mxtrcf) auxiliaire d'entiers
3118 c larmin : tableau ( mxtrcf ) auxiliaire d'entiers
3119 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3120 c auteur : alain perronnet  analyse numerique paris upmc       juin 1997
3121 c....................................................................012
3122       double precision  ampli,ampli2
3123       parameter        (ampli=1.34d0,ampli2=ampli/2d0)
3124       parameter        (lchain=6)
3125       common / unites / lecteu, imprim, nunite(30)
3126       double precision  pxyd(3,*)
3127       double precision  ponder, ponde1, xbar, ybar, x, y, surtd2
3128       double precision  d, dmoy
3129       double precision  d2d3(3,3)
3130       real              origin(3), xyz(3)
3131       integer           noartr(moartr,*),
3132      %                  nosoar(mosoar,*),
3133      %                  noarst(*),
3134      %                  notrcf(mxtrcf),
3135      %                  nslign(*),
3136      %                  nostbo(*),
3137      %                  n1arcf(0:mxtrcf),
3138      %                  noarcf(3,mxtrcf),
3139      %                  larmin(mxtrcf)
3140       double precision  comxmi(3,2)
3141       integer           nosotr(3)
3142 c
3143 c     le nombre d'iterations pour ameliorer la qualite
3144       nbitaq = 4
3145       ier    = 0
3146 c
3147 c     initialisation du parcours
3148       nbs1 = nbsomm
3149       nbs2 = nbarpi + 1
3150       nbs3 = -1
3151 c
3152       do 5000 iter=1,nbitaq
3153 c
3154 c        le nombre de sommets supprimes
3155          nbstsu = 0
3156          nbbaaj = 0
3157 c
3158 c        coefficient de ponderation croissant avec les iterations
3159          ponder = min( 1d0, ( 50 + (50*iter)/nbitaq ) * 0.01d0 )
3160          ponde1 = 1d0 - ponder
3161 c
3162 c        l'ordre du parcours dans le sens croissant ou decroissant
3163          nt   = nbs1
3164          nbs1 = nbs2
3165          nbs2 = nt
3166 c        alternance du parcours
3167          nbs3 = -nbs3
3168 c
3169          do 1000 ns = nbs1, nbs2, nbs3
3170 c
3171 c           le sommet est il interne au domaine?
3172             if( nslign(ns) .ne. 0 ) goto 1000
3173 c
3174 c           existe-t-il une arete de sommet ns ?
3175  10         noar = noarst( ns )
3176             if( noar .le. 0 ) goto 1000
3177 c
3178 c           le 1-er triangle de l'arete noar
3179             nt = nosoar( 4, noar )
3180             if( nt .le. 0 ) goto 1000
3181 c
3182 c           recherche des triangles de sommet ns
3183 c           ils doivent former un contour ferme de type etoile
3184             call trp1st( ns, noarst, mosoar, nosoar, moartr, noartr,
3185      %                   mxtrcf, nbtrcf, notrcf )
3186             if( nbtrcf .le. 0 ) goto 1000
3187 c
3188 c           mise a jour de la distance souhaitee
3189             if( nutysu .gt. 0 ) then
3190 c              la fonction taille_ideale(x,y,z) existe
3191 c              calcul de pxyzd(3,ns) dans le repere initial => xyz(1:3)
3192                call tetaid( nutysu, pxyd(1,ns), pxyd(2,ns),
3193      %                      pxyd(3,ns), ier )
3194             endif
3195 c
3196 c           boucle sur les triangles qui forment une boule autour du sommet ns
3197             nbstbo = 0
3198 c           chainage des aretes simples de la boule a rendre delaunay
3199             noar0  = 0
3200             do 40 i=1,nbtrcf
3201 c
3202 c              le numero de l'arete du triangle nt ne contenant pas le sommet ns
3203                nt = notrcf(i)
3204                do 20 na=1,3
3205 c                 le numero de l'arete na dans le tableau nosoar
3206                   noar = abs( noartr(na,nt) )
3207                   if( nosoar(1,noar) .ne. ns   .and.
3208      %                nosoar(2,noar) .ne. ns ) goto 25
3209  20            continue
3210 c
3211 c              construction de la liste des sommets des aretes simples
3212 c              de la boule des triangles de sommet ns
3213 c              -------------------------------------------------------
3214  25            do 35 na=1,2
3215                   ns1 = nosoar(na,noar)
3216                   do 30 j=nbstbo,1,-1
3217                      if( ns1 .eq. nostbo(j) ) goto 35
3218  30               continue
3219 c                 ns1 est un nouveau sommet a ajouter
3220                   nbstbo = nbstbo + 1
3221                   nostbo(nbstbo) = ns1
3222  35            continue
3223 c
3224 c              noar est une arete potentielle a rendre delaunay
3225                if( nosoar(3,noar) .eq. 0 ) then
3226 c                 arete non frontaliere
3227                   nosoar(lchain,noar) = noar0
3228                   noar0 = noar
3229                endif
3230 c
3231  40         continue
3232 c
3233 c           calcul des 2 coordonnees du barycentre de la boule du sommet ns
3234 c           calcul de la longueur moyenne des aretes issues du sommet ns
3235 c           ---------------------------------------------------------------
3236             xbar = 0d0
3237             ybar = 0d0
3238             dmoy = 0d0
3239             do 50 i=1,nbstbo
3240                x    = pxyd(1,nostbo(i))
3241                y    = pxyd(2,nostbo(i))
3242                xbar = xbar + x
3243                ybar = ybar + y
3244                dmoy = dmoy + sqrt( (x-pxyd(1,ns))**2+(y-pxyd(2,ns))**2 )
3245  50         continue
3246             dmoy = dmoy / nbstbo
3247 c
3248 c           pas de modification de la topologie lors de la derniere iteration
3249 c           =================================================================
3250             if( iter .eq. nbitaq ) goto 200
3251 c
3252 c           si la taille de l'arete moyenne est >ampli*taille souhaitee
3253 c           alors ajout d'un sommet barycentre du plus grand triangle
3254 c                 de sommet ns
3255 c           ===========================================================
3256             if( dmoy .gt. ampli*pxyd(3,ns) ) then
3257 c
3258                dmoy = 0d0
3259                do 150 i=1,nbtrcf
3260 c                 recherche du plus grand triangle en surface
3261                   call nusotr( notrcf(i), mosoar, nosoar,
3262      %                         moartr, noartr, nosotr )
3263                   d  = surtd2( pxyd(1,nosotr(1)),
3264      %                         pxyd(1,nosotr(2)),
3265      %                         pxyd(1,nosotr(3)) )
3266                   if( d .gt. dmoy ) then
3267                      dmoy = d
3268                      imax = i
3269                   endif
3270  150           continue
3271 c
3272 c              ajout du barycentre du triangle notrcf(imax)
3273                nt = notrcf( imax )
3274                call nusotr( nt, mosoar, nosoar,
3275      %                      moartr, noartr, nosotr )
3276                if( nbsomm .ge. mxsomm ) then
3277                   write(imprim,*) 'saturation du tableau pxyd'
3278 c                 abandon de l'amelioration du sommet ns
3279                   goto 9999
3280                endif
3281                nbsomm = nbsomm + 1
3282                do 160 i=1,3
3283                   pxyd(i,nbsomm) = ( pxyd(i,nosotr(1))
3284      %                             + pxyd(i,nosotr(2))
3285      %                             + pxyd(i,nosotr(3)) ) / 3d0
3286  160           continue
3287 c
3288                if( nutysu .gt. 0 ) then
3289 c                 la fonction taille_ideale(x,y,z) existe
3290 c                 calcul de pxyzd(3,nbsomm) dans le repere initial => xyz(1:3)
3291                   call tetaid( nutysu, pxyd(1,nbsomm), pxyd(2,nbsomm),
3292      %                         pxyd(3,nbsomm), ier )
3293                endif
3294 c
3295 c              sommet interne a la triangulation
3296                nslign(nbsomm) = 0
3297 c
3298 c              les 3 aretes du triangle nt sont a rendre delaunay
3299                do 170 i=1,3
3300                   noar = abs( noartr(i,nt) )
3301                   if( nosoar(3,noar) .eq. 0 ) then
3302 c                    arete non frontaliere
3303                      if( nosoar(lchain,noar) .lt. 0 ) then
3304 c                       arete non encore chainee
3305                         nosoar(lchain,noar) = noar0
3306                         noar0 = noar
3307                      endif
3308                   endif
3309  170           continue
3310 c
3311 c              triangulation du triangle de barycentre nbsomm
3312 c              protection a ne pas modifier sinon erreur!
3313                call tr3str( nbsomm, nt,
3314      %                      mosoar, mxsoar, n1soar, nosoar,
3315      %                      moartr, mxartr, n1artr, noartr,
3316      %                      noarst,
3317      %                      nosotr, ierr )
3318                if( ierr .ne. 0 ) goto 9999
3319 c
3320 c              un barycentre ajoute de plus
3321                nbbaaj = nbbaaj + 1
3322 c
3323 c              les aretes chainees de la boule sont rendues delaunay
3324                goto 900
3325 c
3326             endif
3327 c
3328 c           si la taille de l'arete moyenne est <ampli/2*taille souhaitee
3329 c           alors suppression du sommet ns
3330 c           =============================================================
3331             if( dmoy .lt. ampli2*pxyd(3,ns) ) then
3332 c              remise a -1 du chainage des aretes peripheriques de la boule ns
3333                noar = noar0
3334  90            if( noar .gt. 0 ) then
3335 c                 protection du no de l'arete suivante
3336                   na = nosoar(lchain,noar)
3337 c                 l'arete interne est remise a -1
3338                   nosoar(lchain,noar) = -1
3339 c                 l'arete suivante
3340                   noar = na
3341                   goto 90
3342                endif
3343                call te1stm( ns,     pxyd,   noarst,
3344      %                      mosoar, mxsoar, n1soar, nosoar,
3345      %                      moartr, mxartr, n1artr, noartr,
3346      %                      mxtrcf, n1arcf, noarcf,
3347      %                      larmin, notrcf, nostbo,
3348      %                      ierr )
3349                if( ierr .lt. 0 ) then
3350 c                 le sommet ns est externe donc non supprime
3351 c                 ou bien le sommet ns est le centre d'un cf dont toutes
3352 c                 les aretes simples sont frontalieres
3353 c                 dans les 2 cas le sommet ns n'est pas supprime
3354                   ierr = 0
3355                   goto 200
3356                else if( ierr .gt. 0 ) then
3357 c                 erreur irrecuperable
3358                   goto 9999
3359                endif
3360                nbstsu = nbstsu + 1
3361                goto 1000
3362 c
3363             endif
3364 c
3365 c           les 2 coordonnees du barycentre des sommets des aretes
3366 c           simples de la boule du sommet ns
3367 c           ======================================================
3368  200        xbar = xbar / nbstbo
3369             ybar = ybar / nbstbo
3370 c
3371 c           ponderation pour eviter les degenerescenses
3372             pxyd(1,ns) = ponde1 * pxyd(1,ns) + ponder * xbar
3373             pxyd(2,ns) = ponde1 * pxyd(2,ns) + ponder * ybar
3374 c
3375             if( nutysu .gt. 0 ) then
3376 c              la fonction taille_ideale(x,y,z) existe
3377 c              calcul de pxyzd(3,ns) dans le repere initial => xyz(1:3)
3378                call tetaid( nutysu, pxyd(1,ns), pxyd(2,ns),
3379      %                      pxyd(3,ns), ier )
3380             endif
3381 c
3382 c           les aretes chainees de la boule sont rendues delaunay
3383  900        call tedela( pxyd,   noarst,
3384      %                   mosoar, mxsoar, n1soar, nosoar, noar0,
3385      %                   moartr, mxartr, n1artr, noartr, modifs )
3386 c
3387  1000    continue
3388 c
3389 ccc         write(imprim,11000) nbstsu, nbbaaj
3390 ccc11000 format( i6,' sommets supprimes ' ,
3391 ccc     %        i6,' barycentres ajoutes' )
3392 c
3393 c        mise a jour pour ne pas oublier les nouveaux sommets
3394          if( nbs1 .gt. nbs2 ) then
3395             nbs1 = nbsomm
3396          else
3397             nbs2 = nbsomm
3398          endif
3399 c
3400  5000 continue
3401 c
3402  9999 return
3403       end
3404
3405
3406       subroutine teamsf( nutysu,
3407      %                   noarst, mosoar, mxsoar, n1soar, nosoar,
3408      %                   moartr, mxartr, n1artr, noartr,
3409      %                   mxtrcf, notrcf, nostbo,
3410      %                   n1arcf, noarcf, larmin,
3411      %                   comxmi, nbarpi, nbsomm, mxsomm, pxyd, nslign,
3412      %                   ierr )
3413 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3414 c but :    modification de la topologie des triangles autour des
3415 c -----    sommets frontaliers et mise en triangulation delaunay locale
3416 c
3417 c entrees:
3418 c --------
3419 c nutysu : numero de traitement de areteideale() selon le type de surface
3420 c          0 pas d'emploi de la fonction areteideale() => aretmx active
3421 c          1 il existe une fonction areteideale()
3422 c            dont seules les 2 premieres composantes de uv sont actives
3423 c          autres options a definir...
3424 c noarst : noarst(i) numero d'une arete de sommet i
3425 c mosoar : nombre maximal d'entiers par arete et
3426 c          indice dans nosoar de l'arete suivante dans le hachage
3427 c mxsoar : nombre maximal d'aretes frontalieres declarables
3428 c n1soar : numero de la premiere arete vide dans le tableau nosoar
3429 c nosoar : numero des 2 sommets , no ligne, 2 triangles de l'arete,
3430 c          chainage des aretes frontalieres, chainage du hachage des aretes
3431 c moartr : nombre maximal d'entiers par arete du tableau noartr
3432 c mxartr : nombre maximal de triangles declarables dans noartr
3433 c n1artr : numero du premier triangle vide dans le tableau noartr
3434 c          le chainage des triangles vides se fait sur noartr(2,.)
3435 c noartr : les 3 aretes des triangles +-arete1, +-arete2, +-arete3
3436 c mxtrcf : nombre maximal de triangles empilables
3437 c nbarpi : numero du dernier sommet frontalier ou interne impose
3438 c nslign : >0 => ns numero du point dans le lexique point si interne impose
3439 c          ou => 1 000 000 * n + ns1
3440 c              ou n   est le numero (1 a nblftr) de la ligne de ce point
3441 c                 ns1 est le numero du point dans sa ligne
3442 c          = 0 si le point est interne non impose par l'utilisateur
3443 c          =-1 si le sommet est externe au domaine
3444 c comxmi : min et max des coordonneees des sommets du maillage
3445 c
3446 c modifies :
3447 c ----------
3448 c nbsomm : nombre actuel de sommets de la triangulation
3449 c          (certains sommets internes ont ete desactives ou ajoutes)
3450 c pxyd   : tableau des coordonnees 2d des points
3451 c
3452 c auxiliaires:
3453 c ------------