src/MEFISTO2/aptrte.cxx

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