src/MEFISTO2/aptrte.cxx

   1 //  MEFISTO2: a library to compute 2D triangulation from segmented boundaries
   2 //
   3 //
   4 //  Copyright (C) 2006  Laboratoire J.-L. Lions UPMC Paris
   5 //
   6 //  This library is free software; you can redistribute it and/or
   7 //  modify it under the terms of the GNU Lesser General Public
   8 //  License as published by the Free Software Foundation; either
   9 //  version 2.1 of the License.
  10 //
  11 //  This library is distributed in the hope that it will be useful,
  12 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 //  Lesser General Public License for more details.
  15 //
  16 //  You should have received a copy of the GNU Lesser General Public
  17 //  License along with this library; if not, write to the Free Software
  18 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 //  See http://www.ann.jussieu.fr/~perronne or email Perronnet@ann.jussieu.fr
  21 //
  22 //  File   : aptrte.cxx   le C++ de l'appel du trianguleur plan
  23 //  Module : SMESH
  24 //  Author : Alain PERRONNET
  25 //  Date   : 13 novembre 2006
  26
  27 #include "Rn.h"
  28 #include "aptrte.h"
  29 #include "utilities.h"
  30
  31 using namespace std;
  32
  33 extern "C"
  34 {
  35   R aretemaxface_;
  36   MEFISTO2D_EXPORT
  37     R
  38   #ifdef WIN32
  39       __stdcall
  40   #endif
  41       areteideale()//( R3 xyz, R3 direction )
  42   {
  43     return aretemaxface_;
  44   }
  45 }
  46 //calcul de la longueur ideale de l'arete au sommet xyz (z ici inactif)
  47 //dans la direction donnee
  48 //a ajuster pour chaque surface plane et selon l'entier notysu (voir plus bas)
  49
  50
  51 static double cpunew, cpuold=0;
  52
  53 void
  54 #ifdef WIN32
  55               __stdcall
  56 #endif
  57 tempscpu_( double & tempsec )
  58 //Retourne le temps CPU utilise en secondes
  59 {
  60   tempsec = ( (double) clock() ) / CLOCKS_PER_SEC;
  61   //MESSAGE( "temps cpu=" << tempsec );
  62 }
  63
  64
  65 void
  66 #ifdef WIN32
  67               __stdcall
  68 #endif
  69 deltacpu_( R & dtcpu )
  70 //Retourne le temps CPU utilise en secondes depuis le precedent appel
  71 {
  72   tempscpu_( cpunew );
  73   dtcpu  = R( cpunew - cpuold );
  74   cpuold = cpunew;
  75   //MESSAGE( "delta temps cpu=" << dtcpu );
  76   return;
  77 }
  78
  79
  80 void  aptrte( Z   nutysu, R      aretmx,
  81               Z   nblf,   Z  *   nudslf,  R2 * uvslf,
  82               Z   nbpti,  R2 *   uvpti,
  83               Z & nbst,   R2 * & uvst,
  84               Z & nbt,    Z  * & nust,
  85               Z & ierr )
  86 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  87 // but : appel de la triangulation par un arbre-4 recouvrant
  88 // ----- de triangles equilateraux
  89 //       le contour du domaine plan est defini par des lignes fermees
  90 //       la premiere ligne etant l'enveloppe de toutes les autres
  91 //       la fonction areteideale(s,d) donne la taille d'arete
  92 //       au point s dans la direction (actuellement inactive) d
  93 //       des lors toute arete issue d'un sommet s devrait avoir une longueur
  94 //       comprise entre 0.65 areteideale_(s,d) et 1.3 areteideale_(s,d)
  95 //
  96 //Attention:
  97 //  Les tableaux uvslf et uvpti sont supposes ne pas avoir de sommets identiques!
  98 //  De meme, un sommet d'une ligne fermee ne peut appartenir a une autre ligne fermee
  99 //
 100 // entrees:
 101 // --------
 102 // nutysu : numero de traitement de areteideale_(s,d) selon le type de surface
 103 //          0 pas d'emploi de la fonction areteideale_() et aretmx est active
 104 //          1 il existe une fonction areteideale_(s,d)
 105 //            dont seules les 2 premieres composantes de uv sont actives
 106 //          ... autres options a definir ...
 107 // aretmx : longueur maximale des aretes de la future triangulation
 108 // nblf   : nombre de lignes fermees de la surface
 109 // nudslf : numero du dernier sommet de chacune des nblf lignes fermees
 110 //          nudslf(0)=0 pour permettre la difference sans test
 111 //          Attention le dernier sommet de chaque ligne est raccorde au premier
 112 //          tous les sommets et les points internes ont des coordonnees
 113 //          UV differentes <=> Pas de point double!
 114 // uvslf  : uv des nudslf(nblf) sommets des lignes fermees
 115 // nbpti  : nombre de points internes futurs sommets de la triangulation
 116 // uvpti  : uv des points internes futurs sommets de la triangulation
 117 //
 118 // sorties:
 119 // --------
 120 // nbst   : nombre de sommets de la triangulation finale
 121 // uvst   : coordonnees uv des nbst sommets de la triangulation
 122 // nbt    : nombre de triangles de la triangulation finale
 123 // nust   : 4 numeros dans uvst des sommets des nbt triangles
 124 //          s1, s2, s3, 0: no dans uvst des 3 sommets et 0 car quadrangle!
 125 // ierr   : 0 si pas d'erreur
 126 //        > 0 sinon
 127 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 128 // auteur : Alain Perronnet  Laboratoire J.-L. LIONS Paris UPMC mars 2006
 129 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 130 {
 131   Z  nbsttria=4; //Attention: 4 sommets stockes par triangle
 132                  //no st1, st2, st3, 0 (non quadrangle)
 133
 134   R  d, tcpu=0;
 135   R3 direction=R3(0,0,0);  //direction pour areteideale() inactive ici!
 136   Z  nbarfr=nudslf[nblf];  //nombre total d'aretes des lignes fermees
 137   Z  mxtrou = Max( 1024, nblf );  //nombre maximal de trous dans la surface
 138
 139   R3 *mnpxyd=NULL;
 140   Z  *mnsoar=NULL, mosoar=7, mxsoar, n1soar; //le hachage des aretes
 141   Z  *mnartr=NULL, moartr=3, mxartr, n1artr; //le no des 3 aretes des triangles
 142   Z  *mntree=NULL, motree=9, mxtree; //L'arbre 4 de TE et nombre d'entiers par TE
 143   Z  *mnqueu=NULL, mxqueu;
 144   Z  *mn1arcf=NULL;
 145   Z  *mnarcf=NULL, mxarcf;
 146   Z  *mnarcf1=NULL;
 147   Z  *mnarcf2=NULL;
 148   Z  *mnarcf3=NULL;
 149   Z  *mntrsu=NULL;
 150   Z  *mnslig=NULL;
 151   Z  *mnarst=NULL;
 152   Z  *mnlftr=NULL;
 153
 154   R3 comxmi[2];            //coordonnees UV Min et Maximales
 155   R  aremin, aremax;       //longueur minimale et maximale des aretes
 156   R  airemx;               //aire maximale souhaitee d'un triangle
 157   R  quamoy, quamin;
 158
 159   Z  noar0, noar, na;
 160   Z  i, l, n, ns, ns0, ns1, ns2, nosotr[3], nt;
 161   Z  mxsomm, nbsomm, nbarpi, nbarli, ndtri0, mn;
 162   Z  moins1=-1;
 163   Z  nuds = 0;
 164
 165   // initialisation du temps cpu
 166   deltacpu_( d );
 167   ierr = 0;
 168
 169   // quelques reservations de tableaux pour faire les calculs
 170   // ========================================================
 171   // declaration du tableau des coordonnees des sommets de la frontiere
 172   // puis des sommets internes ajoutes
 173   // majoration empirique du nombre de sommets de la triangulation
 174   i = 4*nbarfr/10;
 175   mxsomm = Max( 20000, 64*nbpti+i*i );
 176   MESSAGE( "APTRTE: Debut de la triangulation plane avec " );
 177   MESSAGE( "nutysu=" << nutysu << "  aretmx=" << aretmx
 178            << "  mxsomm=" << mxsomm );
 179   MESSAGE( nbarfr << " sommets sur la frontiere et " << nbpti << " points internes");
 180
 181  NEWDEPART:
 182   //mnpxyd( 3, mxsomm ) les coordonnees UV des sommets et la taille d'arete aux sommets
 183   if( mnpxyd!=NULL ) delete [] mnpxyd;
 184   mnpxyd = new R3[mxsomm];
 185   if( mnpxyd==NULL ) goto ERREUR;
 186
 187   // le tableau mnsoar des aretes des triangles
 188   // 1: sommet 1 dans pxyd,
 189   // 2: sommet 2 dans pxyd,
 190   // 3: numero de 1 a nblf de la ligne qui supporte l'arete
 191   // 4: numero dans mnartr du triangle 1 partageant cette arete,
 192   // 5: numero dans mnartr du triangle 2 partageant cette arete,
 193   // 6: chainage des aretes frontalieres ou internes ou
 194   //    des aretes simples des etoiles de triangles,
 195   // 7: chainage du hachage des aretes
 196   // nombre d'aretes = 3 ( nombre de sommets - 1 + nombre de trous )
 197   // pour le hachage des aretes mxsoar doit etre > 3*mxsomm!
 198   // h(ns1,ns2) = min( ns1, ns2 )
 199   if( mnsoar!=NULL ) delete [] mnsoar;
 200   mxsoar = 3 * ( mxsomm + mxtrou );
 201   mnsoar = new Z[mosoar*mxsoar];
 202   if( mnsoar==NULL ) goto ERREUR;
 203   //initialiser le tableau mnsoar pour le hachage des aretes
 204   insoar( mxsomm, mosoar, mxsoar, n1soar, mnsoar );
 205
 206   // mnarst( mxsomm ) numero mnsoar d'une arete pour chacun des sommets
 207   if( mnarst!=NULL ) delete [] mnarst;
 208   mnarst = new Z[1+mxsomm];
 209   if( mnarst==NULL ) goto ERREUR;
 210   n = 1+mxsomm;
 211   azeroi( n, mnarst );
 212
 213   // mnslig( mxsomm ) no de sommet dans sa ligne pour chaque sommet frontalier
 214   //               ou no du point si interne forc'e par l'utilisateur
 215   //               ou  0 si interne cree par le module
 216   if( mnslig!=NULL ) delete [] mnslig;
 217   mnslig = new Z[mxsomm];
 218   if( mnslig==NULL ) goto ERREUR;
 219   azeroi( mxsomm, mnslig );
 220
 221   // initialisation des aretes frontalieres de la triangulation future
 222   // renumerotation des sommets des aretes des lignes pour la triangulation
 223   // mise a l'echelle des coordonnees des sommets pour obtenir une
 224   // meilleure precision lors des calculs + quelques verifications
 225   // boucle sur les lignes fermees qui forment la frontiere
 226   // ======================================================================
 227   noar = 0;
 228   aremin = 1e100;
 229   aremax = 0;
 230
 231   for (n=1; n<=nblf; n++)
 232   {
 233     //l'initialisation de la premiere arete de la ligne n dans la triangulation
 234     //-------------------------------------------------------------------------
 235     //le sommet ns0 est le numero de l'origine de la ligne
 236     ns0 = nudslf[n-1];
 237     mnpxyd[ns0].x = uvslf[ns0].x;
 238     mnpxyd[ns0].y = uvslf[ns0].y;
 239     mnpxyd[ns0].z = areteideale();//( mnpxyd[ns0], direction );
 240 //     MESSAGE("Sommet " << ns0 << ": " << mnpxyd[ns0].x
 241 //       << " " << mnpxyd[ns0].y << " longueur arete=" << mnpxyd[ns0].z);
 242
 243     //carre de la longueur de l'arete 1 de la ligne fermee n
 244     d = pow( uvslf[ns0+1].x - uvslf[ns0].x, 2 )
 245       + pow( uvslf[ns0+1].y - uvslf[ns0].y, 2 ) ;
 246     aremin = Min( aremin, d );
 247     aremax = Max( aremax, d );
 248
 249     //le numero des 2 sommets (ns1,ns2) de la premiere arete de la ligne
 250     //initialisation de la 1-ere arete ns1-ns1+1 de cette ligne fermee n
 251     //le numero des 2 sommets ns1 ns2 de la 1-ere arete
 252     //Attention: les numeros ns debutent a 1 (ils ont >0)
 253     //           les tableaux c++ demarrent a zero!
 254     //           les tableaux fortran demarrent ou l'on veut!
 255     ns0++;
 256     ns1 = ns0;
 257     ns2 = ns1+1;
 258
 259      //le numero n de la ligne du sommet et son numero ns1 dans la ligne
 260     mnslig[ns0-1] = 1000000 * n + ns1-nudslf[n-1];
 261     fasoar( ns1, ns2, moins1, moins1, n,
 262              mosoar, mxsoar, n1soar, mnsoar, mnarst,
 263              noar0,  ierr );
 264     //pas de test sur ierr car pas de saturation possible a ce niveau
 265
 266     //le pointeur dans le hachage sur la premiere arete de la ligne fermee n
 267     //mndalf[n] = noar0;
 268
 269     //la nouvelle arete est la suivante de l'arete definie juste avant
 270     if( noar > 0 )
 271       mnsoar[mosoar * noar - mosoar + 5] = noar0;
 272
 273     //l'initialisation des aretes suivantes de la ligne dans la triangulation
 274     //-----------------------------------------------------------------------
 275     nbarli = nudslf[n] - nudslf[n-1];  //nombre d'aretes=sommets de la ligne n
 276     for (i=2; i<=nbarli; i++)
 277     {
 278       ns1 = ns2; //le numero de l'arete et le numero du premier sommet de l'arete
 279       if( i < nbarli )
 280         //nbs+1 est le 2-eme sommet de l'arete i de la ligne fermee n
 281         ns2 = ns1+1;
 282       else
 283         //le 2-eme sommet de la derniere arete est le premier sommet de la ligne
 284         ns2 = ns0;
 285
 286       //l'arete precedente est dotee de sa suivante:celle cree ensuite
 287       //les 2 coordonnees du sommet ns2 de la ligne
 288       ns = ns1 - 1;
 289 //debut ajout  5/10/2006  ................................................
 290       nuds = Max( nuds, ns );   //le numero du dernier sommet traite
 291 //fin   ajout  5/10/2006  ................................................
 292       mnpxyd[ns].x = uvslf[ns].x;
 293       mnpxyd[ns].y = uvslf[ns].y;
 294       mnpxyd[ns].z = areteideale();//( mnpxyd[ns], direction );
 295 //       MESSAGE("Sommet " << ns << ": " << mnpxyd[ns].x
 296 //         << " " << mnpxyd[ns].y << " longueur arete=" << mnpxyd[ns].z);
 297
 298       //carre de la longueur de l'arete
 299       d = pow( uvslf[ns2-1].x - uvslf[ns1-1].x, 2)
 300         + pow( uvslf[ns2-1].y - uvslf[ns1-1].y, 2);
 301       aremin = Min( aremin, d );
 302       aremax = Max( aremax, d );
 303
 304 //debut ajout du 5/10/2006  .............................................
 305       //la longueur de l'arete ns1-ns2
 306       d = sqrt( d );
 307       //longueur arete = Min ( aretmx, aretes incidentes )
 308       mnpxyd[ns   ].z = Min( mnpxyd[ns   ].z, d );
 309       mnpxyd[ns2-1].z = Min( mnpxyd[ns2-1].z, d );
 310 //fin ajout du 5/10/2006  ...............................................
 311
 312       //le numero n de la ligne du sommet et son numero ns1 dans la ligne
 313       mnslig[ns] = 1000000 * n + ns1-nudslf[n-1];
 314
 315       //ajout de l'arete dans la liste
 316       fasoar( ns1, ns2, moins1, moins1, n,
 317                mosoar, mxsoar, n1soar, mnsoar,
 318                mnarst, noar, ierr );
 319       //pas de test sur ierr car pas de saturation possible a ce niveau
 320
 321       //chainage des aretes frontalieres en position 6 du tableau mnsoar
 322       //la nouvelle arete est la suivante de l'arete definie juste avant
 323       mnsoar[ mosoar * noar0 - mosoar + 5 ] = noar;
 324       noar0 = noar;
 325    }
 326     //attention: la derniere arete de la ligne fermee enveloppe
 327     //           devient en fait la premiere arete de cette ligne
 328     //           dans le chainage des aretes de la frontiere!
 329   }
 330   if( ierr != 0 ) goto ERREUR;
 331
 332   aremin = sqrt( aremin );  //longueur minimale d'une arete des lignes fermees
 333   aremax = sqrt( aremax );  //longueur maximale d'une arete
 334
 335 //debut ajout  9/11/2006  ................................................
 336   // devenu un commentaire aretmx = Min( aretmx, aremax ); //pour homogeneiser
 337
 338   // protection contre une arete max desiree trop grande ou trop petite
 339   if( aretmx > aremax*2.05 ) aretmx = aremax;
 340
 341   // protection contre une arete max desiree trop petite
 342   if( (aremax-aremin) > (aremin+aremax)*0.05 && aretmx < aremin*0.5 )
 343     aretmx =(aremin+aremax*2)/3.0;
 344
 345   if( aretmx < aremin  && aremin > 0 )
 346     aretmx = aremin;
 347
 348   //sauvegarde pour la fonction areteideale_
 349   aretemaxface_ = aretmx;
 350
 351   //aire maximale souhaitee des triangles
 352   airemx = aretmx * aretmx * sqrt(3.0) / 2.0;  //Aire triangle equilateral
 353
 354   for(i=0; i<=nuds; i++ )
 355     mnpxyd[i].z = Min( mnpxyd[i].z, aretmx );
 356   //MESSAGE("Numero du dernier sommet frontalier=" << nuds+1);
 357 //fin  ajout 9/11/2006  .................................................
 358
 359
 360   MESSAGE("Sur  le  bord: arete min=" << aremin << " arete max=" << aremax );
 361   MESSAGE("Triangulation: arete mx=" << aretmx
 362           << " triangle aire mx=" << airemx );
 363
 364   //chainage des aretes frontalieres : la derniere arete frontaliere
 365   mnsoar[ mosoar * noar - mosoar + 5 ] = 0;
 366
 367   //tous les sommets et aretes frontaliers sont numerotes de 1 a nbarfr
 368   //reservation du tableau des numeros des 3 aretes de chaque triangle
 369   //mnartr( moartr, mxartr )
 370   //En nombre: Triangles = Aretes Internes + Aretes Frontalieres - Sommets + 1-Trous
 371   //          3Triangles = 2 Aretes internes + Aretes frontalieres
 372   //       d'ou 3T/2 < AI + AF => T < 3T/2  - Sommets + 1-Trous
 373   //nombre de triangles < 2 ( nombre de sommets - 1 + nombre de trous )
 374   if( mnartr!=NULL ) delete [] mnartr;
 375   mxartr = 2 * ( mxsomm + mxtrou );
 376   mnartr = new Z[moartr*mxartr];
 377   if( mnartr==NULL ) goto ERREUR;
 378
 379   //Ajout des points internes
 380   ns1 = nudslf[ nblf ];
 381   for (i=0; i<nbpti; i++)
 382   {
 383     //les 2 coordonnees du point i de sommet nbs
 384     mnpxyd[ns1].x = uvpti[i].x;
 385     mnpxyd[ns1].y = uvpti[i].y;
 386     mnpxyd[ns1].z = areteideale();//( mnpxyd[ns1], direction );
 387     //le numero i du point interne
 388     mnslig[ns1] = i+1;
 389     ns1++;
 390   }
 391
 392   //nombre de sommets de la frontiere et internes
 393   nbarpi = ns1;
 394
 395   // creation de l'arbre-4 des te (tableau letree)
 396   // ajout dans les te des sommets des lignes et des points internes imposes
 397   // =======================================================================
 398   // premiere estimation de mxtree
 399   mxtree = 2 * mxsomm;
 400
 401  NEWTREE:  //en cas de saturation de l'un des tableaux, on boucle
 402   MESSAGE( "Debut triangulation avec mxsomm=" << mxsomm );
 403   if( mntree != NULL ) delete [] mntree;
 404   nbsomm = nbarpi;
 405   mntree = new Z[motree*(1+mxtree)];
 406   if( mntree==NULL ) goto ERREUR;
 407
 408   //initialisation du tableau letree et ajout dans letree des sommets 1 a nbsomm
 409   teajte( mxsomm, nbsomm, mnpxyd, comxmi, aretmx, mxtree, mntree, ierr );
 410   comxmi[0].z=0;
 411   comxmi[1].z=0;
 412
 413   if( ierr == 51 )
 414   {
 415     //saturation de letree => sa taille est augmentee et relance
 416     mxtree = mxtree * 2;
 417     ierr   = 0;
 418     MESSAGE( "Nouvelle valeur de mxtree=" << mxtree );
 419     goto NEWTREE;
 420   }
 421
 422   deltacpu_( d );
 423   tcpu += d;
 424   MESSAGE( "Temps de l'ajout arbre-4 des Triangles Equilateraux=" << d << " secondes" );
 425   if( ierr != 0 ) goto ERREUR;
 426   //ici le tableau mnpxyd contient les sommets des te et les points frontaliers et internes
 427
 428   // homogeneisation de l'arbre des te a un saut de taille au plus
 429   // prise en compte des tailles d'aretes souhaitees autour des sommets initiaux
 430   // ===========================================================================
 431   // reservation de la queue pour parcourir les te de l'arbre
 432   if( mnqueu != NULL ) delete [] mnqueu;
 433   mxqueu = mxtree;
 434   mnqueu = new Z[mxqueu];
 435   if( mnqueu==NULL) goto ERREUR;
 436
 437   tehote( nutysu, nbarpi, mxsomm, nbsomm, mnpxyd,
 438            comxmi, aretmx,
 439            mntree, mxqueu, mnqueu,
 440            ierr );
 441
 442   deltacpu_( d );
 443   tcpu += d;
 444   MESSAGE("Temps de l'adaptation et l'homogeneisation de l'arbre-4 des TE="
 445        << d << " secondes");
 446   if( ierr != 0 )
 447   {
 448     //destruction du tableau auxiliaire et de l'arbre
 449     if( ierr == 51 )
 450     {
 451       //letree sature
 452       mxtree = mxtree * 2;
 453       MESSAGE( "Redemarrage avec la valeur de mxtree=" << mxtree );
 454       ierr = 0;
 455       goto NEWTREE;
 456     }
 457     else
 458       goto ERREUR;
 459   }
 460
 461   // trianguler les triangles equilateraux feuilles a partir de leurs 3 sommets
 462   // et des points de la frontiere, des points internes imposes interieurs
 463   // ==========================================================================
 464   tetrte( comxmi, aretmx, nbarpi, mxsomm, mnpxyd,
 465            mxqueu, mnqueu, mntree, mosoar, mxsoar, n1soar, mnsoar,
 466            moartr, mxartr, n1artr, mnartr, mnarst,
 467            ierr );
 468
 469   // destruction de la queue et de l'arbre devenus inutiles
 470   delete [] mnqueu;  mnqueu=NULL;
 471   delete [] mntree;  mntree=NULL;
 472
 473   //Temps calcul
 474   deltacpu_( d );
 475   tcpu += d;
 476   MESSAGE( "Temps de la triangulation des TE=" << d << " secondes" );
 477
 478   // ierr =0 si pas d'erreur
 479   //      =1 si le tableau mnsoar est sature
 480   //      =2 si le tableau mnartr est sature
 481   //      =3 si aucun des triangles ne contient l'un des points internes
 482   //      =5 si saturation de la queue de parcours de l'arbre des te
 483   if( ierr != 0 ) goto ERREUR;
 484
 485   //qualites de la triangulation actuelle
 486   qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
 487                 nbt, quamoy, quamin );
 488
 489   // boucle sur les aretes internes (non sur une ligne de la frontiere)
 490   // avec echange des 2 diagonales afin de rendre la triangulation delaunay
 491   // ======================================================================
 492   // formation du chainage 6 des aretes internes a echanger eventuellement
 493   aisoar( mosoar, mxsoar, mnsoar, na );
 494   tedela( mnpxyd, mnarst,
 495            mosoar, mxsoar, n1soar, mnsoar, na,
 496            moartr, mxartr, n1artr, mnartr, n );
 497
 498   MESSAGE( "Nombre d'echanges des diagonales de 2 triangles=" << n );
 499   deltacpu_( d );
 500   tcpu += d;
 501   MESSAGE("Temps de la triangulation Delaunay par echange des diagonales="
 502        << d << " secondes");
 503
 504   //qualites de la triangulation actuelle
 505   qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
 506                 nbt, quamoy, quamin );
 507
 508   // detection des aretes frontalieres initiales perdues
 509   // triangulation frontale pour les restaurer
 510   // ===================================================
 511   mxarcf = mxsomm/5;
 512   if( mn1arcf != NULL ) delete [] mn1arcf;
 513   if( mnarcf  != NULL ) delete [] mnarcf;
 514   if( mnarcf1 != NULL ) delete [] mnarcf1;
 515   if( mnarcf2 != NULL ) delete [] mnarcf2;
 516   mn1arcf = new Z[1+mxarcf];
 517   if( mn1arcf == NULL ) goto ERREUR;
 518   mnarcf  = new Z[3*mxarcf];
 519   if( mnarcf == NULL ) goto ERREUR;
 520   mnarcf1 = new Z[mxarcf];
 521   if( mnarcf1 == NULL ) goto ERREUR;
 522   mnarcf2 = new Z[mxarcf];
 523   if( mnarcf2 == NULL ) goto ERREUR;
 524
 525   terefr( nbarpi, mnpxyd,
 526            mosoar, mxsoar, n1soar, mnsoar,
 527            moartr, mxartr, n1artr, mnartr, mnarst,
 528            mxarcf, mn1arcf, mnarcf, mnarcf1, mnarcf2,
 529            n, ierr );
 530
 531   MESSAGE( "Restauration de " << n << " aretes perdues de la frontiere  ierr=" << ierr );
 532   deltacpu_( d );
 533   tcpu += d;
 534   MESSAGE("Temps de la recuperation des aretes perdues de la frontiere="
 535        << d << " secondes");
 536
 537   if( ierr != 0 ) goto ERREUR;
 538
 539   //qualites de la triangulation actuelle
 540   qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
 541                 nbt, quamoy, quamin );
 542
 543   // fin de la triangulation avec respect des aretes initiales frontalieres
 544
 545   // suppression des triangles externes a la surface
 546   // ===============================================
 547   // recherche du dernier triangle utilise
 548   mn = mxartr * moartr;
 549   for ( ndtri0=mxartr; ndtri0<=1; ndtri0-- )
 550   {
 551     mn -= moartr;
 552     if( mnartr[mn] != 0 ) break;
 553   }
 554
 555   if( mntrsu != NULL ) delete [] mntrsu;
 556   mntrsu = new Z[ndtri0];
 557   if( mntrsu == NULL ) goto ERREUR;
 558
 559   if( mnlftr != NULL ) delete [] mnlftr;
 560   mnlftr = new Z[nblf];
 561   if( mnlftr == NULL ) goto ERREUR;
 562
 563   for (n=0; n<nblf; n++)  //numero de la ligne fermee de 1 a nblf
 564     mnlftr[n] = n+1;
 565
 566   tesuex( nblf,   mnlftr,
 567            ndtri0, nbsomm, mnpxyd, mnslig,
 568            mosoar, mxsoar, mnsoar,
 569            moartr, mxartr, n1artr, mnartr, mnarst,
 570            nbt, mntrsu, ierr );
 571
 572   delete [] mnlftr; mnlftr=NULL;
 573   delete [] mntrsu; mntrsu=NULL;
 574
 575   deltacpu_( d );
 576   tcpu += d;
 577   MESSAGE( "Temps de la suppression des triangles externes=" << d << "ierr=" << ierr );
 578   if( ierr != 0 ) goto ERREUR;
 579
 580   //qualites de la triangulation actuelle
 581   qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
 582                 nbt, quamoy, quamin );
 583
 584   // amelioration de la qualite de la triangulation par
 585   // barycentrage des sommets internes a la triangulation
 586   // suppression des aretes trop longues ou trop courtes
 587   // modification de la topologie des groupes de triangles
 588   // mise en delaunay de la triangulation
 589   // =====================================================
 590   mnarcf3 = new Z[mxarcf];
 591   if( mnarcf3 == NULL )
 592   {
 593     MESSAGE ( "aptrte: MC saturee mnarcf3=" << mnarcf3 );
 594     goto ERREUR;
 595   }
 596   teamqt( nutysu,  aretmx,  airemx,
 597            mnarst,  mosoar,  mxsoar, n1soar, mnsoar,
 598            moartr,  mxartr,  n1artr, mnartr,
 599            mxarcf,  mnarcf2, mnarcf3,
 600            mn1arcf, mnarcf,  mnarcf1,
 601            nbarpi,  nbsomm, mxsomm, mnpxyd, mnslig,
 602            ierr );
 603   if( mnarcf3 != NULL ) {delete [] mnarcf3; mnarcf3=NULL;}
 604   if( mn1arcf != NULL ) {delete [] mn1arcf; mn1arcf=NULL;}
 605   if( mnarcf  != NULL ) {delete [] mnarcf;  mnarcf =NULL;}
 606   if( mnarcf1 != NULL ) {delete [] mnarcf1; mnarcf1=NULL;}
 607   if( mnarcf2 != NULL ) {delete [] mnarcf2; mnarcf2=NULL;}
 608
 609   deltacpu_( d );
 610   tcpu += d;
 611   MESSAGE( "Temps de l'amelioration de la qualite de la triangulation=" << d );
 612   if( ierr == -13 ) ierr=0; //6/10/2006 arret de l'amelioration apres boucle infinie dans caetoi
 613   if( ierr !=   0 ) goto ERREUR;
 614
 615   //qualites de la triangulation finale
 616   qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
 617                 nbt, quamoy, quamin );
 618
 619   // renumerotation des sommets internes: mnarst(i)=numero final du sommet
 620   // ===================================
 621   for (i=0; i<=nbsomm; i++)
 622     mnarst[i] = 0;
 623
 624   for (nt=1; nt<=mxartr; nt++)
 625   {
 626     if( mnartr[nt*moartr-moartr] != 0 )
 627     {
 628       //le numero des 3 sommets du triangle nt
 629       nusotr( nt, mosoar, mnsoar, moartr, mnartr, nosotr );
 630       //les 3 sommets du triangle sont actifs
 631       mnarst[ nosotr[0] ] = 1;
 632       mnarst[ nosotr[1] ] = 1;
 633       mnarst[ nosotr[2] ] = 1;
 634     }
 635   }
 636   nbst = 0;
 637   for (i=1; i<=nbsomm; i++)
 638   {
 639     if( mnarst[i] >0 )
 640       mnarst[i] = ++nbst;
 641   }
 642
 643   // generation du tableau uvst de la surface triangulee
 644   // ---------------------------------------------------
 645   if( uvst != NULL ) delete [] uvst;
 646   uvst = new R2[nbst];
 647   if( uvst == NULL ) goto ERREUR;
 648
 649   nbst=-1;
 650   for (i=0; i<nbsomm; i++ )
 651   {
 652     if( mnarst[i+1]>0 )
 653     {
 654       nbst++;
 655       uvst[nbst].x = mnpxyd[i].x;
 656       uvst[nbst].y = mnpxyd[i].y;
 657
 658       //si le sommet est un point ou appartient a une ligne
 659       //ses coordonnees initiales sont restaurees
 660       n = mnslig[i];
 661       if( n > 0 )
 662       {
 663         if( n >= 1000000 )
 664         {
 665           //sommet d'une ligne
 666           //retour aux coordonnees initiales dans uvslf
 667           l = n / 1000000;
 668           n = n - 1000000 * l + nudslf[l-1] - 1;
 669           uvst[nbst].x = uvslf[n].x;
 670           uvst[nbst].y = uvslf[n].y;
 671         }
 672         else
 673         {
 674           //point utilisateur n interne impose
 675           //retour aux coordonnees initiales dans uvpti
 676           uvst[nbst].x = uvpti[n-1].x;
 677           uvst[nbst].y = uvpti[n-1].y;
 678         }
 679       }
 680     }
 681   }
 682   nbst++;
 683
 684   // generation du tableau 'nsef' de la surface triangulee
 685   // -----------------------------------------------------
 686   // boucle sur les triangles occupes (internes et externes)
 687   if( nust != NULL ) delete [] nust;
 688   nust = new Z[nbsttria*nbt];
 689   if( nust == NULL ) goto ERREUR;
 690   nbt = 0;
 691   for (i=1; i<=mxartr; i++)
 692   {
 693     //le triangle i de mnartr
 694     if( mnartr[i*moartr-moartr] != 0 )
 695     {
 696       //le triangle i est interne => nosotr numero de ses 3 sommets
 697       nusotr( i, mosoar, mnsoar, moartr, mnartr,  nosotr );
 698       nust[nbt++] = mnarst[ nosotr[0] ];
 699       nust[nbt++] = mnarst[ nosotr[1] ];
 700       nust[nbt++] = mnarst[ nosotr[2] ];
 701       nust[nbt++] = 0;
 702     }
 703   }
 704   nbt /= nbsttria;  //le nombre final de triangles de la surface
 705   MESSAGE( "APTRTE: Fin de la triangulation plane avec "<<nbst<<" sommets et "
 706            << nbt << " triangles" );
 707   deltacpu_( d );
 708   tcpu += d;
 709   MESSAGE( "APTRTE: Temps total de la triangulation plane=" << tcpu << " secondes" );
 710
 711   // destruction des tableaux auxiliaires
 712   // ------------------------------------
 713  NETTOYAGE:
 714   if( mnarst != NULL ) delete [] mnarst;
 715   if( mnartr != NULL ) delete [] mnartr;
 716   if( mnslig != NULL ) delete [] mnslig;
 717   if( mnsoar != NULL ) delete [] mnsoar;
 718   if( mnpxyd != NULL ) delete [] mnpxyd;
 719   if( mntree != NULL ) delete [] mntree;
 720   if( mnqueu != NULL ) delete [] mnqueu;
 721   if( mntrsu != NULL ) delete [] mntrsu;
 722   if( mnlftr != NULL ) delete [] mnlftr;
 723   if( mn1arcf != NULL ) delete [] mn1arcf;
 724   if( mnarcf  != NULL ) delete [] mnarcf;
 725   if( mnarcf1 != NULL ) delete [] mnarcf1;
 726   if( mnarcf2 != NULL ) delete [] mnarcf2;
 727   if( mnarcf3 != NULL ) delete [] mnarcf3;
 728   return;
 729
 730  ERREUR:
 731   if( ierr == 51 || ierr == 52 )
 732   {
 733     //saturation des sommets => redepart avec 2 fois plus de sommets
 734     mxsomm = 2 * mxsomm;
 735     ierr   = 0;
 736     goto NEWDEPART;
 737   }
 738   else
 739   {
 740     MESSAGE( "APTRTE: Triangulation NON REALISEE  avec erreur=" << ierr );
 741     if( ierr == 0 ) ierr=1;
 742     goto NETTOYAGE;
 743   }
 744 }
 745 void
 746 #ifdef WIN32
 747               __stdcall
 748 #endif
 749  qualitetrte( R3 *mnpxyd,
 750                    Z & mosoar, Z & mxsoar, Z *mnsoar,
 751                    Z & moartr, Z & mxartr, Z *mnartr,
 752                    Z & nbtria, R & quamoy, R & quamin )
 753 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 754 // but :    calculer la qualite moyenne et minimale de la triangulation
 755 // -----    actuelle definie par les tableaux mnsoar et mnartr
 756 // entrees:
 757 // --------
 758 // mnpxyd : tableau des coordonnees 2d des points
 759 //          par point : x  y  distance_souhaitee
 760 // mosoar : nombre maximal d'entiers par arete et
 761 //          indice dans mnsoar de l'arete suivante dans le hachage
 762 // mxsoar : nombre maximal d'aretes stockables dans le tableau mnsoar
 763 //          attention: mxsoar>3*mxsomm obligatoire!
 764 // mnsoar : numero des 2 sommets , no ligne, 2 triangles de l'arete,
 765 //          chainage des aretes frontalieres, chainage du hachage des aretes
 766 //          hachage des aretes = mnsoar(1)+mnsoar(2)*2
 767 //          avec mxsoar>=3*mxsomm
 768 //          une arete i de mnsoar est vide <=> mnsoar(1,i)=0 et
 769 //          mnsoar(2,arete vide)=l'arete vide qui precede
 770 //          mnsoar(3,arete vide)=l'arete vide qui suit
 771 // moartr : nombre maximal d'entiers par arete du tableau mnartr
 772 // mxartr : nombre maximal de triangles declarables
 773 // mnartr : les 3 aretes des triangles +-arete1, +-arete2, +-arete3
 774 //          arete1 = 0 si triangle vide => arete2 = triangle vide suivant
 775 // sorties:
 776 // --------
 777 // nbtria : nombre de triangles internes au domaine
 778 // quamoy : qualite moyenne  des triangles actuels
 779 // quamin : qualite minimale des triangles actuels
 780 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 781 {
 782   R  d, aire, qualite;
 783   Z  nosotr[3], mn, nbtrianeg, nt, ntqmin;
 784
 785   aire   = 0;
 786   quamoy = 0;
 787   quamin = 2.0;
 788   nbtria = 0;
 789   nbtrianeg = 0;
 790   ntqmin = 0;
 791
 792   mn = -moartr;
 793   for ( nt=1; nt<=mxartr; nt++ )
 794   {
 795     mn += moartr;
 796     if( mnartr[mn]!=0 )
 797     {
 798       //un triangle occupe de plus
 799       nbtria++;
 800
 801       //le numero des 3 sommets du triangle nt
 802       nusotr( nt, mosoar, mnsoar, moartr, mnartr,  nosotr );
 803
 804       //la qualite du triangle ns1 ns2 ns3
 805       qutr2d( mnpxyd[nosotr[0]-1], mnpxyd[nosotr[1]-1], mnpxyd[nosotr[2]-1],
 806                qualite );
 807
 808       //la qualite moyenne
 809       quamoy += qualite;
 810
 811       //la qualite minimale
 812       if( qualite < quamin )
 813       {
 814          quamin = qualite;
 815          ntqmin = nt;
 816       }
 817
 818       //aire signee du triangle nt
 819       d = surtd2( mnpxyd[nosotr[0]-1], mnpxyd[nosotr[1]-1], mnpxyd[nosotr[2]-1] );
 820       if( d<0 )
 821       {
 822         //un triangle d'aire negative de plus
 823         nbtrianeg++;
 824         MESSAGE("ATTENTION: le triangle " << nt << " de sommets:"
 825              << nosotr[0] << " " << nosotr[1] << " " << nosotr[2]
 826              << " a une aire " << d <<"<=0");
 827       }
 828
 829       //aire des triangles actuels
 830       aire += Abs(d);
 831     }
 832   }
 833
 834   //les affichages
 835   quamoy /= nbtria;
 836   MESSAGE("Qualite moyenne=" << quamoy
 837        << "  Qualite minimale=" << quamin
 838        << " des " << nbtria << " triangles de surface plane totale="
 839        << aire);
 840
 841   if( quamin<0.3 )
 842   {
 843     //le numero des 3 sommets du triangle ntqmin de qualite minimale
 844     nusotr(ntqmin, mosoar, mnsoar, moartr, mnartr,  nosotr );
 845     MESSAGE("Triangle de qualite minimale "<<quamin<<" de sommets:"
 846             <<nosotr[0]<<" "<<nosotr[1]<<" "<<nosotr[2]<<" ");
 847     for (int i=0;i<3;i++)
 848       MESSAGE("Sommet "<<nosotr[i]<<": x="<< mnpxyd[nosotr[i]-1].x
 849               <<" y="<< mnpxyd[nosotr[i]-1].y);
 850   }
 851
 852   if( nbtrianeg>0 )
 853     MESSAGE( "ATTENTION: "<< nbtrianeg << " TRIANGLES d'AIRE NEGATIVE" );
 854
 855   MESSAGE(" ");
 856   return;
 857 }