1 c MEFISTO : library to compute 2D triangulation from segmented boundaries
3 c Copyright (C) 2003 Laboratoire J.-L. Lions UPMC Paris
5 c This library is free software; you can redistribute it and/or
6 c modify it under the terms of the GNU Lesser General Public
7 c License as published by the Free Software Foundation; either
8 c version 2.1 of the License.
10 c This library is distributed in the hope that it will be useful,
11 c but WITHOUT ANY WARRANTY; without even the implied warranty of
12 c MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 c Lesser General Public License for more details.
15 c You should have received a copy of the GNU Lesser General Public
16 c License along with this library; if not, write to the Free Software
17 c Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 c See http://www.ann.jussieu.fr/~perronne or email Perronnet@ann.jussieu.fr
24 c Author: Alain PERRONNET
26 subroutine qutr2d( p1, p2, p3, qualite )
27 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
28 c but : calculer la qualite d'un triangle de r**2
29 c ----- 2 coordonnees des 3 sommets en double precision
33 c p1,p2,p3 : les 3 coordonnees des 3 sommets du triangle
34 c sens direct pour une surface et qualite >0
37 c qualite: valeur de la qualite du triangle entre 0 et 1 (equilateral)
38 c 1 etant la qualite optimale
39 c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
40 c auteur : alain perronnet analyse numerique upmc paris janvier 1995
41 c2345x7..............................................................012
42 parameter ( d2uxr3 = 3.4641016151377544d0 )
43 c d2uxr3 = 2 * sqrt(3)
44 double precision p1(2), p2(2), p3(2), qualite, a, b, c, p
46 c la longueur des 3 cotes
47 a = sqrt( (p2(1)-p1(1))**2 + (p2(2)-p1(2))**2 )
48 b = sqrt( (p3(1)-p2(1))**2 + (p3(2)-p2(2))**2 )
49 c = sqrt( (p1(1)-p3(1))**2 + (p1(2)-p3(2))**2 )
54 if ( (a*b*c) .ne. 0d0 ) then
55 c critere : 2 racine(3) * rayon_inscrit / plus longue arete
56 qualite = d2uxr3 * sqrt( abs( (p-a) / p * (p-b) * (p-c) ) )
63 c autres criteres possibles:
64 c critere : 2 * rayon_inscrit / rayon_circonscrit
65 c qualite = 8d0 * (p-a) * (p-b) * (p-c) / (a * b * c)
67 c critere : 3*sqrt(3.) * ray_inscrit / demi perimetre
68 c qualite = 3*sqrt(3.) * sqrt ((p-a)*(p-b)*(p-c) / p**3)
70 c critere : 2*sqrt(3.) * ray_inscrit / max( des aretes )
71 c qualite = 2*sqrt(3.) * sqrt( (p-a)*(p-b)*(p-c) / p ) / max(a,b,c)
75 double precision function surtd2( p1 , p2 , p3 )
76 c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
77 c but : calcul de la surface d'un triangle defini par 3 points de R**2
79 c parametres d entree :
80 c ---------------------
81 c p1 p2 p3 : les 3 fois 2 coordonnees des sommets du triangle
83 c parametre resultat :
84 c --------------------
85 c surtd2 : surface du triangle
86 c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
87 c auteur : alain perronnet analyse numerique upmc paris fevrier 1992
88 c2345x7..............................................................012
89 double precision p1(2), p2(2), p3(2)
91 c la surface du triangle
92 surtd2 = ( ( p2(1)-p1(1) ) * ( p3(2)-p1(2) )
93 % - ( p2(2)-p1(2) ) * ( p3(1)-p1(1) ) ) * 0.5d0
96 integer function nopre3( i )
97 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
98 c but : numero precedent i dans le sens circulaire 1 2 3 1 ...
100 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
101 c auteur : alain perronnet analyse numerique paris upmc fevrier 1992
102 c2345x7..............................................................012
110 integer function nosui3( i )
111 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
112 c but : numero suivant i dans le sens circulaire 1 2 3 1 ...
114 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
115 c auteur : alain perronnet analyse numerique paris upmc fevrier 1992
116 c2345x7..............................................................012
124 subroutine provec( v1 , v2 , v3 )
125 c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
126 c but : v3 vecteur = produit vectoriel de 2 vecteurs de r ** 3
130 c v1, v2 : les 2 vecteurs de 3 composantes
134 c v3 : vecteur = v1 produit vectoriel v2
135 cc++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
136 c auteur : perronnet alain upmc analyse numerique paris mars 1987
137 c2345x7..............................................................012
138 double precision v1(3), v2(3), v3(3)
140 v3( 1 ) = v1( 2 ) * v2( 3 ) - v1( 3 ) * v2( 2 )
141 v3( 2 ) = v1( 3 ) * v2( 1 ) - v1( 1 ) * v2( 3 )
142 v3( 3 ) = v1( 1 ) * v2( 2 ) - v1( 2 ) * v2( 1 )
147 subroutine norme1( n, v, ierr )
148 c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
149 c but : normalisation euclidienne a 1 d un vecteur v de n composantes
153 c n : nombre de composantes du vecteur
157 c v : le vecteur a normaliser a 1
161 c ierr : 1 si la norme de v est egale a 0
163 c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
164 c auteur : alain perronnet analyse numerique paris mars 1987
165 c ......................................................................
166 double precision v( n ), s, sqrt
170 s = s + v( i ) * v( i )
173 c test de nullite de la norme du vecteur
174 c --------------------------------------
175 if( s .le. 0.0d0 ) then
176 c norme nulle du vecteur non normalisable a 1
181 s = 1.0d0 / sqrt( s )
190 subroutine insoar( mxsomm, mosoar, mxsoar, n1soar, nosoar )
191 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
192 c but : initialiser le tableau nosoar pour le hachage des aretes
197 c mxsomm : plus grand numero de sommet d'une arete au cours du calcul
198 c mosoar : nombre maximal d'entiers par arete du tableau nosoar
199 c mxsoar : nombre maximal d'aretes stockables dans le tableau nosoar
200 c avec mxsoar>=3*mxsomm
204 c n1soar : numero de la premiere arete vide dans le tableau nosoar
205 c une arete i de nosoar est vide <=> nosoar(1,i)=0
206 c chainage des aretes vides amont et aval
207 c l'arete vide qui precede=nosoar(4,i)
208 c l'arete vide qui suit =nosoar(5,i)
209 c nosoar : numero des 2 sommets, no ligne, 2 triangles de l'arete,
210 c chainage momentan'e d'aretes, chainage du hachage des aretes
211 c hachage des aretes = min( nosoar(1), nosoar(2) )
212 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
213 c auteur : alain perronnet analyse numerique paris upmc mars 1997
214 c2345x7..............................................................012
215 integer nosoar(mosoar,mxsoar)
217 c initialisation des aretes 1 a mxsomm
220 c sommet 1 = 0 <=> temoin d'arete vide pour le hachage
223 c arete sur aucune ligne
226 c la position de l'arete interne ou frontaliere est inconnue
229 c fin de chainage du hachage pas d'arete suivante
230 nosoar( mosoar, i ) = 0
234 c la premiere arete vide chainee est la mxsomm+1 du tableau
235 c car ces aretes ne sont pas atteignables par le hachage direct
238 c initialisation des aretes vides et des chainages
239 do 20 i = n1soar, mxsoar
241 c sommet 1 = 0 <=> temoin d'arete vide pour le hachage
244 c arete sur aucune ligne
247 c chainage sur l'arete vide qui precede
248 c (si arete occupee cela deviendra le no du triangle 1 de l'arete)
251 c chainage sur l'arete vide qui suit
252 c (si arete occupee cela deviendra le no du triangle 2 de l'arete)
255 c chainages des aretes frontalieres ou internes ou ...
258 c fin de chainage du hachage
259 nosoar( mosoar, i ) = 0
263 c la premiere arete vide n'a pas de precedent
264 nosoar( 4, n1soar ) = 0
266 c la derniere arete vide est mxsoar sans arete vide suivante
267 nosoar( 5, mxsoar ) = 0
271 subroutine azeroi ( l , ntab )
272 c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
273 c but : initialisation a zero d un tableau ntab de l variables entieres
275 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
276 c auteur : alain perronnet analyse numerique upmc paris septembre 1988
277 c23456---------------------------------------------------------------012
285 subroutine fasoar( ns1, ns2, nt1, nt2, nolign,
286 % mosoar, mxsoar, n1soar, nosoar, noarst,
288 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
289 c but : former l'arete de sommet ns1-ns2 dans le hachage du tableau
290 c ----- nosoar des aretes de la triangulation
294 c ns1 ns2: numero pxyd des 2 sommets de l'arete
295 c nt1 : numero du triangle auquel appartient l'arete
296 c nt1=-1 si numero inconnu
297 c nt2 : numero de l'eventuel second triangle de l'arete si connu
298 c nt2=-1 si numero inconnu
299 c nolign : numero de la ligne de l'arete dans ladefi(wulftr-1+nolign)
300 c =0 si l'arete n'est une arete de ligne
301 c ce numero est ajoute seulement si l'arete est creee
302 c mosoar : nombre maximal d'entiers par arete du tableau nosoar
303 c mxsoar : nombre maximal d'aretes stockables dans le tableau nosoar
307 c n1soar : numero de la premiere arete vide dans le tableau nosoar
308 c une arete i de nosoar est vide <=> nosoar(1,i)=0
309 c chainage des aretes vides amont et aval
310 c l'arete vide qui precede=nosoar(4,i)
311 c l'arete vide qui suit =nosoar(5,i)
312 c nosoar : numero des 2 sommets, no ligne, 2 triangles de l'arete,
313 c chainage momentan'e d'aretes, chainage du hachage des aretes
314 c hachage des aretes = min( nosoar(1), nosoar(2) )
315 c noarst : noarst(np) numero d'une arete du sommet np
317 c ierr : si < 0 en entree pas d'affichage en cas d'erreur du type
318 c "arete appartenant a plus de 2 triangles et a creer!"
319 c si >=0 en entree affichage de ce type d'erreur
323 c noar : >0 numero de l'arete retrouvee ou ajoutee
324 c ierr : =0 si pas d'erreur
325 c =1 si le tableau nosoar est sature
326 c =2 si arete a creer et appartenant a 2 triangles distincts
327 c des triangles nt1 et nt2
328 c =3 si arete appartenant a 2 triangles distincts
329 c differents des triangles nt1 et nt2
330 c =4 si arete appartenant a 2 triangles distincts
331 c dont le second n'est pas le triangle nt2
332 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
333 c auteur : alain perronnet analyse numerique paris upmc mars 1997
334 c2345x7..............................................................012
335 common / unites / lecteu, imprim, nunite(30)
336 integer nosoar(mosoar,mxsoar), noarst(*)
339 c ajout eventuel de l'arete s1 s2 dans nosoar
343 c hachage de l'arete de sommets nu2sar
344 call hasoar( mosoar, mxsoar, n1soar, nosoar, nu2sar, noar )
345 c en sortie: noar>0 => no arete retrouvee
346 c <0 => no arete ajoutee
347 c =0 => saturation du tableau nosoar
349 if( noar .eq. 0 ) then
351 c saturation du tableau nosoar
352 write(imprim,*) 'fasoar: tableau nosoar sature'
356 else if( noar .lt. 0 ) then
358 c l'arete a ete ajoutee. initialisation des autres informations
360 c le numero de la ligne de l'arete
361 nosoar(3,noar) = nolign
362 c le triangle 1 de l'arete => le triangle nt1
364 c le triangle 2 de l'arete => le triangle nt2
367 c le sommet appartient a l'arete noar
368 noarst( nu2sar(1) ) = noar
369 noarst( nu2sar(2) ) = noar
373 c l'arete a ete retrouvee.
374 c si elle appartient a 2 triangles differents de nt1 et nt2
375 c alors il y a une erreur
376 if( nosoar(4,noar) .gt. 0 .and.
377 % nosoar(5,noar) .gt. 0 ) then
378 if( nosoar(4,noar) .ne. nt1 .and.
379 % nosoar(4,noar) .ne. nt2 .or.
380 % nosoar(5,noar) .ne. nt1 .and.
381 % nosoar(5,noar) .ne. nt2 ) then
382 c arete appartenant a plus de 2 triangles => erreur
383 if( ierr .ge. 0 ) then
384 write(imprim,*) 'erreur fasoar: arete ',noar,
385 % ' dans 2 triangles et a creer!'
392 c mise a jour du numero des triangles de l'arete noar
393 c le triangle 2 de l'arete => le triangle nt1
394 if( nosoar(4,noar) .lt. 0 ) then
395 c pas de triangle connu pour cette arete
398 c deja un triangle connu. ce nouveau est le second
399 if( nosoar(5,noar) .gt. 0 .and. nt1 .gt. 0 .and.
400 % nosoar(5,noar) .ne. nt1 ) then
401 c arete appartenant a plus de 2 triangles => erreur
402 write(imprim,*) 'erreur fasoar: arete ',noar,
403 % ' dans plus de 2 triangles'
411 c cas de l'arete frontaliere retrouvee comme diagonale d'un quadrangle
412 if( nt2 .gt. 0 ) then
413 c l'arete appartient a 2 triangles
414 if( nosoar(5,noar) .gt. 0 .and.
415 % nosoar(5,noar) .ne. nt2 ) then
416 c arete appartenant a plus de 2 triangles => erreur
417 write(imprim,*) 'erreur fasoar: arete ',noar,
418 % ' dans plus de 2 triangles'
431 subroutine fq1inv( x, y, s, xc, yc, ierr )
432 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
433 c but : calcul des 2 coordonnees (xc,yc) dans le carre (0,1)
434 c ----- image par f:carre unite-->quadrangle appartenant a q1**2
435 c par une resolution directe due a nicolas thenault
439 c x,y : coordonnees du point image dans le quadrangle de sommets s
440 c s : les 2 coordonnees des 4 sommets du quadrangle
444 c xc,yc : coordonnees dans le carre dont l'image par f vaut (x,y)
445 c ierr : 0 si calcul sans erreur, 1 si quadrangle degenere
446 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
447 c auteurs: thenault tulenew analyse numerique paris janvier 1998
448 c modifs : perronnet alain analyse numerique paris janvier 1998
449 c234567..............................................................012
450 real s(1:2,1:4), dist(2)
451 double precision a,b,c,d,alpha,beta,gamma,delta,x0,y0,t(2),u,v,w
456 d = s(1,1) - s(1,2) + s(1,3) - s(1,4)
459 beta = s(2,2) - s(2,1)
460 gamma = s(2,4) - s(2,1)
461 delta = s(2,1) - s(2,2) + s(2,3) - s(2,4)
463 u = beta * c - b * gamma
465 c quadrangle degenere
469 v = delta * c - d * gamma
470 w = b * delta - beta * d
472 x0 = c * (y-alpha) - gamma * (x-a)
473 y0 = b * (y-alpha) - beta * (x-a)
476 b = u * u - w * x0 - v * y0
481 delta = sqrt( b*b-4*a*c )
482 if( b .ge. 0.0 ) then
487 c la racine de plus grande valeur absolue
488 c (elle donne le plus souvent le point exterieur au carre unite
489 c donc a tester en second pour reduire les calculs)
490 t(2) = t(2) / ( 2 * a )
491 c calcul de la seconde racine a partir de la somme => plus stable
496 c la solution i donne t elle un point interne au carre unite?
497 xc = ( x0 - v * t(i) ) / u
498 yc = ( w * t(i) - y0 ) / u
499 if( 0.0 .le. xc .and. xc .le. 1.0 ) then
500 if( 0.0 .le. yc .and. yc .le. 1.0 ) goto 9000
503 c le point (xc,yc) n'est pas dans le carre unite
504 c cela peut etre du aux erreurs d'arrondi
505 c => choix par le minimum de la distance aux bords du carre
506 dist(i) = max( 0.0, -xc, xc-1.0, -yc, yc-1.0 )
510 if( dist(1) .gt. dist(2) ) then
511 c f(xc,yc) pour la racine 2 est plus proche de x,y
512 c xc yc sont deja calcules
516 else if ( b .ne. 0 ) then
522 c les 2 coordonnees du point dans le carre unite
523 xc = ( x0 - v * t(1) ) / u
524 yc = ( w * t(1) - y0 ) / u
531 subroutine ptdatr( point, pxyd, nosotr, nsigne )
532 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
533 c but : le point est il dans le triangle de sommets nosotr
538 c point : les 2 coordonnees du point
539 c pxyd : les 2 coordonnees et distance souhaitee des points du maillage
540 c nosotr : le numero des 3 sommets du triangle
544 c nsigne : >0 si le point est dans le triangle ou sur une des 3 aretes
545 c =0 si le triangle est degenere ou indirect ou ne contient pas le poin
546 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
547 c auteur : alain perronnet analyse numerique paris upmc mars 1997
548 c....................................................................012
550 double precision point(2), pxyd(3,*)
551 double precision xp,yp, x1,x2,x3, y1,y2,y3, d,dd, cb1,cb2,cb3
568 c 2 fois la surface du triangle = determinant de la matrice
569 c de calcul des coordonnees barycentriques du point p
570 d = ( x2 - x1 ) * ( y3 - y1 ) - ( x3 - x1 ) * ( y2 - y1 )
574 c triangle non degenere
575 c =====================
576 c calcul des 3 coordonnees barycentriques du
577 c point xp yp dans le triangle
578 cb1 = ( ( x2-xp ) * ( y3-yp ) - ( x3-xp ) * ( y2-yp ) ) / d
579 cb2 = ( ( x3-xp ) * ( y1-yp ) - ( x1-xp ) * ( y3-yp ) ) / d
581 ccc cb3 = ( ( x1-xp ) * ( y2-yp ) - ( x2-xp ) * ( y1-yp ) ) / d
583 ccc if( cb1 .ge. -0.00005d0 .and. cb1 .le. 1.00005d0 .and.
584 if( cb1 .ge. 0d0 .and. cb1 .le. 1d0 .and.
585 % cb2 .ge. 0d0 .and. cb2 .le. 1d0 .and.
586 % cb3 .ge. 0d0 .and. cb3 .le. 1d0 ) then
588 c le triangle nosotr contient le point
598 c le point est il du meme cote que le sommet oppose de chaque arete?
601 c le sinus de l'angle p1 p2-p1 point
604 d = ( pxyd(1,n2) - x1 ) * ( point(2) - y1 )
605 % - ( pxyd(2,n2) - y1 ) * ( point(1) - x1 )
606 dd = ( pxyd(1,n2) - x1 ) * ( pxyd(2,n3) - y1 )
607 % - ( pxyd(2,n2) - y1 ) * ( pxyd(1,n3) - x1 )
608 cb1 = ( pxyd(1,n2) - x1 ) ** 2
609 % + ( pxyd(2,n2) - y1 ) ** 2
610 cb2 = ( point(1) - x1 ) ** 2
611 % + ( point(2) - y1 ) ** 2
612 cb3 = ( pxyd(1,n3) - x1 ) ** 2
613 % + ( pxyd(2,n3) - y1 ) ** 2
614 if( abs( dd ) .le. 1e-4 * sqrt( cb1 * cb3 ) ) then
615 c le point 3 est sur l'arete 1-2
616 c le point doit y etre aussi
617 if( abs( d ) .le. 1e-4 * sqrt( cb1 * cb2 ) ) then
622 c le point 3 n'est pas sur l'arete . test des signes
623 if( d * dd .ge. 0 ) then
627 c permutation circulaire des 3 sommets et aretes
633 if( nsigne .ne. 3 ) nsigne = 0
637 integer function nosstr( p, pxyd, nt, letree )
638 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
639 c but : calculer le numero 0 a 3 du sous-triangle te contenant
644 c p : point de r**2 contenu dans le te nt de letree
645 c pxyd : x y distance des points
646 c nt : numero letree du te de te voisin a calculer
647 c letree : arbre-4 des triangles equilateraux (te) fond de la triangulation
648 c letree(0,0) no du 1-er te vide dans letree
649 c letree(0,1) : maximum du 1-er indice de letree (ici 8)
650 c letree(0,2) : maximum declare du 2-eme indice de letree (ici mxtree)
651 c letree(0:8,1) : racine de l'arbre (triangle sans sur triangle)
652 c si letree(0,.)>0 alors
653 c letree(0:3,j) : no (>0) letree des 4 sous-triangles du triangle j
655 c letree(0:3,j) :-no pxyd des 1
\85a 4 points internes au triangle j
657 c ( j est alors une feuille de l'arbre )
658 c letree(4,j) : no letree du sur-triangle du triangle j
659 c letree(5,j) : 0 1 2 3 no du sous-triangle j pour son sur-triangle
660 c letree(6:8,j) : no pxyd des 3 sommets du triangle j
664 c nosstr : 0 si le sous-triangle central contient p
665 c i =1,2,3 numero du sous-triangle contenant p
666 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
667 c auteur : alain perronnet analyse numerique paris upmc fevrier 1992
668 c2345x7..............................................................012
669 integer letree(0:8,0:*)
670 double precision pxyd(3,*), p(2),
671 % x1, y1, x21, y21, x31, y31, d, xe, ye
673 c le numero des 3 sommets du triangle
674 ns1 = letree( 6, nt )
675 ns2 = letree( 7, nt )
676 ns3 = letree( 8, nt )
678 c les coordonnees entre 0 et 1 du point p
682 x21 = pxyd(1,ns2) - x1
683 y21 = pxyd(2,ns2) - y1
685 x31 = pxyd(1,ns3) - x1
686 y31 = pxyd(2,ns3) - y1
688 d = 1.0 / ( x21 * y31 - x31 * y21 )
690 xe = ( ( p(1) - x1 ) * y31 - ( p(2) - y1 ) * x31 ) * d
691 ye = ( ( p(2) - y1 ) * x21 - ( p(1) - x1 ) * y21 ) * d
693 if( xe .gt. 0.5d0 ) then
694 c sous-triangle droit
696 else if( ye .gt. 0.5d0 ) then
699 else if( xe+ye .lt. 0.5d0 ) then
700 c sous-triangle gauche
703 c sous-triangle central
709 integer function notrpt( p, pxyd, notrde, letree )
710 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
711 c but : calculer le numero letree du sous-triangle feuille contenant
712 c ----- le point p a partir du te notrde de letree
716 c p : point de r**2 contenu dans le te nt de letree
717 c pxyd : x y distance des points
718 c notrde : numero letree du triangle depart de recherche (1=>racine)
719 c letree : arbre-4 des triangles equilateraux (te) fond de la triangulation
720 c letree(0,0) no du 1-er te vide dans letree
721 c letree(0,1) : maximum du 1-er indice de letree (ici 8)
722 c letree(0,2) : maximum declare du 2-eme indice de letree (ici mxtree)
723 c letree(0:8,1) : racine de l'arbre (triangle sans sur triangle)
724 c si letree(0,.)>0 alors
725 c letree(0:3,j) : no (>0) letree des 4 sous-triangles du triangle j
727 c letree(0:3,j) :-no pxyd des 1
\85 4 points internes au triangle j
729 c ( j est alors une feuille de l'arbre )
730 c letree(4,j) : no letree du sur-triangle du triangle j
731 c letree(5,j) : 0 1 2 3 no du sous-triangle j pour son sur-triangle
732 c letree(6:8,j) : no pxyd des 3 sommets du triangle j
736 c notrpt : numero letree du triangle contenant le point p
737 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
738 c auteur : alain perronnet analyse numerique paris upmc fevrier 1992
739 c2345x7..............................................................012
740 integer letree(0:8,0:*)
741 double precision pxyd(1:3,*), p(2)
743 c la racine depart de la recherche
746 c tant que la feuille n'est pas atteinte descendre l'arbre
747 10 if( letree(0,notrpt) .gt. 0 ) then
749 c recherche du sous-triangle contenant p
750 nsot = nosstr( p, pxyd, notrpt, letree )
752 c le numero letree du sous-triangle
753 notrpt = letree( nsot, notrpt )
760 subroutine teajpt( ns, nbsomm, mxsomm, pxyd, letree,
762 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
763 c but : ajout du point ns de pxyd dans letree
768 c ns : numero du point a ajouter dans letree
769 c mxsomm : nombre maximal de points declarables dans pxyd
770 c pxyd : tableau des coordonnees des points
771 c par point : x y distance_souhaitee
775 c nbsomm : nombre actuel de points dans pxyd
777 c letree : arbre-4 des triangles equilateraux (te) fond de la triangulation
778 c letree(0,0) : no du 1-er te vide dans letree
779 c letree(0,1) : maximum du 1-er indice de letree (ici 8)
780 c letree(0,2) : maximum declare du 2-eme indice de letree (ici mxtree)
781 c letree(0:8,1) : racine de l'arbre (triangle sans sur triangle)
782 c si letree(0,.)>0 alors
783 c letree(0:3,j) : no (>0) letree des 4 sous-triangles du triangle j
785 c letree(0:3,j) :-no pxyd des 1
\85a 4 points internes au triangle j
787 c ( j est alors une feuille de l'arbre )
788 c letree(4,j) : no letree du sur-triangle du triangle j
789 c letree(5,j) : 0 1 2 3 no du sous-triangle j pour son sur-triangle
790 c letree(6:8,j) : no pxyd des 3 sommets du triangle j
794 c ntrp : numero letree du triangle te ou a ete ajoute le point
795 c ierr : 0 si pas d'erreur, 51 saturation letree, 52 saturation pxyd
796 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
797 c auteur : alain perronnet analyse numerique paris upmc fevrier 1992
798 c2345x7..............................................................012
799 integer letree(0:8,0:*)
800 double precision pxyd(3,mxsomm)
802 c depart de la racine
805 c recherche du triangle contenant le point pxyd(ns)
806 1 ntrp = notrpt( pxyd(1,ns), pxyd, ntrp, letree )
808 c existe t il un point libre
810 if( letree(i,ntrp) .eq. 0 ) then
811 c la place i est libre
817 c pas de place libre => 4 sous-triangles sont crees
818 c a partir des 3 milieux des aretes
819 call te4ste( nbsomm, mxsomm, pxyd, ntrp, letree, ierr )
820 if( ierr .ne. 0 ) return
826 subroutine n1trva( nt, lar, letree, notrva, lhpile )
827 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
828 c but : calculer le numero letree du triangle voisin du te nt
829 c ----- par l'arete lar (1 a 3 ) de nt
830 c attention : notrva n'est pas forcement minimal
834 c nt : numero letree du te de te voisin a calculer
835 c lar : numero 1 a 3 de l'arete du triangle nt
836 c letree : arbre-4 des triangles equilateraux (te) fond de la triangulation
837 c letree(0,0) no du 1-er te vide dans letree
838 c letree(0,1) : maximum du 1-er indice de letree (ici 8)
839 c letree(0,2) : maximum declare du 2-eme indice de letree (ici mxtree)
840 c letree(0:8,1) : racine de l'arbre (triangle sans sur-triangle)
841 c si letree(0,.)>0 alors
842 c letree(0:3,j) : no (>0) letree des 4 sous-triangles du triangle j
844 c letree(0:3,j) :-no pxyd des 1 a 4 points internes au triangle j
846 c ( j est alors une feuille de l'arbre )
847 c letree(4,j) : no letree du sur-triangle du triangle j
848 c letree(5,j) : 0 1 2 3 no du sous-triangle j pour son sur-triangle
849 c letree(6:8,j) : no pxyd des 3 sommets du triangle j
853 c notrva : >0 numero letree du te voisin par l'arete lar
854 c =0 si pas de te voisin (racine , ... )
855 c lhpile : =0 si nt et notrva ont meme taille
856 c >0 nt est 4**lhpile fois plus petit que notrva
857 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
858 c auteur : alain perronnet analyse numerique paris upmc fevrier 1992
859 c2345x7..............................................................012
860 integer letree(0:8,0:*)
863 c initialisation de la pile
864 c le triangle est empile
868 c tant qu'il existe un sur-triangle
869 10 ntr = lapile( lhpile )
870 if( ntr .eq. 1 ) then
871 c racine atteinte => pas de triangle voisin
877 c le type du triangle ntr
878 nty = letree( 5, ntr )
879 c l'eventuel sur-triangle
880 nsut = letree( 4, ntr )
882 if( nty .eq. 0 ) then
884 c triangle de type 0 => triangle voisin de type precedent(lar)
885 c dans le sur-triangle de ntr
886 c ce triangle remplace ntr dans lapile
887 lapile( lhpile ) = letree( nopre3(lar), nsut )
891 c triangle ntr de type nty>0
892 if( nosui3(nty) .eq. lar ) then
894 c le triangle voisin par lar est le triangle 0
895 lapile( lhpile ) = letree( 0, nsut )
899 c triangle sans voisin direct => passage par le sur-triangle
900 if( nsut .eq. 0 ) then
902 c ntr est la racine => pas de triangle voisin par cette arete
907 c le sur-triangle est empile
909 lapile(lhpile) = nsut
913 c descente aux sous-triangles selon la meme arete
914 20 notrva = lapile( lhpile )
916 30 lhpile = lhpile - 1
917 if( letree(0,notrva) .le. 0 ) then
918 c le triangle est une feuille de l'arbre 0 sous-triangle
919 c lhpile = nombre de differences de niveaux dans l'arbre
922 c le triangle a 4 sous-triangles
923 if( lhpile .gt. 0 ) then
925 c bas de pile non atteint
926 nty = letree( 5, lapile(lhpile) )
927 if( nty .eq. lar ) then
928 c l'oppose est suivant(nty) de notrva
929 notrva = letree( nosui3(nty) , notrva )
931 c l'oppose est precedent(nty) de notrva
932 notrva = letree( nopre3(nty) , notrva )
938 c meme niveau dans l'arbre lhpile = 0
942 subroutine cenced( xy1, xy2, xy3, cetria, ierr )
943 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
944 c but : calcul des coordonnees du centre du cercle circonscrit
945 c ----- du triangle defini par ses 3 sommets de coordonnees
946 c xy1 xy2 xy3 ainsi que le carre du rayon de ce cercle
950 c xy1 xy2 xy3 : les 2 coordonnees des 3 sommets du triangle
951 c ierr : <0 => pas d'affichage si triangle degenere
952 c >=0 => affichage si triangle degenere
956 c cetria : cetria(1)=abcisse du centre
957 c cetria(2)=ordonnee du centre
958 c cetria(3)=carre du rayon 1d28 si triangle degenere
959 c ierr : 0 si triangle non degenere
960 c 1 si triangle degenere
961 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
962 c auteur : perronnet alain upmc analyse numerique paris juin 1995
963 c2345x7..............................................................012
964 parameter (epsurf=1d-7)
965 common / unites / lecteu,imprim,nunite(30)
966 double precision x1,y1,x21,y21,x31,y31,
968 % xy1(2),xy2(2),xy3(2),cetria(3)
970 c le calcul de 2 fois l'aire du triangle
971 c attention l'ordre des 3 sommets est direct ou non
980 aire2 = x21 * y31 - x31 * y21
982 c recherche d'un test relatif peu couteux
983 c pour reperer la degenerescence du triangle
985 % epsurf*(abs(x21)+abs(x31))*(abs(y21)+abs(y31)) ) then
986 c triangle de qualite trop faible
987 if( ierr .ge. 0 ) then
989 c kerr(1) = 'erreur cenced: triangle degenere'
991 write(imprim,*) 'erreur cenced: triangle degenere'
992 write(imprim,10000) xy1,xy2,xy3,aire2
994 10000 format( 3(' x=',g24.16,' y=',g24.16/),' aire*2=',g24.16)
1002 c les 2 coordonnees du centre intersection des 2 mediatrices
1003 c x = (x1+x2)/2 + lambda * (y2-y1)
1004 c y = (y1+y2)/2 - lambda * (x2-x1)
1005 c x = (x1+x3)/2 + rot * (y3-y1)
1006 c y = (y1+y3)/2 - rot * (x3-x1)
1007 c ==========================================================
1008 rot = ((xy2(1)-xy3(1))*x21 + (xy2(2)-xy3(2))*y21) / (2 * aire2)
1010 xc = ( x1 + xy3(1) ) * 0.5d0 + rot * y31
1011 yc = ( y1 + xy3(2) ) * 0.5d0 - rot * x31
1017 cetria(3) = (x1-xc) ** 2 + (y1-yc) ** 2
1019 c pas d'erreur rencontree
1024 double precision function angled( p1, p2, p3 )
1025 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1026 c but : calculer l'angle (p1p2,p1p3) en radians
1031 c p1,p2,p3 : les 2 coordonnees des 3 sommets de l'angle
1032 c sens direct pour une surface >0
1035 c angled : angle (p1p2,p1p3) en radians entre [0 et 2pi]
1036 c 0 si p1=p2 ou p1=p3
1037 c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1038 c auteur : alain perronnet analyse numerique upmc paris fevrier 1992
1039 c2345x7..............................................................012
1040 double precision p1(2),p2(2),p3(2),x21,y21,x31,y31,a1,a2,d,c
1048 c longueur des cotes
1049 a1 = x21 * x21 + y21 * y21
1050 a2 = x31 * x31 + y31 * y31
1057 c cosinus de l'angle
1058 c = ( x21 * x31 + y21 * y31 ) / d
1059 if( c .le. -1.d0 ) then
1060 c tilt sur apollo si acos( -1 -eps )
1061 angled = atan( 1.d0 ) * 4.d0
1063 else if( c .ge. 1.d0 ) then
1064 c tilt sur apollo si acos( 1 + eps )
1070 if( x21 * y31 - x31 * y21 .lt. 0 ) then
1071 c demi plan inferieur
1072 angled = 8.d0 * atan( 1.d0 ) - angled
1077 subroutine teajte( mxsomm, nbsomm, pxyd, comxmi,
1078 % aretmx, mxtree, letree,
1080 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1081 c but : initialisation des tableaux letree
1082 c ----- ajout des sommets 1 a nbsomm (valeur en entree) dans letree
1086 c mxsomm : nombre maximal de sommets permis pour la triangulation
1087 c mxtree : nombre maximal de triangles equilateraux (te) declarables
1088 c aretmx : longueur maximale des aretes des triangles equilateraux
1090 c entrees et sorties :
1091 c --------------------
1092 c nbsomm : nombre de sommets apres identification
1093 c pxyd : tableau des coordonnees 2d des points
1094 c par point : x y distance_souhaitee
1095 c tableau reel(3,mxsomm)
1099 c comxmi : coordonnees minimales et maximales des points frontaliers
1100 c letree : arbre-4 des triangles equilateraux (te) fond de la triangulation
1101 c letree(0,0) : no du 1-er te vide dans letree
1102 c letree(0,1) : maximum du 1-er indice de letree (ici 8)
1103 c letree(0,2) : maximum declare du 2-eme indice de letree (ici mxtree)
1104 c letree(0:8,1) : racine de l'arbre (triangle sans sur triangle)
1105 c si letree(0,.)>0 alors
1106 c letree(0:3,j) : no (>0) letree des 4 sous-triangles du triangle j
1108 c letree(0:3,j) :-no pxyd des 1 a 4 points internes au triangle j
1110 c ( j est alors une feuille de l'arbre )
1111 c letree(4,j) : no letree du sur-triangle du triangle j
1112 c letree(5,j) : 0 1 2 3 no du sous-triangle j pour son sur-triangle
1113 c letree(6:8,j) : no pxyd des 3 sommets du triangle j
1115 c ierr : 0 si pas d'erreur
1116 c 51 saturation letree
1117 c 52 saturation pxyd
1118 c 7 tous les points sont alignes
1119 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1120 c auteur : alain perronnet analyse numerique paris upmc juillet 1994
1121 c....................................................................012
1122 integer letree(0:8,0:mxtree)
1123 double precision pxyd(3,mxsomm)
1124 double precision comxmi(3,2)
1125 double precision a(2),s,aretmx,rac3
1127 c protection du nombre de sommets avant d'ajouter ceux de tetree
1130 comxmi(1,1) = min( comxmi(1,1), pxyd(1,i) )
1131 comxmi(1,2) = max( comxmi(1,2), pxyd(1,i) )
1132 comxmi(2,1) = min( comxmi(2,1), pxyd(2,i) )
1133 comxmi(2,2) = max( comxmi(2,2), pxyd(2,i) )
1136 c creation de l'arbre tee
1137 c =======================
1138 c la premiere colonne vide de letree
1140 c chainage des te vides
1144 letree(0,mxtree) = 0
1145 c les maxima des 2 indices de letree
1147 letree(2,0) = mxtree
1150 c aucun point interne au triangle equilateral (te) 1
1155 c pas de sur-triangle
1158 c le numero pxyd des 3 sommets du te 1
1159 letree(6,1) = nbsomm + 1
1160 letree(7,1) = nbsomm + 2
1161 letree(8,1) = nbsomm + 3
1163 c calcul de la largeur et hauteur du rectangle englobant
1164 c ======================================================
1165 a(1) = comxmi(1,2) - comxmi(1,1)
1166 a(2) = comxmi(2,2) - comxmi(2,1)
1167 c la longueur de la diagonale
1168 s = sqrt( a(1)**2 + a(2)**2 )
1170 if( a(k) .lt. 1e-4 * s ) then
1172 write(imprim,*) 'tous les points sont alignes'
1179 c le maximum des ecarts
1182 c le triangle equilateral englobant
1183 c =================================
1184 c ecart du rectangle au triangle equilateral
1185 rac3 = sqrt( 3.0d0 )
1186 arete = a(1) + 2 * aretmx + 2 * ( a(2) + aretmx ) / rac3
1188 c le point nbsomm + 1 en bas a gauche
1190 pxyd(1,nbsomm) = (comxmi(1,1)+comxmi(1,2))*0.5d0 - arete*0.5d0
1191 pxyd(2,nbsomm) = comxmi(2,1) - aretmx
1194 c le point nbsomm + 2 en bas a droite
1196 pxyd(1,nbsomm) = pxyd(1,nbsomm-1) + arete
1197 pxyd(2,nbsomm) = pxyd(2,nbsomm-1)
1200 c le point nbsomm + 3 sommet au dessus
1202 pxyd(1,nbsomm) = pxyd(1,nbsomm-2) + arete * 0.5d0
1203 pxyd(2,nbsomm) = pxyd(2,nbsomm-2) + arete * 0.5d0 * rac3
1206 c ajout des sommets des lignes pour former letree
1207 c ===============================================
1209 c ajout du point i de pxyd a letree
1210 call teajpt( i, nbsomm, mxsomm, pxyd, letree,
1212 if( ierr .ne. 0 ) return
1219 subroutine tetaid( nutysu, dx, dy, longai, ierr )
1220 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1221 c but : calculer la longueur de l'arete ideale en dx,dy
1226 c nutysu : numero de traitement de areteideale() selon le type de surface
1227 c 0 pas d'emploi de la fonction areteideale() => aretmx active
1228 c 1 il existe une fonction areteideale(xyz,xyzdir)
1229 c ... autres options a definir ...
1230 c dx, dy : abscisse et ordonnee dans le plan du point (reel2!)
1234 c longai : longueur de l'areteideale(xyz,xyzdir) autour du point xyz
1235 c ierr : 0 si pas d'erreur, <>0 sinon
1236 c 1 calcul incorrect de areteideale(xyz,xyzdir)
1237 c 2 longueur calculee nulle
1238 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1239 c auteur : alain perronnet analyse numerique paris upmc mars 1997
1240 c2345x7..............................................................012
1241 common / unites / lecteu, imprim, nunite(30)
1243 double precision areteideale
1244 double precision dx, dy, longai
1245 double precision xyz(3), xyzd(3), d0
1248 if( nutysu .gt. 0 ) then
1250 c le point ou se calcule la longueur
1253 c z pour le calcul de la longueur (inactif ici!)
1255 c la direction pour le calcul de la longueur (inactif ici!)
1260 longai = areteideale(xyz,xyzd)
1261 if( longai .lt. 0d0 ) then
1262 write(imprim,10000) xyz
1263 10000 format('attention: longueur de areteideale(',
1264 % g14.6,',',g14.6,',',g14.6,')<=0! => rendue >0' )
1267 if( longai .eq. 0d0 ) then
1268 write(imprim,10001) xyz
1269 10001 format('erreur: longueur de areteideale(',
1270 % g14.6,',',g14.6,',',g14.6,')=0!' )
1278 subroutine tehote( nutysu,
1279 % nbarpi, mxsomm, nbsomm, pxyd,
1281 % letree, mxqueu, laqueu,
1283 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1284 c but : homogeneisation de l'arbre des te a un saut de taille au plus
1285 c ----- prise en compte des distances souhaitees autour des sommets initiaux
1289 c nutysu : numero de traitement de areteideale() selon le type de surface
1290 c 0 pas d'emploi de la fonction areteideale() => aretmx active
1291 c 1 il existe une fonction areteideale()
1292 c dont seules les 2 premieres composantes de uv sont actives
1293 c autres options a definir...
1294 c nbarpi : nombre de sommets de la frontiere + nombre de points internes
1295 c imposes par l'utilisateur
1296 c mxsomm : nombre maximal de sommets permis pour la triangulation et te
1297 c mxqueu : nombre d'entiers utilisables dans laqueu
1298 c comxmi : minimum et maximum des coordonnees de l'objet
1299 c aretmx : longueur maximale des aretes des triangles equilateraux
1300 c permtr : perimetre de la ligne enveloppe dans le plan
1301 c avant mise a l'echelle a 2**20
1305 c nbsomm : nombre de sommets apres identification
1306 c pxyd : tableau des coordonnees 2d des points
1307 c par point : x y distance_souhaitee
1308 c letree : arbre-4 des triangles equilateraux (te) fond de la triangulation
1309 c letree(0,0) : no du 1-er te vide dans letree
1310 c letree(1,0) : maximum du 1-er indice de letree (ici 8)
1311 c letree(2,0) : maximum declare du 2-eme indice de letree (ici mxtree)
1312 c letree(0:8,1) : racine de l'arbre (triangle sans sur triangle)
1313 c si letree(0,.)>0 alors
1314 c letree(0:3,j) : no (>0) letree des 4 sous-triangles du triangle j
1316 c letree(0:3,j) :-no pxyd des 1 a 4 points internes au triangle j
1318 c ( j est alors une feuille de l'arbre )
1319 c letree(4,j) : no letree du sur-triangle du triangle j
1320 c letree(5,j) : 0 1 2 3 no du sous-triangle j pour son sur-triangle
1321 c letree(6:8,j) : no pxyd des 3 sommets du triangle j
1325 c laqueu : mxqueu entiers servant de queue pour le parcours de letree
1329 c ierr : 0 si pas d'erreur
1330 c 51 si saturation letree dans te4ste
1331 c 52 si saturation pxyd dans te4ste
1332 c >0 si autre erreur
1333 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1334 c auteur : alain perronnet analyse numerique paris upmc avril 1997
1335 c2345x7..............................................................012
1336 double precision ampli
1337 parameter (ampli=1.34d0)
1338 common / unites / lecteu, imprim, intera, nunite(29)
1340 double precision pxyd(3,mxsomm), d2, aretm2
1341 double precision comxmi(3,2),aretmx,a,s,xrmin,xrmax,yrmin,yrmax
1342 double precision dmin, dmax
1343 integer letree(0:8,0:*)
1345 integer laqueu(1:mxqueu),lequeu
1346 c lequeu : entree dans la queue
1347 c lhqueu : longueur de la queue
1348 c gestion circulaire
1351 equivalence (nuste(1),ns1),(nuste(2),ns2),(nuste(3),ns3)
1353 c existence ou non de la fonction 'taille_ideale' des aretes
1354 c autour du point. ici la carte est supposee isotrope
1355 c ==========================================================
1356 c attention: si la fonction taille_ideale existe
1357 c alors pxyd(3,*) est la taille_ideale dans l'espace initial
1358 c sinon pxyd(3,*) est la distance calculee dans le plan par
1359 c propagation a partir des tailles des aretes de la frontiere
1361 if( nutysu .gt. 0 ) then
1363 c la fonction taille_ideale(x,y,z) existe
1364 c ---------------------------------------
1365 c initialisation de la distance souhaitee autour des points 1 a nbsomm
1367 c calcul de pxyzd(3,i)
1368 call tetaid( nutysu, pxyd(1,i), pxyd(2,i),
1370 if( ierr .ne. 0 ) goto 9999
1375 c la fonction taille_ideale(x,y,z) n'existe pas
1376 c ---------------------------------------------
1377 c prise en compte des distances souhaitees dans le plan
1378 c autour des points frontaliers et des points internes imposes
1379 c toutes les autres distances souhaitees ont ete mis a aretmx
1380 c lors de l'execution du sp teqini
1382 c le sommet i n'est pas un sommet de letree => sommet frontalier
1383 c recherche du sous-triangle minimal feuille contenant le point i
1385 2 nte = notrpt( pxyd(1,i), pxyd, nte, letree )
1386 c la distance au sommet le plus eloigne est elle inferieure
1387 c a la distance souhaitee?
1391 d2 = max( ( pxyd(1,i)-pxyd(1,ns1) )**2 +
1392 % ( pxyd(2,i)-pxyd(2,ns1) )**2
1393 % , ( pxyd(1,i)-pxyd(1,ns2) )**2 +
1394 % ( pxyd(2,i)-pxyd(2,ns2) )**2
1395 % , ( pxyd(1,i)-pxyd(1,ns3) )**2 +
1396 % ( pxyd(2,i)-pxyd(2,ns3) )**2 )
1397 if( d2 .gt. pxyd(3,i)**2 ) then
1398 c le triangle nte trop grand doit etre subdivise en 4 sous-triangle
1399 call te4ste( nbsomm, mxsomm, pxyd, nte, letree,
1401 if( ierr .ne. 0 ) return
1407 c le sous-triangle central de la racine est decoupe systematiquement
1408 c ==================================================================
1410 if( letree(0,2) .le. 0 ) then
1411 c le sous-triangle central de la racine n'est pas subdivise
1412 c il est donc decoupe en 4 soustriangles
1414 call te4ste( nbsomm, mxsomm, pxyd, nte, letree,
1416 if( ierr .ne. 0 ) return
1417 do 4 i=nbsom0+1,nbsomm
1418 c mise a jour de taille_ideale des nouveaux sommets de te
1419 call tetaid( nutysu, pxyd(1,i), pxyd(2,i), pxyd(3,i), ierr )
1420 if( ierr .ne. 0 ) goto 9999
1424 c le carre de la longueur de l'arete de triangles equilateraux
1425 c souhaitee pour le fond de la triangulation
1426 aretm2 = (aretmx*ampli) ** 2
1428 c tout te contenu dans le rectangle englobant doit avoir un
1429 c cote < aretmx et etre de meme taille que les te voisins
1430 c s'il contient un point; sinon un seul saut de taille est permis
1431 c ===============================================================
1432 c le rectangle englobant pour selectionner les te "internes"
1433 c le numero des 3 sommets du te englobant racine de l'arbre des te
1438 c abscisse du milieu de l'arete gauche du te 1
1439 s = ( pxyd(1,ns1) + pxyd(1,ns3) ) / 2
1440 xrmin = min( s, comxmi(1,1) - aretmx ) - a
1441 c abscisse du milieu de l'arete droite du te 1
1442 s = ( pxyd(1,ns2) + pxyd(1,ns3) ) / 2
1443 xrmax = max( s, comxmi(1,2) + aretmx ) + a
1444 yrmin = comxmi(2,1) - aretmx
1445 c ordonnee de la droite passant par les milieus des 2 aretes
1446 c droite gauche du te 1
1447 s = ( pxyd(2,ns1) + pxyd(2,ns3) ) / 2
1448 yrmax = max( s, comxmi(2,2) + aretmx ) + a
1450 c cas particulier de 3 ou 4 ou peu d'aretes frontalieres
1451 if( nbarpi .le. 8 ) then
1452 c tout le triangle englobant (racine) est a prendre en compte
1453 xrmin = pxyd(1,ns1) - a
1454 xrmax = pxyd(1,ns2) + a
1455 yrmin = pxyd(2,ns1) - a
1456 yrmax = pxyd(2,ns3) + a
1462 c initialisation de la queue
1463 5 nbiter = nbiter + 1
1466 c la racine de letree initialise la queue
1469 c tant que la longueur de la queue est >=0 traiter le debut de queue
1470 10 if( lhqueu .ge. 0 ) then
1472 c le triangle te a traiter
1474 if( i .le. 0 ) i = mxqueu + i
1476 c la longueur de la queue est reduite
1479 c nte est il un sous-triangle feuille minimal ?
1480 15 if( letree(0,nte) .gt. 0 ) then
1482 c non les 4 sous-triangles sont mis dans la queue
1483 if( lhqueu + 4 .ge. mxqueu ) then
1484 write(imprim,*) 'tehote: saturation de la queue'
1489 c ajout du sous-triangle i
1492 if( lequeu .gt. mxqueu ) lequeu = lequeu - mxqueu
1493 laqueu( lequeu ) = letree( i, nte )
1499 c ici nte est un triangle minimal non subdivise
1500 c ---------------------------------------------
1501 c le te est il dans le cadre englobant de l'objet ?
1505 if( pxyd(1,ns1) .gt. pxyd(1,ns2) ) then
1512 if( (xrmin .le. dmin .and. dmin .le. xrmax) .or.
1513 % (xrmin .le. dmax .and. dmax .le. xrmax) ) then
1514 if( pxyd(2,ns1) .gt. pxyd(2,ns3) ) then
1521 if( (yrmin .le. dmin .and. dmin .le. yrmax) .or.
1522 % (yrmin .le. dmax .and. dmax .le. yrmax) ) then
1524 c nte est un te feuille et interne au rectangle englobant
1525 c =======================================================
1526 c le carre de la longueur de l'arete du te de numero nte
1527 d2 = (pxyd(1,ns1)-pxyd(1,ns2)) ** 2 +
1528 % (pxyd(2,ns1)-pxyd(2,ns2)) ** 2
1530 if( nutysu .eq. 0 ) then
1532 c il n'existe pas de fonction 'taille_ideale'
1533 c -------------------------------------------
1534 c si la taille effective de l'arete du te est superieure a aretmx
1535 c alors le te est decoupe
1536 if( d2 .gt. aretm2 ) then
1537 c le triangle nte trop grand doit etre subdivise
1538 c en 4 sous-triangles
1539 call te4ste( nbsomm,mxsomm, pxyd,
1540 % nte, letree, ierr )
1541 if( ierr .ne. 0 ) return
1547 c il existe ici une fonction 'taille_ideale'
1548 c ------------------------------------------
1549 c si la taille effective de l'arete du te est superieure au mini
1550 c des 3 tailles_ideales aux sommets alors le te est decoupe
1552 if( d2 .gt. (pxyd(3,nuste(i))*ampli)**2 ) then
1553 c le triangle nte trop grand doit etre subdivise
1554 c en 4 sous-triangles
1556 call te4ste( nbsomm, mxsomm, pxyd,
1557 & nte, letree, ierr )
1558 if( ierr .ne. 0 ) return
1559 do 27 j=nbsom0+1,nbsomm
1560 c mise a jour de taille_ideale des nouveaux sommets de
1561 call tetaid( nutysu, pxyd(1,j), pxyd(2,j),
1563 if( ierr .ne. 0 ) goto 9999
1570 c recherche du nombre de niveaux entre nte et les te voisins par se
1571 c si la difference de subdivisions excede 1 alors le plus grand des
1572 c =================================================================
1575 c noteva triangle voisin de nte par l'arete i
1576 call n1trva( nte, i, letree, noteva, niveau )
1577 if( noteva .le. 0 ) goto 30
1578 c il existe un te voisin
1579 if( niveau .gt. 0 ) goto 30
1580 c nte a un te voisin plus petit ou egal
1581 if( letree(0,noteva) .le. 0 ) goto 30
1582 c nte a un te voisin noteva subdivise au moins une fois
1584 if( nbiter .gt. 0 ) then
1585 c les 2 sous triangles voisins sont-ils subdivises?
1586 ns2 = letree(i,noteva)
1587 if( letree(0,ns2) .le. 0 ) then
1588 c ns2 n'est pas subdivise
1589 ns2 = letree(nosui3(i),noteva)
1590 if( letree(0,ns2) .le. 0 ) then
1591 c les 2 sous-triangles ne sont pas subdivises
1597 c saut>1 => le triangle nte doit etre subdivise en 4 sous-triang
1598 c --------------------------------------------------------------
1600 call te4ste( nbsomm,mxsomm, pxyd, nte, letree,
1602 if( ierr .ne. 0 ) return
1603 if( nutysu .gt. 0 ) then
1604 do 32 j=nbsom0+1,nbsomm
1605 c mise a jour de taille_ideale des nouveaux sommets de te
1606 call tetaid( nutysu, pxyd(1,j), pxyd(2,j),
1608 if( ierr .ne. 0 ) goto 9999
1618 if( nbs0 .lt. nbsomm ) then
1624 c pb dans le calcul de la fonction taille_ideale
1626 9999 write(imprim,*) 'pb dans le calcul de taille_ideale'
1628 c kerr(1) = 'pb dans le calcul de taille_ideale'
1634 subroutine tetrte( comxmi, aretmx, nbarpi, mxsomm, pxyd,
1635 % mxqueu, laqueu, letree,
1636 % mosoar, mxsoar, n1soar, nosoar,
1637 % moartr, mxartr, n1artr, noartr, noarst,
1639 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1640 c but : trianguler les triangles equilateraux feuilles et
1641 c ----- les points de la frontiere et les points internes imposes
1643 c attention: la triangulation finale n'est pas de type delaunay!
1647 c comxmi : minimum et maximum des coordonnees de l'objet
1648 c aretmx : longueur maximale des aretes des triangles equilateraux
1649 c nbarpi : nombre de sommets de la frontiere + nombre de points internes
1650 c imposes par l'utilisateur
1651 c mxsomm : nombre maximal de sommets declarables dans pxyd
1652 c pxyd : tableau des coordonnees 2d des points
1653 c par point : x y distance_souhaitee
1655 c mxqueu : nombre d'entiers utilisables dans laqueu
1656 c mosoar : nombre maximal d'entiers par arete du tableau nosoar
1657 c mxsoar : nombre maximal d'aretes stockables dans le tableau nosoar
1658 c moartr : nombre maximal d'entiers par arete du tableau noartr
1659 c mxartr : nombre maximal de triangles stockables dans le tableau noartr
1660 c letree : arbre-4 des triangles equilateraux (te) fond de la triangulation
1661 c letree(0,0) : no du 1-er te vide dans letree
1662 c letree(0,1) : maximum du 1-er indice de letree (ici 8)
1663 c letree(0,2) : maximum declare du 2-eme indice de letree (ici mxtree)
1664 c letree(0:8,1) : racine de l'arbre (triangle sans sur triangle)
1665 c si letree(0,.)>0 alors
1666 c letree(0:3,j) : no (>0) letree des 4 sous-triangles du triangle j
1668 c letree(0:3,j) :-no pxyd des 1 a 4 points internes au triangle j
1670 c ( j est alors une feuille de l'arbre )
1671 c letree(4,j) : no letree du sur-triangle du triangle j
1672 c letree(5,j) : 0 1 2 3 no du sous-triangle j pour son sur-triangle
1673 c letree(6:8,j) : no pxyd des 3 sommets du triangle j
1677 c n1soar : numero de la premiere arete vide dans le tableau nosoar
1678 c une arete i de nosoar est vide <=> nosoar(1,i)=0
1679 c nosoar : numero des 2 sommets , no ligne, 2 triangles de l'arete,
1680 c chainage des aretes frontalieres, chainage du hachage des aretes
1681 c hachage des aretes = nosoar(1)+nosoar(2)*2
1682 c noarst : noarst(i) numero d'une arete de sommet i
1686 c laqueu : mxqueu entiers servant de queue pour le parcours de letree
1690 c n1artr : numero du premier triangle vide dans le tableau noartr
1691 c le chainage des triangles vides se fait sur noartr(2,.)
1692 c noartr : les 3 aretes des triangles +-arete1, +-arete2, +-arete3
1693 c arete1 = 0 si triangle vide => arete2 = triangle vide suivant
1694 c ierr : =0 si pas d'erreur
1695 c =1 si le tableau nosoar est sature
1696 c =2 si le tableau noartr est sature
1697 c =3 si aucun des triangles ne contient l'un des points internes d'un t
1698 c =5 si saturation de la queue de parcours de l'arbre des te
1699 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1700 c auteur : alain perronnet analyse numerique paris upmc mars 1997
1701 c2345x7..............................................................012
1702 common / unites / lecteu, imprim, intera, nunite(29)
1704 double precision pxyd(3,mxsomm)
1705 double precision comxmi(3,2),aretmx,a,s,xrmin,xrmax,yrmin,yrmax
1706 double precision dmin, dmax
1708 integer nosoar(mosoar,mxsoar),
1709 % noartr(moartr,mxartr),
1712 integer letree(0:8,0:*)
1713 integer laqueu(1:mxqueu)
1714 c lequeu:entree dans la queue en gestion circulaire
1715 c lhqueu:longueur de la queue en gestion circulaire
1717 integer milieu(3), nutr(1:13)
1719 c le rectangle englobant pour selectionner les te "internes"
1720 c le numero des 3 sommets du te englobant racine de l'arbre des te
1725 c abscisse du milieu de l'arete gauche du te 1
1726 s = ( pxyd(1,ns1) + pxyd(1,ns3) ) / 2
1727 xrmin = min( s, comxmi(1,1) - aretmx ) - a
1728 c abscisse du milieu de l'arete droite du te 1
1729 s = ( pxyd(1,ns2) + pxyd(1,ns3) ) / 2
1730 xrmax = max( s, comxmi(1,2) + aretmx ) + a
1731 yrmin = comxmi(2,1) - aretmx
1732 c ordonnee de la droite passant par les milieus des 2 aretes
1733 c droite gauche du te 1
1734 s = ( pxyd(2,ns1) + pxyd(2,ns3) ) / 2
1735 yrmax = max( s, comxmi(2,2) + aretmx ) + a
1737 c cas particulier de 3 ou 4 ou peu d'aretes frontalieres
1738 if( nbarpi .le. 8 ) then
1739 c tout le triangle englobant (racine) est a prendre en compte
1740 xrmin = pxyd(1,ns1) - a
1741 xrmax = pxyd(1,ns2) + a
1742 yrmin = pxyd(2,ns1) - a
1743 yrmax = pxyd(2,ns3) + a
1746 c initialisation du tableau noartr
1748 c le numero de l'arete est inconnu
1750 c le chainage sur le triangle vide suivant
1753 noartr(2,mxartr) = 0
1756 c parcours des te jusqu'a trianguler toutes les feuilles (triangles eq)
1757 c =====================================================================
1758 c initialisation de la queue sur les te
1762 c la racine de letree initialise la queue
1765 c tant que la longueur de la queue est >=0 traiter le debut de queue
1766 10 if( lhqueu .ge. 0 ) then
1768 c le triangle te a traiter
1770 if( i .le. 0 ) i = mxqueu + i
1772 c la longueur est reduite
1775 c nte est il un sous-triangle feuille (minimal) ?
1776 15 if( letree(0,nte) .gt. 0 ) then
1777 c non les 4 sous-triangles sont mis dans la queue
1778 if( lhqueu + 4 .ge. mxqueu ) then
1779 write(imprim,*) 'tetrte: saturation de la queue'
1784 c ajout du sous-triangle i
1787 if( lequeu .gt. mxqueu ) lequeu = lequeu - mxqueu
1788 laqueu( lequeu ) = letree( i, nte )
1793 c ici nte est un triangle minimal non subdivise
1794 c ---------------------------------------------
1795 c le te est il dans le cadre englobant de l'objet ?
1799 if( pxyd(1,ns1) .gt. pxyd(1,ns2) ) then
1806 if( (xrmin .le. dmin .and. dmin .le. xrmax) .or.
1807 % (xrmin .le. dmax .and. dmax .le. xrmax) ) then
1808 if( pxyd(2,ns1) .gt. pxyd(2,ns3) ) then
1815 if( (yrmin .le. dmin .and. dmin .le. yrmax) .or.
1816 % (yrmin .le. dmax .and. dmax .le. yrmax) ) then
1818 c te minimal et interne au rectangle englobant
1819 c --------------------------------------------
1820 c recherche du nombre de niveaux entre nte et les te voisins
1825 c a priori pas de milieu de l'arete i du te nte
1828 c recherche de noteva te voisin de nte par l'arete i
1829 call n1trva( nte, i, letree, noteva, niveau )
1830 c noteva : >0 numero letree du te voisin par l'arete i
1831 c =0 si pas de te voisin (racine , ... )
1832 c niveau : =0 si nte et noteva ont meme taille
1833 c >0 nte est 4**niveau fois plus petit que noteva
1834 if( noteva .gt. 0 ) then
1835 c il existe un te voisin
1836 if( letree(0,noteva) .gt. 0 ) then
1837 c noteva est plus petit que nte
1838 c => recherche du numero du milieu du cote=sommet du te no
1839 c le sous-te 0 du te noteva
1840 nsot = letree(0,noteva)
1841 c le numero dans pxyd du milieu de l'arete i de nte
1842 milieu( i ) = letree( 5+nopre3(i), nsot )
1849 c triangulation du te nte en fonction du nombre de ses milieux
1850 goto( 50, 100, 200, 300 ) , nbmili + 1
1852 c 0 milieu => 1 triangle = le te nte
1853 c ----------------------------------
1854 50 call f0trte( letree(0,nte), pxyd,
1855 % mosoar, mxsoar, n1soar, nosoar,
1856 % moartr, mxartr, n1artr, noartr,
1858 % nbtr, nutr, ierr )
1859 if( ierr .ne. 0 ) return
1862 c 1 milieu => 2 triangles = 2 demi te
1863 c -----------------------------------
1864 100 call f1trte( letree(0,nte), pxyd, milieu,
1865 % mosoar, mxsoar, n1soar, nosoar,
1866 % moartr, mxartr, n1artr, noartr,
1868 % nbtr, nutr, ierr )
1869 if( ierr .ne. 0 ) return
1872 c 2 milieux => 3 triangles
1873 c -----------------------------------
1874 200 call f2trte( letree(0,nte), pxyd, milieu,
1875 % mosoar, mxsoar, n1soar, nosoar,
1876 % moartr, mxartr, n1artr, noartr,
1878 % nbtr, nutr, ierr )
1879 if( ierr .ne. 0 ) return
1882 c 3 milieux => 4 triangles = 4 quart te
1883 c -------------------------------------
1884 300 call f3trte( letree(0,nte), pxyd, milieu,
1885 % mosoar, mxsoar, n1soar, nosoar,
1886 % moartr, mxartr, n1artr, noartr,
1888 % nbtr, nutr, ierr )
1889 if( ierr .ne. 0 ) return
1898 subroutine aisoar( mosoar, mxsoar, nosoar, na1 )
1899 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1900 c but : chainer en colonne lchain les aretes non vides et
1901 c ----- non frontalieres du tableau nosoar
1905 c mosoar : nombre maximal d'entiers par arete dans le tableau nosoar
1906 c mxsoar : nombre maximal d'aretes frontalieres declarables
1910 c nosoar : numero des 2 sommets , no ligne, 2 triangles, chainages en +
1911 c nosoar(lchain,i)=arete interne suivante
1915 c na1 : numero dans nosoar de la premiere arete interne
1916 c les suivantes sont nosoar(lchain,na1), ...
1917 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1918 c auteur : alain perronnet analyse numerique paris upmc mars 1997
1919 c....................................................................012
1920 parameter (lchain=6)
1921 integer nosoar(mosoar,mxsoar)
1923 c formation du chainage des aretes internes a echanger eventuellement
1924 c recherche de la premiere arete non vide et non frontaliere
1926 if( nosoar(1,na1) .gt. 0 .and. nosoar(3,na1) .le. 0 ) goto 15
1929 c protection de la premiere arete non vide et non frontaliere
1931 do 20 na=na1+1,mxsoar
1932 if( nosoar(1,na) .gt. 0 .and. nosoar(3,na) .le. 0 ) then
1933 c arete interne => elle est chainee a partir de la precedente
1934 nosoar(lchain,na0) = na
1939 c la derniere arete interne n'a pas de suivante
1940 nosoar(lchain,na0) = 0
1944 subroutine tedela( pxyd, noarst,
1945 % mosoar, mxsoar, n1soar, nosoar, n1ardv,
1946 % moartr, mxartr, n1artr, noartr, modifs )
1947 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1948 c but : pour toutes les aretes chainees dans nosoar(lchain,*)
1949 c ----- du tableau nosoar
1950 c echanger la diagonale des 2 triangles si le sommet oppose
1951 c a un triangle ayant en commun une arete appartient au cercle
1952 c circonscrit de l'autre (violation boule vide delaunay)
1956 c pxyd : tableau des x y distance_souhaitee de chaque sommet
1960 c noarst : noarst(i) numero d'une arete de sommet i
1961 c mosoar : nombre maximal d'entiers par arete dans le tableau nosoar
1962 c mxsoar : nombre maximal d'aretes frontalieres declarables
1963 c n1soar : numero de la premiere arete vide dans le tableau nosoar
1964 c nosoar : numero des 2 sommets , no ligne, 2 triangles, chainages en +
1965 c n1ardv : numero dans nosoar de la premiere arete du chainage
1966 c des aretes a rendre delaunay
1968 c moartr : nombre d'entiers par triangle dans le tableau noartr
1969 c mxartr : nombre maximal de triangles declarables dans noartr
1970 c n1artr : numero du premier triangle vide dans le tableau noartr
1971 c le chainage des triangles vides se fait sur noartr(2,.)
1972 c noartr : les 3 aretes des triangles +-arete1, +-arete2, +-arete3
1973 c arete1 = 0 si triangle vide => arete2 = triangle vide suivant
1974 c modifs : nombre d'echanges de diagonales pour maximiser la qualite
1975 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1976 c auteur : alain perronnet analyse numerique paris upmc mars 1997
1977 c....................................................................012
1978 parameter (lchain=6)
1979 common / unites / lecteu, imprim, nunite(30)
1980 double precision pxyd(3,*), surtd2, s123, s142, s143, s234,
1981 % s12, s34, a12, cetria(3), r0
1982 integer nosoar(mosoar,mxsoar),
1983 % noartr(moartr,mxartr),
1986 c le nombre d'echanges de diagonales pour minimiser l'aire
1990 c la premiere arete du chainage des aretes a rendre delaunay
1993 c tant que la pile des aretes a echanger eventuellement est non vide
1994 c ==================================================================
1995 20 if( na0 .gt. 0 ) then
1999 c la prochaine arete a traiter
2000 na0 = nosoar(lchain,na0)
2002 c l'arete est marquee traitee avec le numero -1
2003 nosoar(lchain,na) = -1
2005 c l'arete est elle active?
2006 if( nosoar(1,na) .eq. 0 ) goto 20
2008 c si arete frontaliere pas d'echange possible
2009 if( nosoar(3,na) .gt. 0 ) goto 20
2011 c existe-t-il 2 triangles ayant cette arete commune?
2012 if( nosoar(4,na) .le. 0 .or. nosoar(5,na) .le. 0 ) goto 20
2014 c aucun des 2 triangles est-il desactive?
2015 if( noartr(1,nosoar(4,na)) .eq. 0 .or.
2016 % noartr(1,nosoar(5,na)) .eq. 0 ) goto 20
2018 c l'arete appartient a deux triangles actifs
2019 c le numero des 4 sommets du quadrangle des 2 triangles
2020 call mt4sqa( na, moartr, noartr, mosoar, nosoar,
2021 % ns1, ns2, ns3, ns4 )
2022 if( ns4 .eq. 0 ) goto 20
2024 c carre de la longueur de l'arete ns1 ns2
2025 a12 = (pxyd(1,ns2)-pxyd(1,ns1))**2+(pxyd(2,ns2)-pxyd(2,ns1))**2
2027 c comparaison de la somme des aires des 2 triangles
2028 c -------------------------------------------------
2029 c calcul des surfaces des triangles 123 et 142 de cette arete
2030 s123=surtd2( pxyd(1,ns1), pxyd(1,ns2), pxyd(1,ns3) )
2031 s142=surtd2( pxyd(1,ns1), pxyd(1,ns4), pxyd(1,ns2) )
2032 s12 = abs( s123 ) + abs( s142 )
2033 if( s12 .le. 0.001*a12 ) goto 20
2035 c calcul des surfaces des triangles 143 et 234 de cette arete
2036 s143=surtd2( pxyd(1,ns1), pxyd(1,ns4), pxyd(1,ns3) )
2037 s234=surtd2( pxyd(1,ns2), pxyd(1,ns3), pxyd(1,ns4) )
2038 s34 = abs( s234 ) + abs( s143 )
2040 if( abs(s34-s12) .gt. 1d-15*s34 ) goto 20
2042 c quadrangle convexe : le critere de delaunay intervient
2043 c ------------------ ---------------------------------
2044 c calcul du centre et rayon de la boule circonscrite a 123
2045 c pas d'affichage si le triangle est degenere
2047 call cenced( pxyd(1,ns1), pxyd(1,ns2), pxyd(1,ns3), cetria,
2049 if( ierr .gt. 0 ) then
2050 c ierr=1 si triangle degenere => abandon
2054 if( (cetria(1)-pxyd(1,ns4))**2+(cetria(2)-pxyd(2,ns4))**2
2055 % .lt. cetria(3) ) then
2057 c protection contre une boucle infinie sur le meme cercle
2058 if( r0 .eq. cetria(3) ) goto 20
2060 c oui: ns4 est dans le cercle circonscrit a ns1 ns2 ns3
2061 c => ns3 est aussi dans le cercle circonscrit de ns1 ns2 ns4
2063 cccc les 2 triangles d'arete na sont effaces
2065 ccc nt = nosoar(j,na)
2066 cccc trace du triangle nt
2067 ccc call mttrtr( pxyd, nt, moartr, noartr, mosoar, nosoar,
2068 ccc % ncnoir, ncjaun )
2071 c echange de la diagonale 12 par 34 des 2 triangles
2072 call te2t2t( na, mosoar, n1soar, nosoar, noarst,
2073 % moartr, noartr, na34 )
2074 if( na34 .eq. 0 ) goto 20
2077 c l'arete na34 est marquee traitee
2078 nosoar(lchain,na34) = -1
2081 c les aretes internes peripheriques des 2 triangles sont enchainees
2084 cccc trace du triangle nt
2085 ccc call mttrtr( pxyd, nt, moartr, noartr, mosoar, nosoar,
2086 ccc % ncoran, ncgric )
2088 n = abs( noartr(i,nt) )
2089 if( n .ne. na34 ) then
2090 if( nosoar(3,n) .eq. 0 .and.
2091 % nosoar(lchain,n) .eq. -1 ) then
2092 c cette arete marquee est chainee pour etre traitee
2093 nosoar(lchain,n) = na0
2102 c retour en haut de la pile des aretes a traiter
2108 subroutine terefr( nbarpi, pxyd,
2109 % mosoar, mxsoar, n1soar, nosoar,
2110 % moartr, n1artr, noartr, noarst,
2111 % mxarcf, n1arcf, noarcf, larmin, notrcf,
2113 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2114 c but : recherche des aretes de la frontiere non dans la triangulation
2115 c ----- triangulation frontale pour les reobtenir
2117 c attention: le chainage lchain de nosoar devient celui des cf
2122 c nbarpi : numero du dernier point interne impose par l'utilisateur
2123 c pxyd : tableau des coordonnees 2d des points
2124 c par point : x y distance_souhaitee
2125 c mosoar : nombre maximal d'entiers par arete et
2126 c indice dans nosoar de l'arete suivante dans le hachage
2127 c mxsoar : nombre maximal d'aretes stockables dans le tableau nosoar
2128 c attention: mxsoar>3*mxsomm obligatoire!
2129 c moartr : nombre maximal d'entiers par arete du tableau noartr
2130 c mxarcf : nombre de variables des tableaux n1arcf, noarcf, larmin, notrcf
2134 c n1soar : no de l'eventuelle premiere arete libre dans le tableau nosoar
2135 c chainage des vides suivant en 3 et precedant en 2 de nosoar
2136 c nosoar : numero des 2 sommets , no ligne, 2 triangles de l'arete,
2137 c chainage des aretes frontalieres, chainage du hachage des aretes
2138 c hachage des aretes = nosoar(1)+nosoar(2)*2
2139 c avec mxsoar>=3*mxsomm
2140 c une arete i de nosoar est vide <=> nosoar(1,i)=0 et
2141 c nosoar(2,arete vide)=l'arete vide qui precede
2142 c nosoar(3,arete vide)=l'arete vide qui suit
2143 c n1artr : numero du premier triangle vide dans le tableau noartr
2144 c le chainage des triangles vides se fait sur noartr(2,.)
2145 c noartr : les 3 aretes des triangles +-arete1, +-arete2, +-arete3
2146 c arete1 = 0 si triangle vide => arete2 = triangle vide suivant
2147 c noarst : noarst(i) numero d'une arete de sommet i
2152 c n1arcf : tableau (0:mxarcf) auxiliaire d'entiers
2153 c noarcf : tableau (3,mxarcf) auxiliaire d'entiers
2154 c larmin : tableau (mxarcf) auxiliaire d'entiers
2155 c notrcf : tableau (mxarcf) auxiliaire d'entiers
2159 c nbarpe : nombre d'aretes perdues puis retrouvees
2160 c ierr : =0 si pas d'erreur
2161 c >0 si une erreur est survenue
2162 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2163 c auteur : alain perronnet analyse numerique paris upmc mars 1997
2164 c....................................................................012
2165 parameter (lchain=6)
2166 common / unites / lecteu,imprim,intera,nunite(29)
2167 double precision pxyd(3,*)
2168 integer nosoar(mosoar,mxsoar),
2176 c le nombre d'aretes de la frontiere non arete de la triangulation
2179 c initialisation du chainage des aretes des cf => 0 arete de cf
2180 do 10 narete=1,mxsoar
2181 nosoar( lchain, narete) = -1
2184 c boucle sur l'ensemble des aretes actuelles
2185 c ==========================================
2186 do 30 narete=1,mxsoar
2188 if( nosoar(3,narete) .gt. 0 ) then
2189 c arete appartenant a une ligne => frontaliere
2191 if(nosoar(4,narete) .le. 0 .or. nosoar(5,narete) .le. 0)then
2192 c l'arete narete frontaliere n'appartient pas a 2 triangles
2193 c => elle est perdue
2196 c le numero des 2 sommets de l'arete frontaliere perdue
2197 ns1 = nosoar( 1, narete )
2198 ns2 = nosoar( 2, narete )
2199 c write(imprim,10000) ns1,(pxyd(j,ns1),j=1,2),
2200 c % ns2,(pxyd(j,ns2),j=1,2)
2201 10000 format(' arete perdue a forcer',
2202 % (t24,'sommet=',i6,' x=',g13.5,' y=',g13.5))
2204 c traitement de cette arete perdue ns1-ns2
2205 call tefoar( narete, nbarpi, pxyd,
2206 % mosoar, mxsoar, n1soar, nosoar,
2207 % moartr, n1artr, noartr, noarst,
2208 % mxarcf, n1arcf, noarcf, larmin, notrcf,
2210 if( ierr .ne. 0 ) return
2212 c fin du traitement de cette arete perdue et retrouvee
2220 subroutine tesuex( nblftr, nulftr,
2221 % ndtri0, nbsomm, pxyd, nslign,
2222 % mosoar, mxsoar, nosoar,
2223 % moartr, mxartr, n1artr, noartr, noarst,
2224 % nbtria, letrsu, ierr )
2225 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2226 c but : supprimer du tableau noartr les triangles externes au domaine
2227 c ----- en annulant le numero de leur 1-ere arete dans noartr
2228 c et en les chainant comme triangles vides
2232 c nblftr : nombre de lignes fermees definissant la surface
2233 c nulftr : numero des lignes fermees definissant la surface
2234 c ndtri0 : plus grand numero dans noartr d'un triangle
2235 c pxyd : tableau des coordonnees 2d des points
2236 c par point : x y distance_souhaitee
2237 c nslign : tableau du numero de sommet dans sa ligne pour chaque
2239 c numero du point dans le lexique point si interne impose
2240 c 0 si le point est interne non impose par l'utilisateur
2241 c -1 si le sommet est externe au domaine
2242 c mosoar : nombre maximal d'entiers par arete et
2243 c indice dans nosoar de l'arete suivante dans le hachage
2244 c mxsoar : nombre maximal d'aretes stockables dans le tableau nosoar
2245 c attention: mxsoar>3*mxsomm obligatoire!
2246 c nosoar : numero des 2 sommets , no ligne, 2 triangles de l'arete,
2247 c chainage des aretes frontalieres, chainage du hachage des aretes
2248 c hachage des aretes = nosoar(1)+nosoar(2)*2
2249 c avec mxsoar>=3*mxsomm
2250 c une arete i de nosoar est vide <=> nosoar(1,i)=0 et
2251 c nosoar(2,arete vide)=l'arete vide qui precede
2252 c nosoar(3,arete vide)=l'arete vide qui suit
2253 c moartr : nombre maximal d'entiers par arete du tableau noartr
2254 c mxartr : nombre maximal de triangles declarables
2255 c n1artr : numero du premier triangle vide dans le tableau noartr
2256 c le chainage des triangles vides se fait sur noartr(2,.)
2257 c noartr : les 3 aretes des triangles +-arete1, +-arete2, +-arete3
2258 c arete1 = 0 si triangle vide => arete2 = triangle vide suivant
2259 c noarst : noarst(i) numero nosoar d'une arete de sommet i
2263 c nbtria : nombre de triangles internes au domaine
2264 c letrsu : letrsu(nt)=numero du triangle interne, 0 sinon
2265 c noarst : noarst(i) numero nosoar d'une arete du sommet i (modifi'e)
2266 c ierr : 0 si pas d'erreur, >0 sinon
2267 cc++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2268 c auteur : alain perronnet analyse numerique paris upmc mai 1999
2269 c2345x7..............................................................012
2270 double precision pxyd(3,*)
2271 integer nulftr(nblftr),nslign(nbsomm),
2272 % nosoar(mosoar,mxsoar),
2273 % noartr(moartr,mxartr),
2275 integer letrsu(1:ndtri0)
2276 double precision dmin
2278 c les triangles sont a priori non marques
2283 c les aretes sont marquees non chainees
2284 do 10 noar1=1,mxsoar
2285 nosoar(6,noar1) = -2
2288 c recherche du sommet de la triangulation de plus petite abscisse
2289 c ===============================================================
2293 if( pxyd(1,i) .lt. dmin ) then
2294 c le nouveau minimum
2296 if( noar1 .gt. 0 ) then
2297 c le sommet appartient a une arete de triangle
2298 if( nosoar(4,noar1) .gt. 0 ) then
2299 c le nouveau minimum
2307 c une arete de sommet ntmin
2308 noar1 = noarst( ntmin )
2309 c un triangle d'arete noar1
2310 ntmin = nosoar( 4, noar1 )
2311 if( ntmin .le. 0 ) then
2313 c kerr(1) = 'pas de triangle d''abscisse minimale'
2315 write(imprim,*) 'pas de triangle d''abscisse minimale'
2320 c chainage des 3 aretes du triangle ntmin
2321 c =======================================
2322 c la premiere arete du chainage des aretes traitees
2323 noar1 = abs( noartr(1,ntmin) )
2324 na0 = abs( noartr(2,ntmin) )
2325 c elle est chainee sur la seconde arete du triangle ntmin
2326 nosoar(6,noar1) = na0
2327 c les 2 autres aretes du triangle ntmin sont chainees
2328 na1 = abs( noartr(3,ntmin) )
2329 c la seconde est chainee sur la troisieme arete
2331 c la troisieme n'a pas de suivante
2334 c le triangle ntmin est a l'exterieur du domaine
2335 c tous les triangles externes sont marques -123 456 789
2336 c les triangles de l'autre cote d'une arete sur une ligne
2337 c sont marques: no de la ligne de l'arete * signe oppose
2338 c =======================================================
2340 ligne = -123 456 789
2342 40 if( noar1 .ne. 0 ) then
2344 c l'arete noar1 du tableau nosoar est a traiter
2345 c ---------------------------------------------
2347 c l'arete suivante devient la premiere a traiter ensuite
2348 noar1 = nosoar(6,noar1)
2349 c l'arete noar est traitee
2354 c l'un des 2 triangles de l'arete
2356 if( nt .gt. 0 ) then
2358 c triangle deja traite pour une ligne anterieure?
2359 if( letrsu(nt) .ne. 0 .and.
2360 % abs(letrsu(nt)) .ne. ligne ) goto 60
2362 cccc trace du triangle nt en couleur ligne0
2363 ccc call mttrtr( pxyd, nt, moartr, noartr, mosoar, nosoar,
2364 ccc % ligne0, ncnoir )
2366 c le triangle est marque avec la valeur de ligne
2369 c chainage eventuel des autres aretes de ce triangle
2370 c si ce n'est pas encore fait
2373 c le numero na de l'arete j du triangle nt dans nosoar
2374 na = abs( noartr(j,nt) )
2375 if( nosoar(6,na) .ne. -2 ) goto 50
2377 c le numero de 1 a nblftr dans nulftr de la ligne de l'arete
2380 c si l'arete est sur une ligne fermee differente de celle envelo
2381 c et non marquee alors examen du triangle oppose
2382 if( nl .gt. 0 ) then
2384 if( nl .eq. ligne0 ) goto 50
2386 c arete frontaliere de ligne non traitee
2387 c => passage de l'autre cote de la ligne
2388 c le triangle de l'autre cote de la ligne est recherche
2389 if( nt .eq. abs( nosoar(4,na) ) ) then
2394 nt2 = abs( nosoar(nt2,na) )
2395 if( nt2 .gt. 0 ) then
2397 c le triangle nt2 de l'autre cote est marque avec le
2398 c avec le signe oppose de celui de ligne
2399 if( ligne .ge. 0 ) then
2404 letrsu(nt2) = lsigne * nl
2406 c temoin de ligne a traiter ensuite dans nulftr
2407 nulftr(nl) = -abs( nulftr(nl) )
2409 cccc trace du triangle nt2 en jaune borde de magenta
2410 ccc call mttrtr( pxyd,nt2,
2411 ccc % moartr,noartr,mosoar,nosoar,
2412 ccc % ncjaun, ncmage )
2414 c l'arete est traitee
2419 c l'arete est traitee
2424 c arete non traitee => elle est chainee
2425 nosoar(6,na) = noar1
2435 c les triangles de la ligne fermee ont tous ete marques
2436 c plus d'arete chainee
2438 c recherche d'une nouvelle ligne fermee a traiter
2439 c ===============================================
2440 65 do 70 nl=1,nblftr
2441 if( nulftr(nl) .lt. 0 ) goto 80
2443 c plus de ligne fermee a traiter
2446 c tous les triangles de cette composante connexe
2447 c entre ligne et ligne0 vont etre marques
2448 c ==============================================
2449 c remise en etat du numero de ligne
2450 c nl est le numero de la ligne dans nulftr a traiter
2451 80 nulftr(nl) = -nulftr(nl)
2453 if( abs(letrsu(nt2)) .eq. nl ) goto 92
2456 c recherche de l'arete j du triangle nt2 avec ce numero de ligne nl
2459 c le numero de l'arete j du triangle dans nosoar
2461 na0 = abs( noartr(j,nt2) )
2462 if( nl .eq. nosoar(3,na0) ) then
2464 c na0 est l'arete de ligne nl
2465 c l'arete suivante du triangle nt2
2467 c le numero dans nosoar de l'arete i de nt2
2468 na1 = abs( noartr(i,nt2) )
2469 if( nosoar(6,na1) .eq. -2 ) then
2470 c arete non traitee => elle est la premiere du chainage
2472 c pas de suivante dans ce chainage
2478 c l'eventuelle seconde arete suivante
2480 na = abs( noartr(i,nt2) )
2481 if( nosoar(6,na) .eq. -2 ) then
2482 if( na1 .eq. 0 ) then
2483 c 1 arete non traitee et seule a chainer
2487 c 2 aretes a chainer
2493 if( noar1 .gt. 0 ) then
2495 c il existe au moins une arete a visiter pour ligne
2496 c marquage des triangles internes a la ligne nl
2503 c nt2 est le seul triangle de la ligne fermee
2510 c reperage des sommets internes ou externes dans nslign
2511 c nslign(sommet externe au domaine)=-1
2512 c nslign(sommet interne au domaine)= 0
2513 c =====================================================
2514 110 do 170 ns1=1,nbsomm
2515 c tout sommet non sur la frontiere ou interne impose
2516 c est suppose externe
2517 if( nslign(ns1) .eq. 0 ) nslign(ns1) = -1
2520 c les triangles externes sont marques vides dans le tableau noartr
2521 c ================================================================
2525 if( letrsu(nt) .le. 0 ) then
2527 c triangle nt externe
2528 if( noartr(1,nt) .ne. 0 ) then
2529 c la premiere arete est annulee
2531 c le triangle nt est considere comme etant vide
2532 noartr(2,nt) = n1artr
2538 c triangle nt interne
2542 c marquage des 3 sommets du triangle nt
2544 c le numero nosoar de l'arete i du triangle nt
2545 noar = abs( noartr(i,nt) )
2546 c le numero des 2 sommets
2547 ns1 = nosoar(1,noar)
2548 ns2 = nosoar(2,noar)
2549 c mise a jour du numero d'une arete des 2 sommets de l'arete
2550 noarst( ns1 ) = noar
2551 noarst( ns2 ) = noar
2552 c ns1 et ns2 sont des sommets de la triangulation du domaine
2553 if( nslign(ns1) .lt. 0 ) nslign(ns1)=0
2554 if( nslign(ns2) .lt. 0 ) nslign(ns2)=0
2560 c ici tout sommet externe ns verifie nslign(ns)=-1
2562 c les triangles externes sont mis a zero dans nosoar
2563 c ==================================================
2564 do 300 noar=1,mxsoar
2566 if( nosoar(1,noar) .gt. 0 ) then
2568 c le second triangle de l'arete noar
2570 if( nt .gt. 0 ) then
2571 c si le triangle nt est externe
2572 c alors il est supprime pour l'arete noar
2573 if( letrsu(nt) .le. 0 ) nosoar(5,noar)=0
2576 c le premier triangle de l'arete noar
2578 if( nt .gt. 0 ) then
2579 if( letrsu(nt) .le. 0 ) then
2580 c si le triangle nt est externe
2581 c alors il est supprime pour l'arete noar
2582 c et l'eventuel triangle oppose prend sa place
2583 c en position 4 de nosoar
2584 if( nosoar(5,noar) .gt. 0 ) then
2585 nosoar(4,noar)=nosoar(5,noar)
2596 c remise en etat pour eviter les modifications de ladefi
2597 9990 do 9991 nl=1,nblftr
2598 if( nulftr(nl) .lt. 0 ) nulftr(nl)=-nulftr(nl)
2605 subroutine trp1st( ns, noarst, mosoar, nosoar, moartr, noartr,
2606 % mxpile, lhpile, lapile )
2607 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2608 c but : recherche des triangles de noartr partageant le sommet ns
2610 c limite: un camembert de centre ns entame 2 fois
2611 c ne donne que l'une des parties
2615 c ns : numero du sommet
2616 c noarst : noarst(i) numero d'une arete de sommet i
2617 c mosoar : nombre maximal d'entiers par arete et
2618 c indice dans nosoar de l'arete suivante dans le hachage
2619 c nosoar : numero des 2 sommets , no ligne, 2 triangles de l'arete,
2620 c chainage des aretes frontalieres, chainage du hachage des aretes
2621 c moartr : nombre maximal d'entiers par arete du tableau noartr
2622 c noartr : les 3 aretes des triangles +-arete1, +-arete2, +-arete3
2623 c mxpile : nombre maximal de triangles empilables
2627 c lhpile : >0 nombre de triangles empiles
2628 c =0 si impossible de tourner autour du point
2629 c =-lhpile si apres butee sur la frontiere il y a a nouveau
2630 c butee sur la frontiere . a ce stade on ne peut dire si tous
2631 c les triangles ayant ce sommet ont ete recenses
2632 c ce cas arrive seulement si le sommet est sur la frontiere
2633 c lapile : numero dans noartr des triangles de sommet ns
2634 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2635 c auteur : alain perronnet analyse numerique paris upmc mars 1997
2636 c....................................................................012
2637 common / unites / lecteu, imprim, nunite(30)
2638 integer noartr(moartr,*),
2641 integer lapile(1:mxpile)
2644 c la premiere arete de sommet ns
2646 if( nar .le. 0 ) then
2647 write(imprim,*) 'trp1st: sommet',ns,' sans arete'
2651 c l'arete nar est elle active?
2652 if( nosoar(1,nar) .le. 0 ) then
2653 ccc write(imprim,*) 'trp1st: arete vide',nar,
2654 ccc % ' st1:', nosoar(1,nar),' st2:',nosoar(2,nar)
2658 c le premier triangle de sommet ns
2659 nt0 = abs( nosoar(4,nar) )
2660 if( nt0 .le. 0 ) then
2661 write(imprim,*) 'trp1st: sommet',ns,' dans aucun triangle'
2665 c le triangle est il interne?
2666 if( noartr(1,nt0) .eq. 0 ) goto 9999
2668 c le numero des 3 sommets du triangle nt0 dans le sens direct
2669 call nusotr( nt0, mosoar, nosoar, moartr, noartr, nosotr )
2671 c reperage du sommet ns dans le triangle nt0
2673 if( nosotr(nar) .eq. ns ) goto 10
2678 c ns retrouve : le triangle nt0 est empile
2683 c recherche dans le sens des aiguilles d'une montre
2684 c (sens indirect) du triangle nt1 de l'autre cote de l'arete
2685 c nar du triangle et en tournant autour du sommet ns
2686 c ==========================================================
2687 noar = abs( noartr(nar,nt0) )
2688 c le triangle nt1 oppose du triangle nt0 par l'arete noar
2689 if( nosoar(4,noar) .eq. nt0 ) then
2690 nt1 = nosoar(5,noar)
2692 nt1 = nosoar(4,noar)
2695 c la boucle sur les triangles nt1 de sommet ns dans le sens indirect
2696 c ==================================================================
2697 if( nt1 .gt. 0 ) then
2699 if( noartr(1,nt1) .eq. 0 ) goto 30
2701 c le triangle nt1 n'a pas ete detruit. il est actif
2702 c le triangle oppose par l'arete noar existe
2703 c le numero des 3 sommets du triangle nt1 dans le sens direct
2704 15 call nusotr( nt1, mosoar, nosoar, moartr, noartr, nosotr )
2706 c reperage du sommet ns dans nt1
2708 if( nosotr(nar) .eq. ns ) goto 25
2714 25 if( lhpile .ge. mxpile ) goto 9990
2716 lapile(lhpile) = nt1
2718 c le triangle nt1 de l'autre cote de l'arete de sommet ns
2719 c sauvegarde du precedent triangle dans nta
2721 noar = abs( noartr(nar,nt1) )
2722 if( nosoar(4,noar) .eq. nt1 ) then
2723 nt1 = nosoar(5,noar)
2725 nt1 = nosoar(4,noar)
2727 if( nt1 .le. 0 ) goto 30
2728 c le triangle suivant est a l'exterieur
2729 if( nt1 .ne. nt0 ) goto 15
2731 c recherche terminee par arrivee sur nt0
2732 c les triangles forment un "cercle" de "centre" ns
2737 c pas de triangle voisin a nt1
2738 c ============================
2739 c le parcours passe par 1 des triangles exterieurs
2740 c le parcours est inverse par l'arete de gauche
2741 c le triangle nta est le premier triangle empile
2743 lapile(lhpile) = nta
2745 c le numero des 3 sommets du triangle nta dans le sens direct
2746 call nusotr( nta, mosoar, nosoar, moartr, noartr, nosotr )
2748 if( nosotr(nar) .eq. ns ) goto 33
2752 c l'arete qui precede (rotation / ns dans le sens direct)
2753 33 if( nar .eq. 1 ) then
2759 c le triangle voisin de nta dans le sens direct
2760 noar = abs( noartr(nar,nta) )
2761 if( nosoar(4,noar) .eq. nta ) then
2762 nt1 = nosoar(5,noar)
2764 nt1 = nosoar(4,noar)
2766 if( nt1 .le. 0 ) then
2767 c un seul triangle contient ns
2771 c boucle sur les triangles de sommet ns dans le sens direct
2772 c ==========================================================
2773 40 if( noartr(1,nt1) .eq. 0 ) goto 70
2775 c le triangle nt1 n'a pas ete detruit. il est actif
2776 c le numero des 3 sommets du triangle nt1 dans le sens direct
2777 call nusotr( nt1, mosoar, nosoar, moartr, noartr, nosotr )
2779 c reperage du sommet ns dans nt1
2781 if( nosotr(nar) .eq. ns ) goto 60
2787 60 if( lhpile .ge. mxpile ) goto 9990
2789 lapile(lhpile) = nt1
2791 c l'arete qui precede dans le sens direct
2792 if( nar .eq. 1 ) then
2798 c l'arete de sommet ns dans nosoar
2799 noar = abs( noartr(nar,nt1) )
2801 c le triangle voisin de nta dans le sens direct
2803 if( nosoar(4,noar) .eq. nt1 ) then
2804 nt1 = nosoar(5,noar)
2806 nt1 = nosoar(4,noar)
2809 if( nt1 .gt. 0 ) goto 40
2811 c butee sur le trou => fin des triangles de sommet ns
2812 c ----------------------------------------------------
2814 c impossible ici de trouver les autres triangles de sommet ns
2815 c les triangles de sommet ns ne forment pas une boule de centre ns
2818 c saturation de la pile des triangles
2819 c -----------------------------------
2820 9990 write(imprim,*) 'trp1st:saturation pile des triangles autour ',
2824 c erreur triangle ne contenant pas le sommet ns
2825 c ----------------------------------------------
2826 9995 write(imprim,*) 'trp1st:triangle ',nta,' st=',
2827 % (nosotr(nar),nar=1,3),' sans le sommet' ,ns
2835 subroutine nusotr( nt, mosoar, nosoar, moartr, noartr, nosotr )
2836 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2837 c but : calcul du numero des 3 sommets du triangle nt de noartr
2838 c ----- dans le sens direct (aire>0 si non degenere)
2842 c nt : numero du triangle dans le tableau noartr
2843 c mosoar : nombre maximal d'entiers par arete
2844 c nosoar : numero des 2 sommets , no ligne, 2 triangles, chainages en +
2845 c sommet 1 = 0 si arete vide => sommet 2 = arete vide suivante
2846 c moartr : nombre maximal d'entiers par arete du tableau noartr
2847 c noartr : les 3 aretes des triangles +-arete1, +-arete2, +-arete3
2848 c arete1=0 si triangle vide => arete2=triangle vide suivant
2852 c nosotr : numero (dans le tableau pxyd) des 3 sommets du triangle
2853 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2854 c auteur : alain perronnet analyse numerique paris upmc mars 1997
2855 c2345x7..............................................................012
2856 integer nosoar(mosoar,*), noartr(moartr,*), nosotr(3)
2858 c les 2 sommets de l'arete 1 du triangle nt dans le sens direct
2859 na = noartr( 1, nt )
2860 if( na .gt. 0 ) then
2868 nosotr(1) = nosoar( nosotr(1), na )
2869 nosotr(2) = nosoar( nosotr(2), na )
2872 na = abs( noartr(2,nt) )
2874 c le sommet nosotr(3 du triangle 123
2875 nosotr(3) = nosoar( 1, na )
2876 if( nosotr(3) .eq. nosotr(1) .or. nosotr(3) .eq. nosotr(2) ) then
2877 nosotr(3) = nosoar(2,na)
2882 subroutine tesusp( nbarpi, pxyd, noarst,
2883 % mosoar, mxsoar, n1soar, nosoar,
2884 % moartr, mxartr, n1artr, noartr,
2885 % mxarcf, n1arcf, noarcf, larmin, notrcf, liarcf,
2887 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2888 c but : supprimer de la triangulation les sommets de te trop proches
2889 c ----- soit d'un sommet frontalier ou point interne impose
2890 c soit d'une arete frontaliere
2892 c attention: le chainage lchain de nosoar devient celui des cf
2896 c nbarpi : numero du dernier point interne impose par l'utilisateur
2897 c pxyd : tableau des coordonnees 2d des points
2898 c par point : x y distance_souhaitee
2899 c mosoar : nombre maximal d'entiers par arete et
2900 c indice dans nosoar de l'arete suivante dans le hachage
2901 c mxsoar : nombre maximal d'aretes stockables dans le tableau nosoar
2902 c attention: mxsoar>3*mxsomm obligatoire!
2903 c moartr : nombre maximal d'entiers par arete du tableau noartr
2904 c mxarcf : nombre de variables des tableaux n1arcf, noarcf, larmin, notrcf
2908 c noarst : noarst(i) numero d'une arete de sommet i
2909 c n1soar : no de l'eventuelle premiere arete libre dans le tableau nosoar
2910 c chainage des vides suivant en 3 et precedant en 2 de nosoar
2911 c nosoar : numero des 2 sommets , no ligne, 2 triangles de l'arete,
2912 c chainage des aretes frontalieres, chainage du hachage des aretes
2913 c hachage des aretes = nosoar(1)+nosoar(2)*2
2914 c avec mxsoar>=3*mxsomm
2915 c une arete i de nosoar est vide <=> nosoar(1,i)=0 et
2916 c nosoar(2,arete vide)=l'arete vide qui precede
2917 c nosoar(3,arete vide)=l'arete vide qui suit
2918 c n1artr : numero du premier triangle vide dans le tableau noartr
2919 c le chainage des triangles vides se fait sur noartr(2,.)
2920 c noartr : les 3 aretes des triangles +-arete1, +-arete2, +-arete3
2921 c arete1 = 0 si triangle vide => arete2 = triangle vide suivant
2926 c n1arcf : tableau (0:mxarcf) auxiliaire d'entiers
2927 c noarcf : tableau (3,mxarcf) auxiliaire d'entiers
2928 c larmin : tableau ( mxarcf ) auxiliaire d'entiers
2929 c notrcf : tableau ( mxarcf ) auxiliaire d'entiers
2930 c liarcf : tableau ( mxarcf ) auxiliaire d'entiers
2934 c nbstsu : nombre de sommets de te supprimes
2935 c ierr : =0 si pas d'erreur
2936 c >0 si une erreur est survenue
2937 c 11 algorithme defaillant
2938 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2939 c auteur : alain perronnet analyse numerique paris upmc mars 1997
2940 c....................................................................012
2941 c parameter ( quamal=0.3 ) => ok
2942 c parameter ( quamal=0.4 ) => pb pour le test ocean
2943 c parameter ( quamal=0.5 ) => pb pour le test ocean
2945 parameter ( quamal=0.333, lchain=6 )
2946 common / unites / lecteu,imprim,intera,nunite(29)
2947 double precision pxyd(3,*), qualit
2948 integer nosoar(mosoar,mxsoar),
2958 equivalence (nosotr(1),ns1), (nosotr(2),ns2),
2961 c le nombre de sommets de te supprimes
2964 c initialisation du chainage des aretes des cf => 0 arete de cf
2965 do 10 narete=1,mxsoar
2966 nosoar( lchain, narete ) = -1
2969 c boucle sur l'ensemble des sommets frontaliers ou points internes
2970 c ================================================================
2971 do 100 ns = 1, nbarpi
2973 cccc le nombre de sommets supprimes pour ce sommet ns
2976 c la qualite minimale au dessous de laquelle le point proche
2977 c interne est supprime
2980 c une arete de sommet ns
2981 15 narete = noarst( ns )
2982 if( narete .le. 0 ) then
2983 c erreur: le point appartient a aucune arete
2984 write(imprim,*) 'sommet ',ns,' dans aucune arete'
2989 c recherche des triangles de sommet ns
2990 c ils doivent former un contour ferme de type etoile
2991 call trp1st( ns, noarst, mosoar, nosoar, moartr, noartr,
2992 % mxarcf, nbtrcf, notrcf )
2993 if( nbtrcf .le. 0 ) then
2994 c erreur: impossible de trouver tous les triangles de sommet ns
2995 c seule une partie est a priori retrouvee
2999 c boucle sur les triangles de l'etoile du sommet ns
3003 c le numero des 3 sommets du triangle nt
3005 call nusotr( nt, mosoar, nosoar, moartr, noartr,
3007 c nosotr(1:3) est en equivalence avec ns1, ns2, ns3
3009 c la qualite du triangle ns1 ns2 ns3
3010 call qutr2d( pxyd(1,ns1), pxyd(1,ns2), pxyd(1,ns3), qualit )
3011 if( qualit .lt. quamin ) then
3017 c bilan sur la qualite des triangles de sommet ns
3018 if( quamin .lt. quaopt ) then
3020 c recherche du sommet de ntqmin le plus proche et non frontalier
3021 c ==============================================================
3022 c le numero des 3 sommets du triangle nt
3023 call nusotr( ntqmin, mosoar, nosoar, moartr, noartr,
3028 if( nosotr(j) .ne. ns .and. nosotr(j) .gt. nbarpi ) then
3029 d = (pxyd(1,nosotr(j))-pxyd(1,ns))**2
3030 % + (pxyd(2,nosotr(j))-pxyd(2,ns))**2
3031 if( d .lt. quamin ) then