src/MEFISTO2/aptrte.cxx

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