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