src/MEFISTO2/aptrte.cxx

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