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