src/MEFISTO2/trte.f

   1 c  MEFISTO : library to compute 2D triangulation from segmented boundaries
   2 c
   3 c  Copyright (C) 2003  Laboratoire J.-L. Lions UPMC Paris
   4 c
   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.
   9 c
  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.
  14 c
  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
  18 c
  19 c  See http://www.ann.jussieu.fr/~perronne or email Perronnet@ann.jussieu.fr
  20 c
  21 c
  22 c  File   : trte.f
  23 c  Module : SMESH
  24 c  Author: Alain PERRONNET
  25
  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
  30 c
  31 c entrees :
  32 c ---------
  33 c p1,p2,p3 : les 3 coordonnees des 3 sommets du triangle
  34 c            sens direct pour une surface et qualite >0
  35 c sorties :
  36 c ---------
  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
  45 c
  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 )
  50 c
  51 c     demi perimetre
  52       p = (a+b+c) * 0.5d0
  53 c
  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) ) )
  57      %          / max(a,b,c)
  58       else
  59          qualite = 0d0
  60       endif
  61 c
  62 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)
  66 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)
  69 c
  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)
  72       end
  73
  74
  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
  78 c -----
  79 c parametres d entree :
  80 c ---------------------
  81 c p1 p2 p3 : les 3 fois 2 coordonnees des sommets du triangle
  82 c
  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)
  90 c
  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
  94       end
  95
  96       integer function nopre3( i )
  97 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  98 c but :   numero precedent i dans le sens circulaire  1 2 3 1 ...
  99 c -----
 100 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 101 c auteur : alain perronnet  analyse numerique paris upmc    fevrier 1992
 102 c2345x7..............................................................012
 103       if( i .eq. 1 ) then
 104          nopre3 = 3
 105       else
 106          nopre3 = i - 1
 107       endif
 108       end
 109
 110       integer function nosui3( i )
 111 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 112 c but :   numero suivant i dans le sens circulaire  1 2 3 1 ...
 113 c -----
 114 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 115 c auteur : alain perronnet  analyse numerique paris upmc    fevrier 1992
 116 c2345x7..............................................................012
 117       if( i .eq. 3 ) then
 118          nosui3 = 1
 119       else
 120          nosui3 = i + 1
 121       endif
 122       end
 123
 124       subroutine provec( v1 , v2 , v3 )
 125 c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 126 c but :    v3 vecteur = produit vectoriel de 2 vecteurs de r ** 3
 127 c -----
 128 c entrees:
 129 c --------
 130 c v1, v2 : les 2 vecteurs de 3 composantes
 131 c
 132 c sortie :
 133 c --------
 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)
 139 c
 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 )
 143 c
 144       return
 145       end
 146
 147       subroutine norme1( n, v, ierr )
 148 c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 149 c but :   normalisation euclidienne a 1 d un vecteur v de n composantes
 150 c -----
 151 c entrees :
 152 c ---------
 153 c n       : nombre de composantes du vecteur
 154 c
 155 c modifie :
 156 c ---------
 157 c v       : le vecteur a normaliser a 1
 158 c
 159 c sortie  :
 160 c ---------
 161 c ierr    : 1 si la norme de v est egale a 0
 162 c           0 si pas d'erreur
 163 c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 164 c auteur : alain perronnet analyse numerique paris             mars 1987
 165 c ......................................................................
 166       double precision  v( n ), s, sqrt
 167 c
 168       s = 0.0d0
 169       do 10 i=1,n
 170          s = s + v( i ) * v( i )
 171    10 continue
 172 c
 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
 177          ierr = 1
 178          return
 179       endif
 180 c
 181       s = 1.0d0 / sqrt( s )
 182       do 20 i=1,n
 183          v( i ) = v ( i ) * s
 184    20 continue
 185 c
 186       ierr = 0
 187       end
 188
 189
 190       subroutine insoar( mxsomm, mosoar, mxsoar, n1soar, nosoar )
 191 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 192 c but :    initialiser le tableau nosoar pour le hachage des aretes
 193 c -----
 194 c
 195 c entrees:
 196 c --------
 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
 201 c
 202 c sorties:
 203 c --------
 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)
 216 c
 217 c     initialisation des aretes 1 a mxsomm
 218       do 10 i=1,mxsomm
 219 c
 220 c        sommet 1 = 0 <=> temoin d'arete vide pour le hachage
 221          nosoar( 1, i ) = 0
 222 c
 223 c        arete sur aucune ligne
 224          nosoar( 3, i ) = 0
 225 c
 226 c        la position de l'arete interne ou frontaliere est inconnue
 227          nosoar( 6, i ) = -2
 228 c
 229 c        fin de chainage du hachage pas d'arete suivante
 230          nosoar( mosoar, i ) = 0
 231 c
 232  10   continue
 233 c
 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
 236       n1soar = mxsomm + 1
 237 c
 238 c     initialisation des aretes vides et des chainages
 239       do 20 i = n1soar, mxsoar
 240 c
 241 c        sommet 1 = 0 <=> temoin d'arete vide pour le hachage
 242          nosoar( 1, i ) = 0
 243 c
 244 c        arete sur aucune ligne
 245          nosoar( 3, i ) = 0
 246 c
 247 c        chainage sur l'arete vide qui precede
 248 c        (si arete occupee cela deviendra le no du triangle 1 de l'arete)
 249          nosoar( 4, i ) = i-1
 250 c
 251 c        chainage sur l'arete vide qui suit
 252 c        (si arete occupee cela deviendra le no du triangle 2 de l'arete)
 253          nosoar( 5, i ) = i+1
 254 c
 255 c        chainages des aretes frontalieres ou internes ou ...
 256          nosoar( 6, i ) = -2
 257 c
 258 c        fin de chainage du hachage
 259          nosoar( mosoar, i ) = 0
 260 c
 261  20   continue
 262 c
 263 c     la premiere arete vide n'a pas de precedent
 264       nosoar( 4, n1soar ) = 0
 265 c
 266 c     la derniere arete vide est mxsoar sans arete vide suivante
 267       nosoar( 5, mxsoar ) = 0
 268       end
 269
 270
 271       subroutine azeroi ( l , ntab )
 272 c ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 273 c but : initialisation a zero d un tableau ntab de l variables entieres
 274 c -----
 275 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 276 c auteur : alain perronnet analyse numerique upmc paris septembre 1988
 277 c23456---------------------------------------------------------------012
 278       integer ntab(l)
 279       do 1 i = 1 , l
 280          ntab( i ) = 0
 281     1 continue
 282       end
 283
 284
 285       subroutine fasoar( ns1,    ns2,    nt1,    nt2,    nolign,
 286      %                   mosoar, mxsoar, n1soar, nosoar, noarst,
 287      %                   noar,   ierr )
 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
 291 c
 292 c entrees:
 293 c --------
 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
 304 c
 305 c modifies:
 306 c ---------
 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
 316 c
 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
 320 c
 321 c sorties:
 322 c --------
 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(*)
 337       integer           nu2sar(2)
 338 c
 339 c     ajout eventuel de l'arete s1 s2 dans nosoar
 340       nu2sar(1) = ns1
 341       nu2sar(2) = ns2
 342 c
 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
 348 c
 349       if( noar .eq. 0 ) then
 350 c
 351 c        saturation du tableau nosoar
 352          write(imprim,*) 'fasoar: tableau nosoar sature'
 353          ierr = 1
 354          return
 355 c
 356       else if( noar .lt. 0 ) then
 357 c
 358 c        l'arete a ete ajoutee. initialisation des autres informations
 359          noar = -noar
 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
 363          nosoar(4,noar) = nt1
 364 c        le triangle 2 de l'arete => le triangle nt2
 365          nosoar(5,noar) = nt2
 366 c
 367 c        le sommet appartient a l'arete noar
 368          noarst( nu2sar(1) ) = noar
 369          noarst( nu2sar(2) ) = noar
 370 c
 371       else
 372 c
 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!'
 386                  endif
 387                  ierr = 2
 388                  return
 389              endif
 390          endif
 391 c
 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
 396              n = 4
 397          else
 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'
 404                 ierr = 3
 405                 return
 406              endif
 407              n = 5
 408          endif
 409          nosoar(n,noar) = nt1
 410 c
 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'
 419                 ierr = 4
 420                 return
 421             endif
 422             nosoar(5,noar) = nt2
 423          endif
 424 c
 425       endif
 426 c
 427 c     pas d'erreur
 428       ierr = 0
 429       end
 430
 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
 436 c
 437 c entrees:
 438 c --------
 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
 441 c
 442 c sorties:
 443 c --------
 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
 452 c
 453       a = s(1,1)
 454       b = s(1,2) - s(1,1)
 455       c = s(1,4) - s(1,1)
 456       d = s(1,1) - s(1,2) + s(1,3) - s(1,4)
 457 c
 458       alpha = s(2,1)
 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)
 462 c
 463       u = beta  * c - b * gamma
 464       if( u .eq. 0 ) then
 465 c        quadrangle degenere
 466          ierr = 1
 467          return
 468       endif
 469       v = delta * c - d * gamma
 470       w = b * delta - beta * d
 471 c
 472       x0 = c * (y-alpha) - gamma * (x-a)
 473       y0 = b * (y-alpha) - beta  * (x-a)
 474 c
 475       a = v  * w
 476       b = u  * u - w * x0 - v * y0
 477       c = x0 * y0
 478 c
 479       if( a .ne. 0 ) then
 480 c
 481          delta = sqrt( b*b-4*a*c )
 482          if( b .ge. 0.0 ) then
 483             t(2) = -b - delta
 484          else
 485             t(2) = -b + delta
 486          endif
 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
 492          t(1) = - b/a - t(2)
 493 c
 494          do 10 i=1,2
 495 c
 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
 501             endif
 502 c
 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 )
 507 c
 508  10      continue
 509 c
 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
 513             goto 9000
 514          endif
 515 c
 516       else if ( b .ne. 0 ) then
 517          t(1) = - c / b
 518       else
 519          t(1) = 0
 520       endif
 521 c
 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
 525 c
 526  9000 ierr = 0
 527       return
 528       end
 529
 530
 531       subroutine ptdatr( point, pxyd, nosotr, nsigne )
 532 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 533 c but :    le point est il dans le triangle de sommets nosotr
 534 c -----
 535 c
 536 c entrees:
 537 c --------
 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
 541 c
 542 c sorties:
 543 c --------
 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
 549       integer           nosotr(3)
 550       double precision  point(2), pxyd(3,*)
 551       double precision  xp,yp, x1,x2,x3, y1,y2,y3, d,dd, cb1,cb2,cb3
 552 c
 553       xp = point( 1 )
 554       yp = point( 2 )
 555 c
 556       n1 = nosotr( 1 )
 557       x1 = pxyd( 1 , n1 )
 558       y1 = pxyd( 2 , n1 )
 559 c
 560       n2 = nosotr( 2 )
 561       x2 = pxyd( 1 , n2 )
 562       y2 = pxyd( 2 , n2 )
 563 c
 564       n3 = nosotr( 3 )
 565       x3 = pxyd( 1 , n3 )
 566       y3 = pxyd( 2 , n3 )
 567 c
 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 )
 571 c
 572       if( d .gt. 0 ) then
 573 c
 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
 580          cb3 = 1d0 - cb1 -cb2
 581 ccc         cb3 = ( ( x1-xp ) * ( y2-yp ) - ( x2-xp ) * ( y1-yp ) ) / d
 582 c
 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
 587 c
 588 c           le triangle nosotr contient le point
 589             nsigne = 1
 590          else
 591             nsigne = 0
 592          endif
 593 c
 594       else
 595 c
 596 c        triangle degenere
 597 c        =================
 598 c        le point est il du meme cote que le sommet oppose de chaque arete?
 599          nsigne = 0
 600          do 10 i=1,3
 601 c           le sinus de l'angle p1 p2-p1 point
 602             x1  = pxyd(1,n1)
 603             y1  = pxyd(2,n1)
 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
 618 c                 point sur l'arete
 619                   nsigne = nsigne + 1
 620                endif
 621             else
 622 c              le point 3 n'est pas sur l'arete . test des signes
 623                if( d * dd .ge. 0 ) then
 624                   nsigne = nsigne + 1
 625                endif
 626             endif
 627 c           permutation circulaire des 3 sommets et aretes
 628             n  = n1
 629             n1 = n2
 630             n2 = n3
 631             n3 = n
 632  10      continue
 633          if( nsigne .ne. 3 ) nsigne = 0
 634       endif
 635       end
 636
 637       integer function nosstr( p, pxyd, nt, letree )
 638 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 639 c but :    calculer le numero 0 a 3 du sous-triangle te contenant
 640 c -----    le point p
 641 c
 642 c entrees:
 643 c --------
 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
 654 c      sinon
 655 c         letree(0:3,j) :-no pxyd des 1 \85a 4 points internes au triangle j
 656 c                         0  si pas de point
 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
 661 c
 662 c sorties :
 663 c ---------
 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
 672 c
 673 c     le numero des 3 sommets du triangle
 674       ns1 = letree( 6, nt )
 675       ns2 = letree( 7, nt )
 676       ns3 = letree( 8, nt )
 677 c
 678 c     les coordonnees entre 0 et 1 du point p
 679       x1  = pxyd(1,ns1)
 680       y1  = pxyd(2,ns1)
 681 c
 682       x21 = pxyd(1,ns2) - x1
 683       y21 = pxyd(2,ns2) - y1
 684 c
 685       x31 = pxyd(1,ns3) - x1
 686       y31 = pxyd(2,ns3) - y1
 687 c
 688       d   = 1.0 / ( x21 * y31 - x31 * y21 )
 689 c
 690       xe  = ( ( p(1) - x1 ) * y31 - ( p(2) - y1 ) * x31 ) * d
 691       ye  = ( ( p(2) - y1 ) * x21 - ( p(1) - x1 ) * y21 ) * d
 692 c
 693       if( xe .gt. 0.5d0 ) then
 694 c        sous-triangle droit
 695          nosstr = 2
 696       else if( ye .gt. 0.5d0 ) then
 697 c        sous-triangle haut
 698          nosstr = 3
 699       else if( xe+ye .lt. 0.5d0 ) then
 700 c        sous-triangle gauche
 701          nosstr = 1
 702       else
 703 c        sous-triangle central
 704          nosstr = 0
 705       endif
 706       end
 707
 708
 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
 713 c
 714 c entrees:
 715 c --------
 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
 726 c      sinon
 727 c         letree(0:3,j) :-no pxyd des 1 \85 4 points internes au triangle j
 728 c                         0  si pas de point
 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
 733 c
 734 c sorties :
 735 c ---------
 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)
 742 c
 743 c     la racine depart de la recherche
 744       notrpt = notrde
 745 c
 746 c     tant que la feuille n'est pas atteinte descendre l'arbre
 747  10   if( letree(0,notrpt) .gt. 0 ) then
 748 c
 749 c        recherche du sous-triangle contenant p
 750          nsot = nosstr( p, pxyd, notrpt, letree )
 751 c
 752 c        le numero letree du sous-triangle
 753          notrpt = letree( nsot, notrpt )
 754          goto 10
 755 c
 756       endif
 757       end
 758
 759
 760       subroutine teajpt( ns,   nbsomm, mxsomm, pxyd, letree,
 761      &                   ntrp, ierr )
 762 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 763 c but :    ajout du point ns de pxyd dans letree
 764 c -----
 765 c
 766 c entrees:
 767 c --------
 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
 772 c
 773 c modifies :
 774 c ----------
 775 c nbsomm : nombre actuel de points dans pxyd
 776 c
 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
 784 c      sinon
 785 c         letree(0:3,j) :-no pxyd des 1 \85a 4 points internes au triangle j
 786 c                         0  si pas de point
 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
 791 c
 792 c sorties :
 793 c ---------
 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)
 801 c
 802 c     depart de la racine
 803       ntrp = 1
 804 c
 805 c     recherche du triangle contenant le point pxyd(ns)
 806  1    ntrp = notrpt( pxyd(1,ns), pxyd, ntrp, letree )
 807 c
 808 c     existe t il un point libre
 809       do 10 i=0,3
 810          if( letree(i,ntrp) .eq. 0 ) then
 811 c           la place i est libre
 812             letree(i,ntrp) = -ns
 813             return
 814          endif
 815  10   continue
 816 c
 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
 821 c
 822 c     ajout du point ns
 823       goto 1
 824       end
 825
 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
 831 c
 832 c entrees:
 833 c --------
 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
 843 c   sinon
 844 c      letree(0:3,j) :-no pxyd des 1 a 4 points internes au triangle j
 845 c                      0  si pas de point
 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
 850 c
 851 c sorties :
 852 c ---------
 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:*)
 861       integer   lapile(1:64)
 862 c
 863 c     initialisation de la pile
 864 c     le triangle est empile
 865       lapile(1) = nt
 866       lhpile = 1
 867 c
 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
 872          notrva = 0
 873          lhpile = lhpile - 1
 874          return
 875       endif
 876 c
 877 c     le type du triangle ntr
 878       nty  = letree( 5, ntr )
 879 c     l'eventuel sur-triangle
 880       nsut = letree( 4, ntr )
 881 c
 882       if( nty .eq. 0 ) then
 883 c
 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 )
 888          goto 20
 889       endif
 890 c
 891 c     triangle ntr de type nty>0
 892       if( nosui3(nty) .eq. lar ) then
 893 c
 894 c        le triangle voisin par lar est le triangle 0
 895          lapile( lhpile ) = letree( 0, nsut )
 896          goto 20
 897       endif
 898 c
 899 c     triangle sans voisin direct => passage par le sur-triangle
 900       if( nsut .eq. 0 ) then
 901 c
 902 c        ntr est la racine => pas de triangle voisin par cette arete
 903          notrva = 0
 904          return
 905       else
 906 c
 907 c        le sur-triangle est empile
 908          lhpile = lhpile + 1
 909          lapile(lhpile) = nsut
 910          goto 10
 911       endif
 912 c
 913 c     descente aux sous-triangles selon la meme arete
 914  20   notrva = lapile( lhpile )
 915 c
 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
 920          return
 921       else
 922 c        le triangle a 4 sous-triangles
 923          if( lhpile .gt. 0 ) then
 924 c
 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 )
 930             else
 931 c              l'oppose est precedent(nty) de notrva
 932                notrva = letree( nopre3(nty) , notrva )
 933             endif
 934             goto 30
 935          endif
 936       endif
 937 c
 938 c     meme niveau dans l'arbre lhpile = 0
 939       end
 940
 941
 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
 947 c
 948 c entrees :
 949 c ---------
 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
 953 c
 954 c sortie :
 955 c --------
 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,
 967      %                  aire2,xc,yc,rot,
 968      %                  xy1(2),xy2(2),xy3(2),cetria(3)
 969 c
 970 c     le calcul de 2 fois l'aire du triangle
 971 c     attention l'ordre des 3 sommets est direct ou non
 972       x1  = xy1(1)
 973       x21 = xy2(1) - x1
 974       x31 = xy3(1) - x1
 975 c
 976       y1  = xy1(2)
 977       y21 = xy2(2) - y1
 978       y31 = xy3(2) - y1
 979 c
 980       aire2  = x21 * y31 - x31 * y21
 981 c
 982 c     recherche d'un test relatif peu couteux
 983 c     pour reperer la degenerescence du triangle
 984       if( abs(aire2) .le.
 985      %    epsurf*(abs(x21)+abs(x31))*(abs(y21)+abs(y31)) ) then
 986 c        triangle de qualite trop faible
 987          if( ierr .ge. 0 ) then
 988 c            nblgrc(nrerr) = 1
 989 c            kerr(1) = 'erreur cenced: triangle degenere'
 990 c            call lereur
 991             write(imprim,*) 'erreur cenced: triangle degenere'
 992             write(imprim,10000)  xy1,xy2,xy3,aire2
 993          endif
 994 10000 format( 3(' x=',g24.16,' y=',g24.16/),' aire*2=',g24.16)
 995          cetria(1) = 0d0
 996          cetria(2) = 0d0
 997          cetria(3) = 1d28
 998          ierr = 1
 999          return
