src/MEFISTO2/aptrte.cxx

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