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