1000       endif
1001 c
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)
1009 c
1010       xc = ( x1 + xy3(1) ) * 0.5d0 + rot * y31
1011       yc = ( y1 + xy3(2) ) * 0.5d0 - rot * x31
1012 c
1013       cetria(1) = xc
1014       cetria(2) = yc
1015 c
1016 c     le carre du rayon
1017       cetria(3) = (x1-xc) ** 2 + (y1-yc) ** 2
1018 c
1019 c     pas d'erreur rencontree
1020       ierr = 0
1021       end
1022
1023
1024       double precision function angled( p1, p2, p3 )
1025 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1026 c but :   calculer l'angle (p1p2,p1p3) en radians
1027 c -----
1028 c
1029 c entrees :
1030 c ---------
1031 c p1,p2,p3 : les 2 coordonnees des 3 sommets de l'angle
1032 c               sens direct pour une surface >0
1033 c sorties :
1034 c ---------
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
1041 c
1042 c     les cotes
1043       x21 = p2(1) - p1(1)
1044       y21 = p2(2) - p1(2)
1045       x31 = p3(1) - p1(1)
1046       y31 = p3(2) - p1(2)
1047 c
1048 c     longueur des cotes
1049       a1 = x21 * x21 + y21 * y21
1050       a2 = x31 * x31 + y31 * y31
1051       d  = sqrt( a1 * a2 )
1052       if( d .eq. 0 ) then
1053          angled = 0
1054          return
1055       endif
1056 c
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
1062          return
1063       else if( c .ge. 1.d0 ) then
1064 c        tilt sur apollo si acos( 1 + eps )
1065          angled = 0
1066          return
1067       endif
1068 c
1069       angled = acos( c )
1070       if( x21 * y31 - x31 * y21 .lt. 0 ) then
1071 c        demi plan inferieur
1072          angled = 8.d0 * atan( 1.d0 ) - angled
1073       endif
1074       end
1075
1076
1077       subroutine teajte( mxsomm, nbsomm, pxyd,   comxmi,
1078      %                   aretmx, mxtree, letree,
1079      %                   ierr )
1080 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1081 c but :    initialisation des tableaux letree
1082 c -----    ajout des sommets 1 a nbsomm (valeur en entree) dans letree
1083 c
1084 c entrees:
1085 c --------
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
1089 c
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)
1096 c
1097 c sorties:
1098 c --------
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
1107 c          sinon
1108 c             letree(0:3,j) :-no pxyd des 1 a 4 points internes au triangle j
1109 c                             0  si pas de point
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
1114 c
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
1126 c
1127 c     protection du nombre de sommets avant d'ajouter ceux de tetree
1128       nbsofr = nbsomm
1129       do 1 i = 1, nbsomm
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) )
1134  1    continue
1135 c
1136 c     creation de l'arbre tee
1137 c     =======================
1138 c     la premiere colonne vide de letree
1139       letree(0,0) = 2
1140 c     chainage des te vides
1141       do 4 i = 2 , mxtree
1142          letree(0,i) = i+1
1143  4    continue
1144       letree(0,mxtree) = 0
1145 c     les maxima des 2 indices de letree
1146       letree(1,0) = 8
1147       letree(2,0) = mxtree
1148 c
1149 c     la racine
1150 c     aucun point interne au triangle equilateral (te) 1
1151       letree(0,1) = 0
1152       letree(1,1) = 0
1153       letree(2,1) = 0
1154       letree(3,1) = 0
1155 c     pas de sur-triangle
1156       letree(4,1) = 0
1157       letree(5,1) = 0
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
1162 c
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 )
1169       do 60 k=1,2
1170          if( a(k) .lt. 1e-4 * s ) then
1171 c            nblgrc(nrerr) = 1
1172             write(imprim,*) 'tous les points sont alignes'
1173 c            call lereur
1174             ierr = 7
1175             return
1176          endif
1177  60   continue
1178 c
1179 c     le maximum des ecarts
1180       s = s + s
1181 c
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
1187 c
1188 c     le point nbsomm + 1 en bas a gauche
1189       nbsomm = nbsomm + 1
1190       pxyd(1,nbsomm) = (comxmi(1,1)+comxmi(1,2))*0.5d0 - arete*0.5d0
1191       pxyd(2,nbsomm) =  comxmi(2,1) - aretmx
1192       pxyd(3,nbsomm) = s
1193 c
1194 c     le point nbsomm + 2 en bas a droite
1195       nbsomm = nbsomm + 1
1196       pxyd(1,nbsomm) = pxyd(1,nbsomm-1) + arete
1197       pxyd(2,nbsomm) = pxyd(2,nbsomm-1)
1198       pxyd(3,nbsomm) = s
1199 c
1200 c     le point nbsomm + 3 sommet au dessus
1201       nbsomm = nbsomm + 1
1202       pxyd(1,nbsomm) = pxyd(1,nbsomm-2) + arete * 0.5d0
1203       pxyd(2,nbsomm) = pxyd(2,nbsomm-2) + arete * 0.5d0 * rac3
1204       pxyd(3,nbsomm) = s
1205 c
1206 c     ajout des sommets des lignes pour former letree
1207 c     ===============================================
1208       do 150 i=1,nbsofr
1209 c        ajout du point i de pxyd a letree
1210          call teajpt(  i, nbsomm, mxsomm, pxyd, letree,
1211      &                nt, ierr )
1212          if( ierr .ne. 0 ) return
1213  150  continue
1214 c
1215       return
1216       end
1217
1218
1219       subroutine tetaid( nutysu, dx, dy, longai, ierr )
1220 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1221 c but :     calculer la longueur de l'arete ideale en dx,dy
1222 c -----
1223 c
1224 c entrees:
1225 c --------
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!)
1231 c
1232 c sorties:
1233 c --------
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)
1242 c
1243       double precision  areteideale
1244       double precision  dx, dy, longai
1245       double precision  xyz(3), xyzd(3), d0
1246 c
1247       ierr = 0
1248       if( nutysu .gt. 0 ) then
1249          d0 = longai
1250 c        le point ou se calcule la longueur
1251          xyz(1) = dx
1252          xyz(2) = dy
1253 c        z pour le calcul de la longueur (inactif ici!)
1254          xyz(3) = 0d0
1255 c        la direction pour le calcul de la longueur (inactif ici!)
1256          xyzd(1) = 0d0
1257          xyzd(2) = 0d0
1258          xyzd(3) = 0d0
1259
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' )
1265             longai = -longai
1266          endif
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!' )
1271             ierr = 2
1272             longai = d0
1273          endif
1274       endif
1275       end
1276
1277
1278       subroutine tehote( nutysu,
1279      %                   nbarpi, mxsomm, nbsomm, pxyd,
1280      %                   comxmi, aretmx,
1281      %                   letree, mxqueu, laqueu,
1282      %                   ierr )
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
1286 c
1287 c entrees:
1288 c --------
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
1302 c
1303 c modifies :
1304 c ----------
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
1315 c          sinon
1316 c             letree(0:3,j) :-no pxyd des 1 a 4 points internes au triangle j
1317 c                             0  si pas de point
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
1322 c
1323 c auxiliaire :
1324 c ------------
1325 c laqueu : mxqueu entiers servant de queue pour le parcours de letree
1326 c
1327 c sorties:
1328 c --------
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)
1339 c
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:*)
1344 c
1345       integer           laqueu(1:mxqueu),lequeu
1346 c                       lequeu : entree dans la queue
1347 c                       lhqueu : longueur de la queue
1348 c                       gestion circulaire
1349 c
1350       integer           nuste(3)
1351       equivalence      (nuste(1),ns1),(nuste(2),ns2),(nuste(3),ns3)
1352 c
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
1360 c
1361       if( nutysu .gt. 0 ) then
1362 c
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
1366          do 1 i=1,nbsomm
1367 c           calcul de pxyzd(3,i)
1368             call tetaid( nutysu, pxyd(1,i), pxyd(2,i),
1369      %                   pxyd(3,i), ierr )
1370             if( ierr .ne. 0 ) goto 9999
1371  1       continue
1372 c
1373       else
1374 c
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
1381          do 3 i=1,nbarpi
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
1384             nte = 1
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?
1388             ns1 = letree(6,nte)
1389             ns2 = letree(7,nte)
1390             ns3 = letree(8,nte)
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,
1400      &                      ierr )
1401                if( ierr .ne. 0 ) return
1402                goto 2
1403             endif
1404  3       continue
1405       endif
1406 c
1407 c     le sous-triangle central de la racine est decoupe systematiquement
1408 c     ==================================================================
1409       nte = 2
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
1413          nbsom0 = nbsomm
1414          call te4ste( nbsomm, mxsomm, pxyd, nte, letree,
1415      %                ierr )
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
1421  4       continue
1422       endif
1423 c
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
1427 c
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
1434       ns1 = letree(6,1)
1435       ns2 = letree(7,1)
1436       ns3 = letree(8,1)
1437       a   = aretmx * 0.01d0
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
1449 c
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
1457       endif
1458 c
1459       nbs0   = nbsomm
1460       nbiter = -1
1461 c
1462 c     initialisation de la queue
1463   5   nbiter = nbiter + 1
1464       lequeu = 1
1465       lhqueu = 0
1466 c     la racine de letree initialise la queue
1467       laqueu(1) = 1
1468 c
1469 c     tant que la longueur de la queue est >=0 traiter le debut de queue
1470  10   if( lhqueu .ge. 0 ) then
1471 c
1472 c        le triangle te a traiter
1473          i   = lequeu - lhqueu
1474          if( i .le. 0 ) i = mxqueu + i
1475          nte = laqueu( i )
1476 c        la longueur de la queue est reduite
1477          lhqueu = lhqueu - 1
1478 c
1479 c        nte est il un sous-triangle feuille minimal ?
1480  15      if( letree(0,nte) .gt. 0 ) then
1481 c
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'
1485                ierr = 7
1486                return
1487             endif
1488             do 20 i=3,0,-1
1489 c              ajout du sous-triangle i
1490                lhqueu = lhqueu + 1
1491                lequeu = lequeu + 1
1492                if( lequeu .gt. mxqueu ) lequeu = lequeu - mxqueu
1493                laqueu( lequeu ) = letree( i, nte )
1494  20         continue
1495             goto 10
1496 c
1497          endif
1498 c
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 ?
1502          ns1 = letree(6,nte)
1503          ns2 = letree(7,nte)
1504          ns3 = letree(8,nte)
1505          if( pxyd(1,ns1) .gt. pxyd(1,ns2) ) then
1506             dmin = pxyd(1,ns2)
1507             dmax = pxyd(1,ns1)
1508          else
1509             dmin = pxyd(1,ns1)
1510             dmax = pxyd(1,ns2)
1511          endif
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
1515                dmin = pxyd(2,ns3)
1516                dmax = pxyd(2,ns1)
1517             else
1518                dmin = pxyd(2,ns1)
1519                dmax = pxyd(2,ns3)
1520             endif
1521             if( (yrmin .le. dmin .and. dmin .le. yrmax) .or.
1522      %          (yrmin .le. dmax .and. dmax .le. yrmax) ) then
1523 c
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
1529 c
1530                if( nutysu .eq. 0 ) then
1531 c
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
1542                      goto 15
1543                   endif
1544 c
1545                else
1546 c
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
1551                   do 28 i=1,3
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
1555                         nbsom0 = nbsomm
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),
1562      %                                  pxyd(3,j), ierr )
1563                            if( ierr .ne. 0 ) goto 9999
1564  27                     continue
1565                         goto 15
1566                      endif
1567  28               continue
1568                endif
1569 c
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              =================================================================
1573  29            do 30 i=1,3
1574 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
1583 c
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
1592                            goto 30
1593                         endif
1594                      endif
1595                   endif
1596 c
1597 c                 saut>1 => le triangle nte doit etre subdivise en 4 sous-triang
1598 c                 --------------------------------------------------------------
1599                   nbsom0 = nbsomm
1600                   call te4ste( nbsomm,mxsomm, pxyd, nte, letree,
1601      &                         ierr )
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),
1607      %                               pxyd(3,j), ierr )
1608                         if( ierr .ne. 0 ) goto 9999
1609  32                  continue
1610                   endif
1611                   goto 15
1612 c
1613  30            continue
1614             endif
1615          endif
1616          goto 10
1617       endif
1618       if( nbs0 .lt. nbsomm ) then
1619          nbs0 = nbsomm
1620          goto 5
1621       endif
1622       return
1623 c
1624 c     pb dans le calcul de la fonction taille_ideale
1625
1626  9999 write(imprim,*) 'pb dans le calcul de taille_ideale'
1627 c      nblgrc(nrerr) = 1
1628 c      kerr(1) = 'pb dans le calcul de taille_ideale'
1629 c      call lereur
1630       return
1631       end
1632
1633
1634       subroutine tetrte( comxmi, aretmx, nbarpi, mxsomm, pxyd,
1635      %                   mxqueu, laqueu, letree,
1636      %                   mosoar, mxsoar, n1soar, nosoar,
1637      %                   moartr, mxartr, n1artr, noartr, noarst,
1638      %                   ierr  )
1639 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1640 c but :    trianguler les triangles equilateraux feuilles et
1641 c -----    les points de la frontiere et les points internes imposes
1642 c
1643 c attention: la triangulation finale n'est pas de type delaunay!
1644 c
1645 c entrees:
1646 c --------
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
1654 c
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
1667 c          sinon
1668 c             letree(0:3,j) :-no pxyd des 1 a 4 points internes au triangle j
1669 c                             0  si pas de point
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
1674 c
1675 c modifies:
1676 c ---------
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
1683 c
1684 c auxiliaire :
1685 c ------------
1686 c laqueu : mxqueu entiers servant de queue pour le parcours de letree
1687 c
1688 c sorties:
1689 c --------
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)
1703 c
1704       double precision  pxyd(3,mxsomm)
1705       double precision  comxmi(3,2),aretmx,a,s,xrmin,xrmax,yrmin,yrmax
1706       double precision  dmin, dmax
1707 c
1708       integer           nosoar(mosoar,mxsoar),
1709      %                  noartr(moartr,mxartr),
1710      %                  noarst(mxsomm)
1711 c
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
1716 c
1717       integer           milieu(3), nutr(1:13)
1718 c
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
1721       ns1 = letree(6,1)
1722       ns2 = letree(7,1)
1723       ns3 = letree(8,1)
1724       a   = aretmx * 0.01d0
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
1736 c
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
1744       endif
1745 c
1746 c     initialisation du tableau noartr
1747       do 5 i=1,mxartr
1748 c        le numero de l'arete est inconnu
1749          noartr(1,i) = 0
1750 c        le chainage sur le triangle vide suivant
1751          noartr(2,i) = i+1
1752  5    continue
1753       noartr(2,mxartr) = 0
1754       n1artr = 1
1755 c
1756 c     parcours des te jusqu'a trianguler toutes les feuilles (triangles eq)
1757 c     =====================================================================
1758 c     initialisation de la queue sur les te
1759       ierr   = 0
1760       lequeu = 1
1761       lhqueu = 0
1762 c     la racine de letree initialise la queue
1763       laqueu(1) = 1
1764 c
1765 c     tant que la longueur de la queue est >=0 traiter le debut de queue
1766  10   if( lhqueu .ge. 0 ) then
1767 c
1768 c        le triangle te a traiter
1769          i   = lequeu - lhqueu
1770          if( i .le. 0 ) i = mxqueu + i
1771          nte = laqueu( i )
1772 c        la longueur est reduite
1773          lhqueu = lhqueu - 1
1774 c
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'
1780                ierr = 5
1781                return
1782             endif
1783             do 20 i=3,0,-1
1784 c              ajout du sous-triangle i
1785                lhqueu = lhqueu + 1
1786                lequeu = lequeu + 1
1787                if( lequeu .gt. mxqueu ) lequeu = lequeu - mxqueu
1788                laqueu( lequeu ) = letree( i, nte )
1789  20         continue
1790             goto 10
1791          endif
1792 c
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 ?
1796          ns1 = letree(6,nte)
1797          ns2 = letree(7,nte)
1798          ns3 = letree(8,nte)
1799          if( pxyd(1,ns1) .gt. pxyd(1,ns2) ) then
1800             dmin = pxyd(1,ns2)
1801             dmax = pxyd(1,ns1)
1802          else
1803             dmin = pxyd(1,ns1)
1804             dmax = pxyd(1,ns2)
1805          endif
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
1809                dmin = pxyd(2,ns3)
1810                dmax = pxyd(2,ns1)
1811             else
1812                dmin = pxyd(2,ns1)
1813                dmax = pxyd(2,ns3)
1814             endif
1815             if( (yrmin .le. dmin .and. dmin .le. yrmax) .or.
1816      %          (yrmin .le. dmax .and. dmax .le. yrmax) ) then
1817 c
1818 c              te minimal et interne au rectangle englobant
1819 c              --------------------------------------------
1820 c              recherche du nombre de niveaux entre nte et les te voisins
1821 c              par ses aretes
1822                nbmili = 0
1823                do 30 i=1,3
1824 c
1825 c                 a priori pas de milieu de l'arete i du te nte
1826                   milieu(i) = 0
1827 c
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 )
1843                         nbmili = nbmili + 1
1844                      endif
1845                   endif
1846 c
1847  30            continue
1848 c
1849 c              triangulation du te nte en fonction du nombre de ses milieux
1850                goto( 50, 100, 200, 300 ) , nbmili + 1
1851 c
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,
1857      %                      noarst,
1858      %                      nbtr,   nutr,   ierr )
1859                if( ierr .ne. 0 ) return
1860                goto 10
1861 c
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,
1867      %                      noarst,
1868      %                      nbtr,   nutr,   ierr )
1869                if( ierr .ne. 0 ) return
1870                goto 10
1871 c
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,
1877      %                      noarst,
1878      %                      nbtr,   nutr,   ierr )
1879                if( ierr .ne. 0 ) return
1880                goto 10
1881 c
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,
1887      %                      noarst,
1888      %                      nbtr,   nutr,   ierr )
1889                if( ierr .ne. 0 ) return
1890                goto 10
1891             endif
1892          endif
1893          goto 10
1894       endif
1895       end
1896
1897
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
1902 c
1903 c entrees:
1904 c --------
1905 c mosoar : nombre maximal d'entiers par arete dans le tableau nosoar
1906 c mxsoar : nombre maximal d'aretes frontalieres declarables
1907 c
1908 c modifies :
1909 c ----------
1910 c nosoar : numero des 2 sommets , no ligne, 2 triangles, chainages en +
1911 c          nosoar(lchain,i)=arete interne suivante
1912 c
1913 c sortie :
1914 c --------
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)
1922 c
1923 c     formation du chainage des aretes internes a echanger eventuellement
1924 c     recherche de la premiere arete non vide et non frontaliere
1925       do 10 na1=1,mxsoar
1926          if( nosoar(1,na1) .gt. 0 .and. nosoar(3,na1) .le. 0 ) goto 15
1927  10   continue
1928 c
1929 c     protection de la premiere arete non vide et non frontaliere
1930  15   na0 = na1
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
1935             na0 = na
1936          endif
1937  20   continue
1938 c
1939 c     la derniere arete interne n'a pas de suivante
1940       nosoar(lchain,na0) = 0
1941       end
1942
1943
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)
1953 c
1954 c entrees:
1955 c --------
1956 c pxyd   : tableau des x  y  distance_souhaitee de chaque sommet
1957 c
1958 c modifies :
1959 c ----------
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
1967 c
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),
1984      %                  noarst(*)
1985 c
1986 c     le nombre d'echanges de diagonales pour minimiser l'aire
1987       modifs = 0
1988       r0     = 0
1989 c
1990 c     la premiere arete du chainage des aretes a rendre delaunay
1991       na0 = n1ardv
1992 c
1993 c     tant que la pile des aretes a echanger eventuellement est non vide
1994 c     ==================================================================
1995  20   if( na0 .gt. 0 ) then
1996 c
1997 c        l'arete a traiter
1998          na  = na0
1999 c        la prochaine arete a traiter
2000          na0 = nosoar(lchain,na0)
2001 c
2002 c        l'arete est marquee traitee avec le numero -1
2003          nosoar(lchain,na) = -1
2004 c
2005 c        l'arete est elle active?
2006          if( nosoar(1,na) .eq. 0 ) goto 20
2007 c
2008 c        si arete frontaliere pas d'echange possible
2009          if( nosoar(3,na) .gt. 0 ) goto 20
2010 c
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
2013 c
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
2017 c
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
2023 c
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
2026 c
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
2034 c
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 )
2039 c
2040          if( abs(s34-s12) .gt. 1d-15*s34 ) goto 20
2041 c
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
2046          ierr = -1
2047          call cenced( pxyd(1,ns1), pxyd(1,ns2), pxyd(1,ns3), cetria,
2048      %                ierr )
2049          if( ierr .gt. 0 ) then
2050 c           ierr=1 si triangle degenere  => abandon
2051             goto 20
2052          endif
2053 c
2054          if( (cetria(1)-pxyd(1,ns4))**2+(cetria(2)-pxyd(2,ns4))**2
2055      %       .lt. cetria(3) ) then
2056 c
2057 c           protection contre une boucle infinie sur le meme cercle
2058             if( r0 .eq. cetria(3) ) goto 20
2059 c
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
2062 c
2063 cccc           les 2 triangles d'arete na sont effaces
2064 ccc            do 25 j=4,5
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 )
2069 ccc 25         continue
2070 c
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
2075             r0 = cetria(3)
2076 c
2077 c           l'arete na34 est marquee traitee
2078             nosoar(lchain,na34) = -1
2079             modifs = modifs + 1
2080 c
2081 c           les aretes internes peripheriques des 2 triangles sont enchainees
2082             do 60 j=4,5
2083                nt = nosoar(j,na34)
2084 cccc              trace du triangle nt
2085 ccc               call mttrtr( pxyd, nt, moartr, noartr, mosoar, nosoar,
2086 ccc     %                      ncoran, ncgric )
2087                do 50 i=1,3
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
2094                          na0 = n
2095                      endif
2096                   endif
2097  50            continue
2098  60         continue
2099             goto 20
2100          endif
2101 c
2102 c        retour en haut de la pile des aretes a traiter
2103          goto 20
2104       endif
2105       end
2106
2107
2108       subroutine terefr( nbarpi, pxyd,
2109      %                   mosoar, mxsoar, n1soar, nosoar,
2110      %                   moartr, n1artr, noartr, noarst,
2111      %                   mxarcf, n1arcf, noarcf, larmin, notrcf,
2112      %                   nbarpe, ierr )
2113 c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2114 c but :   recherche des aretes de la frontiere non dans la triangulation
2115 c -----   triangulation frontale pour les reobtenir
2116 c
2117 c         attention: le chainage lchain de nosoar devient celui des cf
2118 c
2119 c entrees:
2120 c --------
2121 c          le tableau nosoar
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
2131 c
2132 c modifies:
2133 c ---------
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
2148 c
2149 c
2150 c auxiliaires :
2151 c -------------
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
2156 c
2157 c sortie :
2158 c --------
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),
2169      %                  noartr(moartr,*),
2170      %                  noarst(*),
2171      %                  n1arcf(0:mxarcf),
2172      %                  noarcf(3,mxarcf),
2173      %                  larmin(mxarcf),
2174      %                  notrcf(mxarcf)
2175 c
2176 c     le nombre d'aretes de la frontiere non arete de la triangulation
2177       nbarpe = 0
2178 c
2179 c     initialisation du chainage des aretes des cf => 0 arete de cf
2180       do 10 narete=1,mxsoar
2181          nosoar( lchain, narete) = -1
2182  10   continue
2183 c
2184 c     boucle sur l'ensemble des aretes actuelles
2185 c     ==========================================
2186       do 30 narete=1,mxsoar
2187 c
2188          if( nosoar(3,narete) .gt. 0 ) then
2189 c           arete appartenant a une ligne => frontaliere
2190 c
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
2194                nbarpe = nbarpe + 1
2195 c
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))
2203 c
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,
2209      %                      ierr )
2210                if( ierr .ne. 0 ) return
2211 c
2212 c              fin du traitement de cette arete perdue et retrouvee
2213             endif
2214          endif
2215 c
2216  30   continue
2217       end
2218
2219
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
2229 c
2230 c entrees:
2231 c --------
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
2238 c          sommet frontalier
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
2260 c
2261 c sorties:
2262 c --------
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),
2274      %                  noarst(*)
2275       integer           letrsu(1:ndtri0)
2276       double precision  dmin
2277 c
2278 c     les triangles sont a priori non marques
2279       do 5 nt=1,ndtri0
2280          letrsu(nt) = 0
2281  5    continue
2282 c
2283 c     les aretes sont marquees non chainees
2284       do 10 noar1=1,mxsoar
2285          nosoar(6,noar1) = -2
2286  10   continue
2287 c
2288 c     recherche du sommet de la triangulation de plus petite abscisse
2289 c     ===============================================================
2290       ntmin = 0
2291       dmin  = 1d38
2292       do 20 i=1,nbsomm
2293          if( pxyd(1,i) .lt. dmin ) then
2294 c           le nouveau minimum
2295             noar1 = noarst(i)
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
2300                   dmin  = pxyd(1,i)
2301                   ntmin = i
2302                endif
2303             endif
2304          endif
2305  20   continue
2306 c
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
2312 c         nblgrc(nrerr) = 1
2313 c         kerr(1) = 'pas de triangle d''abscisse minimale'
2314 c         call lereur
2315          write(imprim,*) 'pas de triangle d''abscisse minimale'
2316          ierr = 2
2317          goto 9990
2318       endif
2319 c
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
2330       nosoar(6,na0) = na1
2331 c     la troisieme n'a pas de suivante
2332       nosoar(6,na1) = 0
2333 c
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     =======================================================
2339       ligne0 = 0
2340       ligne  = -123 456 789
2341 c
2342  40   if( noar1 .ne. 0 ) then
2343 c
2344 c        l'arete noar1 du tableau nosoar est a traiter
2345 c        ---------------------------------------------
2346          noar = noar1
2347 c        l'arete suivante devient la premiere a traiter ensuite
2348          noar1 = nosoar(6,noar1)
2349 c        l'arete noar est traitee
2350          nosoar(6,noar) = -3
2351 c
2352          do 60 i=4,5
2353 c
2354 c           l'un des 2 triangles de l'arete
2355             nt = nosoar(i,noar)
2356             if( nt .gt. 0 ) then
2357 c
2358 c              triangle deja traite pour une ligne anterieure?
2359                if(     letrsu(nt)  .ne. 0      .and.
2360      %             abs(letrsu(nt)) .ne. ligne ) goto 60
2361 c
2362 cccc              trace du triangle nt en couleur ligne0
2363 ccc               call mttrtr( pxyd,   nt, moartr, noartr, mosoar, nosoar,
2364 ccc     %                      ligne0, ncnoir )
2365 c
2366 c              le triangle est marque avec la valeur de ligne
2367                letrsu(nt) = ligne
2368 c
2369 c              chainage eventuel des autres aretes de ce triangle
2370 c              si ce n'est pas encore fait
2371                do 50 j=1,3
2372 c
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
2376 c
2377 c                 le numero de 1 a nblftr dans nulftr de la ligne de l'arete
2378                   nl = nosoar(3,na)
2379 c
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
2383 c
2384                      if( nl .eq. ligne0 ) goto 50
2385 c
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
2390                         nt2 = 5
2391                      else
2392                         nt2 = 4
2393                      endif
2394                      nt2 = abs( nosoar(nt2,na) )
2395                      if( nt2 .gt. 0 ) then
2396 c
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
2400                            lsigne = -1
2401                         else
2402                            lsigne =  1
2403                         endif
2404                         letrsu(nt2) = lsigne * nl
2405 c
2406 c                       temoin de ligne a traiter ensuite dans nulftr
2407                         nulftr(nl) = -abs( nulftr(nl) )
2408 c
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 )
2413 c
2414 c                       l'arete est traitee
2415                         nosoar(6,na) = -3
2416 c
2417                      endif
2418 c
2419 c                    l'arete est traitee
2420                      goto 50
2421 c
2422                   endif
2423 c
2424 c                 arete non traitee => elle est chainee
2425                   nosoar(6,na) = noar1
2426                   noar1 = na
2427 c
2428  50            continue
2429 c
2430             endif
2431  60      continue
2432 c
2433          goto 40
2434       endif
2435 c     les triangles de la ligne fermee ont tous ete marques
2436 c     plus d'arete chainee
2437 c
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
2442  70   continue
2443 c     plus de ligne fermee a traiter
2444       goto 110
2445 c
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)
2452       do 90 nt2=1,ndtri0
2453          if( abs(letrsu(nt2)) .eq. nl ) goto 92
2454  90   continue
2455 c
2456 c     recherche de l'arete j du triangle nt2 avec ce numero de ligne nl
2457  92   do 95 j=1,3
2458 c
2459 c        le numero de l'arete j du triangle dans nosoar
2460          noar1 = 0
2461          na0   = abs( noartr(j,nt2) )
2462          if( nl .eq. nosoar(3,na0) ) then
2463 c
2464 c           na0 est l'arete de ligne nl
2465 c           l'arete suivante du triangle nt2
2466             i   = mod(j,3) + 1
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
2471                noar1 = na1
2472 c              pas de suivante dans ce chainage
2473                nosoar(6,na1) = 0
2474             else
2475                na1 = 0
2476             endif
2477 c
2478 c           l'eventuelle seconde arete suivante
2479             i  = mod(i,3) + 1
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
2484                   noar1 = na
2485                   nosoar(6,na) = 0
2486                else
2487 c                 2 aretes a chainer
2488                   noar1 = na
2489                   nosoar(6,na) = na1
2490                endif
2491             endif
2492 c
2493             if( noar1 .gt. 0 ) then
2494 c
2495 c              il existe au moins une arete a visiter pour ligne
2496 c              marquage des triangles internes a la ligne nl
2497                ligne  = letrsu(nt2)
2498                ligne0 = nl
2499                goto 40
2500 c
2501             else
2502 c
2503 c              nt2 est le seul triangle de la ligne fermee
2504                goto 65
2505 c
2506             endif
2507          endif
2508  95   continue
2509 c
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
2518  170  continue
2519 c
2520 c     les triangles externes sont marques vides dans le tableau noartr
2521 c     ================================================================
2522       nbtria = 0
2523       do 200 nt=1,ndtri0
2524 c
2525          if( letrsu(nt) .le. 0 ) then
2526 c
2527 c           triangle nt externe
2528             if( noartr(1,nt) .ne. 0 ) then
2529 c              la premiere arete est annulee
2530                noartr(1,nt) = 0
2531 c              le triangle nt est considere comme etant vide
2532                noartr(2,nt) = n1artr
2533                n1artr = nt
2534             endif
2535 c
2536          else
2537 c
2538 c           triangle nt interne
2539             nbtria = nbtria + 1
2540             letrsu(nt) = nbtria
2541 c
2542 c           marquage des 3 sommets du triangle nt
2543             do 190 i=1,3
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
2555  190        continue
2556 c
2557          endif
2558 c
2559  200  continue
2560 c     ici tout sommet externe ns verifie nslign(ns)=-1
2561 c
2562 c     les triangles externes sont mis a zero dans nosoar
2563 c     ==================================================
2564       do 300 noar=1,mxsoar
2565 c
2566          if( nosoar(1,noar) .gt. 0 ) then
2567 c
2568 c           le second triangle de l'arete noar
2569             nt = nosoar(5,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
2574             endif
2575 c
2576 c           le premier triangle de l'arete noar
2577             nt = nosoar(4,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)
2586                      nosoar(5,noar)=0
2587                   else
2588                      nosoar(4,noar)=0
2589                   endif
2590                endif
2591             endif
2592          endif
2593 c
2594  300  continue
2595 c
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)
2599  9991 continue
2600       return
2601       end
2602
2603
2604
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
2609 c -----
2610 c         limite: un camembert de centre ns entame 2 fois
2611 c                 ne donne que l'une des parties
2612 c
2613 c entrees:
2614 c --------
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
2624 c
2625 c sorties :
2626 c --------
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,*),
2639      %                  nosoar(mosoar,*),
2640      %                  noarst(*)
2641       integer           lapile(1:mxpile)
2642       integer           nosotr(3)
2643 c
2644 c     la premiere arete de sommet ns
2645       nar = noarst( ns )
2646       if( nar .le. 0 ) then
2647          write(imprim,*) 'trp1st: sommet',ns,' sans arete'
2648          goto 9999
2649       endif
2650 c
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)
2655          goto 9999
2656       endif
2657 c
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'
2662          goto 9999
2663       endif
2664 c
2665 c     le triangle est il interne?
2666       if( noartr(1,nt0) .eq. 0 ) goto 9999
2667 c
2668 c     le numero des 3 sommets du triangle nt0 dans le sens direct
2669       call nusotr( nt0, mosoar, nosoar, moartr, noartr, nosotr )
2670 c
2671 c     reperage du sommet ns dans le triangle nt0
2672       do 5 nar=1,3
2673          if( nosotr(nar) .eq. ns ) goto 10
2674  5    continue
2675       nta = nt0
2676       goto 9995
2677 c
2678 c     ns retrouve : le triangle nt0 est empile
2679  10   lhpile = 1
2680       lapile(1) = nt0
2681       nta = nt0
2682 c
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)
2691       else
2692          nt1 = nosoar(4,noar)
2693       endif
2694 c
2695 c     la boucle sur les triangles nt1 de sommet ns dans le sens indirect
2696 c     ==================================================================
2697       if( nt1 .gt. 0 ) then
2698 c
2699          if( noartr(1,nt1) .eq. 0 ) goto 30
2700 c
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 )
2705 c
2706 c        reperage du sommet ns dans nt1
2707          do 20 nar=1,3
2708             if( nosotr(nar) .eq. ns ) goto 25
2709  20      continue
2710          nta = nt1
2711          goto 9995
2712 c
2713 c        nt1 est empile
2714  25      if( lhpile .ge. mxpile ) goto 9990
2715          lhpile = lhpile + 1
2716          lapile(lhpile) = nt1
2717 c
2718 c        le triangle nt1 de l'autre cote de l'arete de sommet ns
2719 c        sauvegarde du precedent triangle dans nta
2720          nta  = nt1
2721          noar = abs( noartr(nar,nt1) )
2722          if( nosoar(4,noar) .eq. nt1 ) then
2723             nt1 = nosoar(5,noar)
2724          else
2725             nt1 = nosoar(4,noar)
2726          endif
2727          if( nt1 .le. 0   ) goto 30
2728 c        le triangle suivant est a l'exterieur
2729          if( nt1 .ne. nt0 ) goto 15
2730 c
2731 c        recherche terminee par arrivee sur nt0
2732 c        les triangles forment un "cercle" de "centre" ns
2733          return
2734 c
2735       endif
2736 c
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
2742  30   lhpile = 1
2743       lapile(lhpile) = nta
2744 c
2745 c     le numero des 3 sommets du triangle nta dans le sens direct
2746       call nusotr( nta, mosoar, nosoar, moartr, noartr, nosotr )
2747       do 32 nar=1,3
2748          if( nosotr(nar) .eq. ns ) goto 33
2749  32   continue
2750       goto 9995
2751 c
2752 c     l'arete qui precede (rotation / ns dans le sens direct)
2753  33   if( nar .eq. 1 ) then
2754          nar = 3
2755       else
2756          nar = nar - 1
2757       endif
2758 c
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)
2763       else
2764          nt1 = nosoar(4,noar)
2765       endif
2766       if( nt1 .le. 0 ) then
2767 c        un seul triangle contient ns
2768          goto 70
2769       endif
2770 c
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
2774 c
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 )
2778 c
2779 c     reperage du sommet ns dans nt1
2780       do 50 nar=1,3
2781          if( nosotr(nar) .eq. ns ) goto 60
2782  50   continue
2783       nta = nt1
2784       goto 9995
2785 c
2786 c     nt1 est empile
2787  60   if( lhpile .ge. mxpile ) goto 9990
2788       lhpile = lhpile + 1
2789       lapile(lhpile) = nt1
2790 c
2791 c     l'arete qui precede dans le sens direct
2792       if( nar .eq. 1 ) then
2793          nar = 3
2794       else
2795          nar = nar - 1
2796       endif
2797 c
2798 c     l'arete de sommet ns dans nosoar
2799       noar = abs( noartr(nar,nt1) )
2800 c
2801 c     le triangle voisin de nta dans le sens direct
2802       nta  = nt1
2803       if( nosoar(4,noar) .eq. nt1 ) then
2804          nt1 = nosoar(5,noar)
2805       else
2806          nt1 = nosoar(4,noar)
2807       endif
2808       nta = nt1
2809       if( nt1 .gt. 0 ) goto 40
2810 c
2811 c     butee sur le trou => fin des triangles de sommet ns
2812 c     ----------------------------------------------------
2813  70   lhpile = -lhpile
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
2816       return
2817 c
2818 c     saturation de la pile des triangles
2819 c     -----------------------------------
2820  9990 write(imprim,*) 'trp1st:saturation pile des triangles autour ',
2821      %'sommet',ns
2822       goto 9999
2823 c
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
2828 c
2829  9999 lhpile = 0
2830       return
2831       end
2832
2833
2834
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)
2839 c
2840 c entrees:
2841 c --------
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
2849 c
2850 c sorties:
2851 c --------
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)
2857 c
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
2861          nosotr(1) = 1
2862          nosotr(2) = 2
2863       else
2864          nosotr(1) = 2
2865          nosotr(2) = 1
2866          na = -na
2867       endif
2868       nosotr(1) = nosoar( nosotr(1), na )
2869       nosotr(2) = nosoar( nosotr(2), na )
2870 c
2871 c     l'arete suivante
2872       na = abs( noartr(2,nt) )
2873 c
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)
2878       endif
2879       end
2880
2881
2882       subroutine tesusp( nbarpi, pxyd,   noarst,
2883      %                   mosoar, mxsoar, n1soar, nosoar,
2884      %                   moartr, mxartr, n1artr, noartr,
2885      %                   mxarcf, n1arcf, noarcf, larmin, notrcf, liarcf,
2886      %                   nbstsu, ierr )
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
2891 c
2892 c         attention: le chainage lchain de nosoar devient celui des cf
2893 c
2894 c entrees:
2895 c --------
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
2905 c
2906 c modifies:
2907 c ---------
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
2922 c
2923 c
2924 c auxiliaires :
2925 c -------------
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
2931 c
2932 c sortie :
2933 c --------
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
2944 c
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),
2949      %                  noartr(moartr,*),
2950      %                  noarst(*),
2951      %                  n1arcf(0:mxarcf),
2952      %                  noarcf(3,mxarcf),
2953      %                  larmin(mxarcf),
2954      %                  notrcf(mxarcf),
2955      %                  liarcf(mxarcf)
2956 c
2957       integer           nosotr(3)
2958       equivalence      (nosotr(1),ns1), (nosotr(2),ns2),
2959      %                 (nosotr(3),ns3)
2960 c
2961 c     le nombre de sommets de te supprimes
2962       nbstsu = 0
2963 c
2964 c     initialisation du chainage des aretes des cf => 0 arete de cf
2965       do 10 narete=1,mxsoar
2966          nosoar( lchain, narete ) = -1
2967  10   continue
2968 c
2969 c     boucle sur l'ensemble des sommets frontaliers ou points internes
2970 c     ================================================================
2971       do 100 ns = 1, nbarpi
2972 c
2973 cccc        le nombre de sommets supprimes pour ce sommet ns
2974 ccc         nbsuns = 0
2975 c
2976 c        la qualite minimale au dessous de laquelle le point proche
2977 c        interne est supprime
2978          quaopt = quamal
2979 c
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'
2985             ierr = 11
2986             return
2987          endif
2988 c
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
2996             nbtrcf = -nbtrcf
2997          endif
2998 c
2999 c        boucle sur les triangles de l'etoile du sommet ns
3000          quamin = 2.0
3001          do 20 i=1,nbtrcf
3002 c
3003 c           le numero des 3 sommets du triangle nt
3004             nt = notrcf(i)
3005             call nusotr( nt, mosoar, nosoar, moartr, noartr,
3006      %                   nosotr )
3007 c           nosotr(1:3) est en equivalence avec ns1, ns2, ns3
3008 c
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
3012                quamin = qualit
3013                ntqmin = nt
3014             endif
3015  20      continue
3016 c
3017 c        bilan sur la qualite des triangles de sommet ns
3018          if( quamin .lt. quaopt ) then
3019 c
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,
3024      %                   nosotr )
3025             nste   = 0
3026             quamin = 1e28
3027             do 30 j=1,3
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
3032