src/MEFISTO2/aptrte.cxx

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