src/MEFISTO2/aptrte.cxx

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