src/MEFISTO2/aptrte.cxx

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