src/MEFISTO2/aptrte.cxx

   1 //  MEFISTO2: a library to compute 2D triangulation from segmented boundaries
   2 //
   3 //
   4 //  Copyright (C) 2006  Laboratoire J.-L. Lions UPMC Paris
   5 //
   6 //  This library is free software; you can redistribute it and/or
   7 //  modify it under the terms of the GNU Lesser General Public
   8 //  License as published by the Free Software Foundation; either
   9 //  version 2.1 of the License.
  10 //
  11 //  This library is distributed in the hope that it will be useful,
  12 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  13 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 //  Lesser General Public License for more details.
  15 //
  16 //  You should have received a copy of the GNU Lesser General Public
  17 //  License along with this library; if not, write to the Free Software
  18 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  19 //
  20 //  See http://www.ann.jussieu.fr/~perronne or email Perronnet@ann.jussieu.fr
  21 //
  22 //  File   : aptrte.cxx   le C++ de l'appel du trianguleur plan
  23 //  Module : SMESH
  24 //  Author : Alain PERRONNET
  25 //  Date   : 13 novembre 2006
  26
  27 #include "Rn.h"
  28 #include "aptrte.h"
  29 #include "utilities.h"
  30
  31 using namespace std;
  32
  33 extern "C"
  34 {
  35   R aretemaxface_;
  36   MEFISTO2D_EXPORT
  37     R
  38   #ifdef WIN32
  39   #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 #define _MAXFACT 4.05
 350   if( aretmx > aremin*_MAXFACT ) aretmx = aremin*_MAXFACT;
 351
 352   // protection contre une arete max desiree trop petite
 353 //   if( (aremax-aremin) > (aremin+aremax)*0.05 && aretmx < aremin*0.5 )
 354 //     aretmx =(aremin+aremax*2)/3.0;
 355
 356 //   if( aretmx < aremin  && aremin > 0 )
 357 //     aretmx = aremin;
 358
 359   //sauvegarde pour la fonction areteideale_
 360   aretemaxface_ = aretmx;
 361
 362   //aire maximale souhaitee des triangles
 363   airemx = aretmx * aretmx * sqrt(3.0) / 2.0;  //Aire triangle equilateral
 364
 365   for(i=0; i<=nuds; i++ )
 366     mnpxyd[i].z = Min( mnpxyd[i].z, aretmx );
 367   //MESSAGE("Numero du dernier sommet frontalier=" << nuds+1);
 368 //fin  ajout 9/11/2006  .................................................
 369
 370
 371   MESSAGE("Sur  le  bord: arete min=" << aremin << " arete max=" << aremax );
 372   MESSAGE("Triangulation: arete mx=" << aretmx
 373           << " triangle aire mx=" << airemx );
 374
 375   //chainage des aretes frontalieres : la derniere arete frontaliere
 376   mnsoar[ mosoar * noar - mosoar + 5 ] = 0;
 377
 378   //tous les sommets et aretes frontaliers sont numerotes de 1 a nbarfr
 379   //reservation du tableau des numeros des 3 aretes de chaque triangle
 380   //mnartr( moartr, mxartr )
 381   //En nombre: Triangles = Aretes Internes + Aretes Frontalieres - Sommets + 1-Trous
 382   //          3Triangles = 2 Aretes internes + Aretes frontalieres
 383   //       d'ou 3T/2 < AI + AF => T < 3T/2  - Sommets + 1-Trous
 384   //nombre de triangles < 2 ( nombre de sommets - 1 + nombre de trous )
 385   if( mnartr!=NULL ) delete [] mnartr;
 386   mxartr = 2 * ( mxsomm + mxtrou );
 387   mnartr = new Z[moartr*mxartr];
 388   if( mnartr==NULL ) goto ERREUR;
 389
 390   //Ajout des points internes
 391   ns1 = nudslf[ nblf ];
 392   for (i=0; i<nbpti; i++)
 393   {
 394     //les 2 coordonnees du point i de sommet nbs
 395     mnpxyd[ns1].x = uvpti[i].x;
 396     mnpxyd[ns1].y = uvpti[i].y;
 397     mnpxyd[ns1].z = areteideale();//( mnpxyd[ns1], direction );
 398     //le numero i du point interne
 399     mnslig[ns1] = i+1;
 400     ns1++;
 401   }
 402
 403   //nombre de sommets de la frontiere et internes
 404   nbarpi = ns1;
 405
 406   // creation de l'arbre-4 des te (tableau letree)
 407   // ajout dans les te des sommets des lignes et des points internes imposes
 408   // =======================================================================
 409   // premiere estimation de mxtree
 410   mxtree = 2 * mxsomm;
 411
 412  NEWTREE:  //en cas de saturation de l'un des tableaux, on boucle
 413   MESSAGE( "Debut triangulation avec mxsomm=" << mxsomm );
 414   if( mntree != NULL ) delete [] mntree;
 415   nbsomm = nbarpi;
 416   mntree = new Z[motree*(1+mxtree)];
 417   if( mntree==NULL ) goto ERREUR;
 418
 419   //initialisation du tableau letree et ajout dans letree des sommets 1 a nbsomm
 420   comxmi[0].x = comxmi[1].x = uvslf[0].x;
 421   comxmi[0].y = comxmi[1].y = uvslf[0].y;
 422   teajte( mxsomm, nbsomm, mnpxyd, comxmi, aretmx, mxtree, mntree, ierr );
 423   comxmi[0].z=0;
 424   comxmi[1].z=0;
 425
 426   if( ierr == 51 )
 427   {
 428     //saturation de letree => sa taille est augmentee et relance
 429     mxtree = mxtree * 2;
 430     ierr   = 0;
 431     MESSAGE( "Nouvelle valeur de mxtree=" << mxtree );
 432     goto NEWTREE;
 433   }
 434
 435   deltacpu_( d );
 436   tcpu += d;
 437   MESSAGE( "Temps de l'ajout arbre-4 des Triangles Equilateraux=" << d << " secondes" );
 438   if( ierr != 0 ) goto ERREUR;
 439   //ici le tableau mnpxyd contient les sommets des te et les points frontaliers et internes
 440
 441   // homogeneisation de l'arbre des te a un saut de taille au plus
 442   // prise en compte des tailles d'aretes souhaitees autour des sommets initiaux
 443   // ===========================================================================
 444   // reservation de la queue pour parcourir les te de l'arbre
 445   if( mnqueu != NULL ) delete [] mnqueu;
 446   mxqueu = mxtree;
 447   mnqueu = new Z[mxqueu];
 448   if( mnqueu==NULL) goto ERREUR;
 449
 450   tehote( nutysu, nbarpi, mxsomm, nbsomm, mnpxyd,
 451            comxmi, aretmx,
 452            mntree, mxqueu, mnqueu,
 453            ierr );
 454
 455   deltacpu_( d );
 456   tcpu += d;
 457   MESSAGE("Temps de l'adaptation et l'homogeneisation de l'arbre-4 des TE="
 458        << d << " secondes");
 459   if( ierr != 0 )
 460   {
 461     //destruction du tableau auxiliaire et de l'arbre
 462     if( ierr == 51 )
 463     {
 464       //letree sature
 465       mxtree = mxtree * 2;
 466       MESSAGE( "Redemarrage avec la valeur de mxtree=" << mxtree );
 467       ierr = 0;
 468       goto NEWTREE;
 469     }
 470     else
 471       goto ERREUR;
 472   }
 473
 474   // trianguler les triangles equilateraux feuilles a partir de leurs 3 sommets
 475   // et des points de la frontiere, des points internes imposes interieurs
 476   // ==========================================================================
 477   tetrte( comxmi, aretmx, nbarpi, mxsomm, mnpxyd,
 478            mxqueu, mnqueu, mntree, mosoar, mxsoar, n1soar, mnsoar,
 479            moartr, mxartr, n1artr, mnartr, mnarst,
 480            ierr );
 481
 482   // destruction de la queue et de l'arbre devenus inutiles
 483   delete [] mnqueu;  mnqueu=NULL;
 484   delete [] mntree;  mntree=NULL;
 485
 486   //Temps calcul
 487   deltacpu_( d );
 488   tcpu += d;
 489   MESSAGE( "Temps de la triangulation des TE=" << d << " secondes" );
 490
 491   // ierr =0 si pas d'erreur
 492   //      =1 si le tableau mnsoar est sature
 493   //      =2 si le tableau mnartr est sature
 494   //      =3 si aucun des triangles ne contient l'un des points internes
 495   //      =5 si saturation de la queue de parcours de l'arbre des te
 496   if( ierr != 0 ) goto ERREUR;
 497
 498   //qualites de la triangulation actuelle
 499   qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
 500                 nbt, quamoy, quamin );
 501
 502   // boucle sur les aretes internes (non sur une ligne de la frontiere)
 503   // avec echange des 2 diagonales afin de rendre la triangulation delaunay
 504   // ======================================================================
 505   // formation du chainage 6 des aretes internes a echanger eventuellement
 506   aisoar( mosoar, mxsoar, mnsoar, na );
 507   tedela( mnpxyd, mnarst,
 508            mosoar, mxsoar, n1soar, mnsoar, na,
 509            moartr, mxartr, n1artr, mnartr, n );
 510
 511   MESSAGE( "Nombre d'echanges des diagonales de 2 triangles=" << n );
 512   deltacpu_( d );
 513   tcpu += d;
 514   MESSAGE("Temps de la triangulation Delaunay par echange des diagonales="
 515        << d << " secondes");
 516
 517   //qualites de la triangulation actuelle
 518   qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
 519                 nbt, quamoy, quamin );
 520
 521   // detection des aretes frontalieres initiales perdues
 522   // triangulation frontale pour les restaurer
 523   // ===================================================
 524   mxarcf = mxsomm/5;
 525   if( mn1arcf != NULL ) delete [] mn1arcf;
 526   if( mnarcf  != NULL ) delete [] mnarcf;
 527   if( mnarcf1 != NULL ) delete [] mnarcf1;
 528   if( mnarcf2 != NULL ) delete [] mnarcf2;
 529   mn1arcf = new Z[1+mxarcf];
 530   if( mn1arcf == NULL ) goto ERREUR;
 531   mnarcf  = new Z[3*mxarcf];
 532   if( mnarcf == NULL ) goto ERREUR;
 533   mnarcf1 = new Z[mxarcf];
 534   if( mnarcf1 == NULL ) goto ERREUR;
 535   mnarcf2 = new Z[mxarcf];
 536   if( mnarcf2 == NULL ) goto ERREUR;
 537
 538   terefr( nbarpi, mnpxyd,
 539            mosoar, mxsoar, n1soar, mnsoar,
 540            moartr, mxartr, n1artr, mnartr, mnarst,
 541            mxarcf, mn1arcf, mnarcf, mnarcf1, mnarcf2,
 542            n, ierr );
 543
 544   MESSAGE( "Restauration de " << n << " aretes perdues de la frontiere  ierr=" << ierr );
 545   deltacpu_( d );
 546   tcpu += d;
 547   MESSAGE("Temps de la recuperation des aretes perdues de la frontiere="
 548        << d << " secondes");
 549
 550   if( ierr != 0 ) goto ERREUR;
 551
 552   //qualites de la triangulation actuelle
 553   qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
 554                 nbt, quamoy, quamin );
 555
 556   // fin de la triangulation avec respect des aretes initiales frontalieres
 557
 558   // suppression des triangles externes a la surface
 559   // ===============================================
 560   // recherche du dernier triangle utilise
 561   mn = mxartr * moartr;
 562   for ( ndtri0=mxartr; ndtri0<=1; ndtri0-- )
 563   {
 564     mn -= moartr;
 565     if( mnartr[mn] != 0 ) break;
 566   }
 567
 568   if( mntrsu != NULL ) delete [] mntrsu;
 569   mntrsu = new Z[ndtri0];
 570   if( mntrsu == NULL ) goto ERREUR;
 571
 572   if( mnlftr != NULL ) delete [] mnlftr;
 573   mnlftr = new Z[nblf];
 574   if( mnlftr == NULL ) goto ERREUR;
 575
 576   for (n=0; n<nblf; n++)  //numero de la ligne fermee de 1 a nblf
 577     mnlftr[n] = n+1;
 578
 579   tesuex( nblf,   mnlftr,
 580            ndtri0, nbsomm, mnpxyd, mnslig,
 581            mosoar, mxsoar, mnsoar,
 582            moartr, mxartr, n1artr, mnartr, mnarst,
 583            nbt, mntrsu, ierr );
 584
 585   delete [] mnlftr; mnlftr=NULL;
 586   delete [] mntrsu; mntrsu=NULL;
 587
 588   deltacpu_( d );
 589   tcpu += d;
 590   MESSAGE( "Temps de la suppression des triangles externes=" << d << "ierr=" << ierr );
 591   if( ierr != 0 ) goto ERREUR;
 592
 593   //qualites de la triangulation actuelle
 594   qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
 595                 nbt, quamoy, quamin );
 596
 597   // amelioration de la qualite de la triangulation par
 598   // barycentrage des sommets internes a la triangulation
 599   // suppression des aretes trop longues ou trop courtes
 600   // modification de la topologie des groupes de triangles
 601   // mise en delaunay de la triangulation
 602   // =====================================================
 603   mnarcf3 = new Z[mxarcf];
 604   if( mnarcf3 == NULL )
 605   {
 606     MESSAGE ( "aptrte: MC saturee mnarcf3=" << mnarcf3 );
 607     goto ERREUR;
 608   }
 609   teamqt( nutysu,  aretmx,  airemx,
 610            mnarst,  mosoar,  mxsoar, n1soar, mnsoar,
 611            moartr,  mxartr,  n1artr, mnartr,
 612            mxarcf,  mnarcf2, mnarcf3,
 613            mn1arcf, mnarcf,  mnarcf1,
 614            nbarpi,  nbsomm, mxsomm, mnpxyd, mnslig,
 615            ierr );
 616   if( mnarcf3 != NULL ) {delete [] mnarcf3; mnarcf3=NULL;}
 617   if( mn1arcf != NULL ) {delete [] mn1arcf; mn1arcf=NULL;}
 618   if( mnarcf  != NULL ) {delete [] mnarcf;  mnarcf =NULL;}
 619   if( mnarcf1 != NULL ) {delete [] mnarcf1; mnarcf1=NULL;}
 620   if( mnarcf2 != NULL ) {delete [] mnarcf2; mnarcf2=NULL;}
 621
 622   deltacpu_( d );
 623   tcpu += d;
 624   MESSAGE( "Temps de l'amelioration de la qualite de la triangulation=" << d );
 625   if( ierr == -13 ) ierr=0; //6/10/2006 arret de l'amelioration apres boucle infinie dans caetoi
 626   if( ierr !=   0 ) goto ERREUR;
 627
 628   //qualites de la triangulation finale
 629   qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
 630                 nbt, quamoy, quamin );
 631
 632   // renumerotation des sommets internes: mnarst(i)=numero final du sommet
 633   // ===================================
 634   for (i=0; i<=nbsomm; i++)
 635     mnarst[i] = 0;
 636
 637   for (nt=1; nt<=mxartr; nt++)
 638   {
 639     if( mnartr[nt*moartr-moartr] != 0 )
 640     {
 641       //le numero des 3 sommets du triangle nt
 642       nusotr( nt, mosoar, mnsoar, moartr, mnartr, nosotr );
 643       //les 3 sommets du triangle sont actifs
 644       mnarst[ nosotr[0] ] = 1;
 645       mnarst[ nosotr[1] ] = 1;
 646       mnarst[ nosotr[2] ] = 1;
 647     }
 648   }
 649   nbst = 0;
 650   for (i=1; i<=nbsomm; i++)
 651   {
 652     if( mnarst[i] >0 )
 653       mnarst[i] = ++nbst;
 654   }
 655
 656   // generation du tableau uvst de la surface triangulee
 657   // ---------------------------------------------------
 658   if( uvst != NULL ) delete [] uvst;
 659   uvst = new R2[nbst];
 660   if( uvst == NULL ) goto ERREUR;
 661
 662   nbst=-1;
 663   for (i=0; i<nbsomm; i++ )
 664   {
 665     if( mnarst[i+1]>0 )
 666     {
 667       nbst++;
 668       uvst[nbst].x = mnpxyd[i].x;
 669       uvst[nbst].y = mnpxyd[i].y;
 670
 671       //si le sommet est un point ou appartient a une ligne
 672       //ses coordonnees initiales sont restaurees
 673       n = mnslig[i];
 674       if( n > 0 )
 675       {
 676         if( n >= 1000000 )
 677         {
 678           //sommet d'une ligne
 679           //retour aux coordonnees initiales dans uvslf
 680           l = n / 1000000;
 681           n = n - 1000000 * l + nudslf[l-1] - 1;
 682           uvst[nbst].x = uvslf[n].x;
 683           uvst[nbst].y = uvslf[n].y;
 684         }
 685         else
 686         {
 687           //point utilisateur n interne impose
 688           //retour aux coordonnees initiales dans uvpti
 689           uvst[nbst].x = uvpti[n-1].x;
 690           uvst[nbst].y = uvpti[n-1].y;
 691         }
 692       }
 693     }
 694   }
 695   nbst++;
 696
 697   // generation du tableau 'nsef' de la surface triangulee
 698   // -----------------------------------------------------
 699   // boucle sur les triangles occupes (internes et externes)
 700   if( nust != NULL ) delete [] nust;
 701   nust = new Z[nbsttria*nbt];
 702   if( nust == NULL ) goto ERREUR;
 703   nbt = 0;
 704   for (i=1; i<=mxartr; i++)
 705   {
 706     //le triangle i de mnartr
 707     if( mnartr[i*moartr-moartr] != 0 )
 708     {
 709       //le triangle i est interne => nosotr numero de ses 3 sommets
 710       nusotr( i, mosoar, mnsoar, moartr, mnartr,  nosotr );
 711       nust[nbt++] = mnarst[ nosotr[0] ];
 712       nust[nbt++] = mnarst[ nosotr[1] ];
 713       nust[nbt++] = mnarst[ nosotr[2] ];
 714       nust[nbt++] = 0;
 715     }
 716   }
 717   nbt /= nbsttria;  //le nombre final de triangles de la surface
 718   MESSAGE( "APTRTE: Fin de la triangulation plane avec "<<nbst<<" sommets et "
 719            << nbt << " triangles" );
 720   deltacpu_( d );
 721   tcpu += d;
 722   MESSAGE( "APTRTE: Temps total de la triangulation plane=" << tcpu << " secondes" );
 723
 724   // destruction des tableaux auxiliaires
 725   // ------------------------------------
 726  NETTOYAGE:
 727   if( mnarst != NULL ) delete [] mnarst;
 728   if( mnartr != NULL ) delete [] mnartr;
 729   if( mnslig != NULL ) delete [] mnslig;
 730   if( mnsoar != NULL ) delete [] mnsoar;
 731   if( mnpxyd != NULL ) delete [] mnpxyd;
 732   if( mntree != NULL ) delete [] mntree;
 733   if( mnqueu != NULL ) delete [] mnqueu;
 734   if( mntrsu != NULL ) delete [] mntrsu;
 735   if( mnlftr != NULL ) delete [] mnlftr;
 736   if( mn1arcf != NULL ) delete [] mn1arcf;
 737   if( mnarcf  != NULL ) delete [] mnarcf;
 738   if( mnarcf1 != NULL ) delete [] mnarcf1;
 739   if( mnarcf2 != NULL ) delete [] mnarcf2;
 740   if( mnarcf3 != NULL ) delete [] mnarcf3;
 741   return;
 742
 743  ERREUR:
 744   if( ierr == 51 || ierr == 52 )
 745   {
 746     //saturation des sommets => redepart avec 2 fois plus de sommets
 747     mxsomm = 2 * mxsomm;
 748     ierr   = 0;
 749     goto NEWDEPART;
 750   }
 751   else
 752   {
 753     MESSAGE( "APTRTE: Triangulation NON REALISEE  avec erreur=" << ierr );
 754     if( ierr == 0 ) ierr=1;
 755     goto NETTOYAGE;
 756   }
 757 }
 758 void
 759 #ifdef WIN32
 760 #ifdef F2C_BUILD
 761 #else
 762               __stdcall
 763 #endif
 764 #endif
 765  qualitetrte( R3 *mnpxyd,
 766                    Z & mosoar, Z & mxsoar, Z *mnsoar,
 767                    Z & moartr, Z & mxartr, Z *mnartr,
 768                    Z & nbtria, R & quamoy, R & quamin )
 769 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 770 // but :    calculer la qualite moyenne et minimale de la triangulation
 771 // -----    actuelle definie par les tableaux mnsoar et mnartr
 772 // entrees:
 773 // --------
 774 // mnpxyd : tableau des coordonnees 2d des points
 775 //          par point : x  y  distance_souhaitee
 776 // mosoar : nombre maximal d'entiers par arete et
 777 //          indice dans mnsoar de l'arete suivante dans le hachage
 778 // mxsoar : nombre maximal d'aretes stockables dans le tableau mnsoar
 779 //          attention: mxsoar>3*mxsomm obligatoire!
 780 // mnsoar : numero des 2 sommets , no ligne, 2 triangles de l'arete,
 781 //          chainage des aretes frontalieres, chainage du hachage des aretes
 782 //          hachage des aretes = mnsoar(1)+mnsoar(2)*2
 783 //          avec mxsoar>=3*mxsomm
 784 //          une arete i de mnsoar est vide <=> mnsoar(1,i)=0 et
 785 //          mnsoar(2,arete vide)=l'arete vide qui precede
 786 //          mnsoar(3,arete vide)=l'arete vide qui suit
 787 // moartr : nombre maximal d'entiers par arete du tableau mnartr
 788 // mxartr : nombre maximal de triangles declarables
 789 // mnartr : les 3 aretes des triangles +-arete1, +-arete2, +-arete3
 790 //          arete1 = 0 si triangle vide => arete2 = triangle vide suivant
 791 // sorties:
 792 // --------
 793 // nbtria : nombre de triangles internes au domaine
 794 // quamoy : qualite moyenne  des triangles actuels
 795 // quamin : qualite minimale des triangles actuels
 796 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 797 {
 798   R  d, aire, qualite;
 799   Z  nosotr[3], mn, nbtrianeg, nt, ntqmin;
 800
 801   aire   = 0;
 802   quamoy = 0;
 803   quamin = 2.0;
 804   nbtria = 0;
 805   nbtrianeg = 0;
 806   ntqmin = 0;
 807
 808   mn = -moartr;
 809   for ( nt=1; nt<=mxartr; nt++ )
 810   {
 811     mn += moartr;
 812     if( mnartr[mn]!=0 )
 813     {
 814       //un triangle occupe de plus
 815       nbtria++;
 816
 817       //le numero des 3 sommets du triangle nt
 818       nusotr( nt, mosoar, mnsoar, moartr, mnartr,  nosotr );
 819
 820       //la qualite du triangle ns1 ns2 ns3
 821       qutr2d( mnpxyd[nosotr[0]-1], mnpxyd[nosotr[1]-1], mnpxyd[nosotr[2]-1],
 822                qualite );
 823
 824       //la qualite moyenne
 825       quamoy += qualite;
 826
 827       //la qualite minimale
 828       if( qualite < quamin )
 829       {
 830          quamin = qualite;
 831          ntqmin = nt;
 832       }
 833
 834       //aire signee du triangle nt
 835       d = surtd2( mnpxyd[nosotr[0]-1], mnpxyd[nosotr[1]-1], mnpxyd[nosotr[2]-1] );
 836       if( d<0 )
 837       {
 838         //un triangle d'aire negative de plus
 839         nbtrianeg++;
 840         MESSAGE("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   MESSAGE("Qualite moyenne=" << quamoy
 853        << "  Qualite minimale=" << quamin
 854        << " des " << nbtria << " triangles de surface plane totale="
 855        << aire);
 856
 857   if( quamin<0.3 )
 858   {
 859     //le numero des 3 sommets du triangle ntqmin de qualite minimale
 860     nusotr(ntqmin, mosoar, mnsoar, moartr, mnartr,  nosotr );
 861     MESSAGE("Triangle de qualite minimale "<<quamin<<" de sommets:"
 862             <<nosotr[0]<<" "<<nosotr[1]<<" "<<nosotr[2]<<" ");
 863     for (int i=0;i<3;i++)
 864       MESSAGE("Sommet "<<nosotr[i]<<": x="<< mnpxyd[nosotr[i]-1].x
 865               <<" y="<< mnpxyd[nosotr[i]-1].y);
 866   }
 867
 868   if( nbtrianeg>0 )
 869     MESSAGE( "ATTENTION: "<< nbtrianeg << " TRIANGLES d'AIRE NEGATIVE" );
 870
 871   MESSAGE(" ");
 872   return;
 873 }