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
 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   teajte( mxsomm, nbsomm, mnpxyd, comxmi, aretmx, mxtree, mntree, ierr );
 419   comxmi[0].z=0;
 420   comxmi[1].z=0;
 421
 422   if( ierr == 51 )
 423   {
 424     //saturation de letree => sa taille est augmentee et relance
 425     mxtree = mxtree * 2;
 426     ierr   = 0;
 427     MESSAGE( "Nouvelle valeur de mxtree=" << mxtree );
 428     goto NEWTREE;
 429   }
 430
 431   deltacpu_( d );
 432   tcpu += d;
 433   MESSAGE( "Temps de l'ajout arbre-4 des Triangles Equilateraux=" << d << " secondes" );
 434   if( ierr != 0 ) goto ERREUR;
 435   //ici le tableau mnpxyd contient les sommets des te et les points frontaliers et internes
 436
 437   // homogeneisation de l'arbre des te a un saut de taille au plus
 438   // prise en compte des tailles d'aretes souhaitees autour des sommets initiaux
 439   // ===========================================================================
 440   // reservation de la queue pour parcourir les te de l'arbre
 441   if( mnqueu != NULL ) delete [] mnqueu;
 442   mxqueu = mxtree;
 443   mnqueu = new Z[mxqueu];
 444   if( mnqueu==NULL) goto ERREUR;
 445
 446   tehote( nutysu, nbarpi, mxsomm, nbsomm, mnpxyd,
 447            comxmi, aretmx,
 448            mntree, mxqueu, mnqueu,
 449            ierr );
 450
 451   deltacpu_( d );
 452   tcpu += d;
 453   MESSAGE("Temps de l'adaptation et l'homogeneisation de l'arbre-4 des TE="
 454        << d << " secondes");
 455   if( ierr != 0 )
 456   {
 457     //destruction du tableau auxiliaire et de l'arbre
 458     if( ierr == 51 )
 459     {
 460       //letree sature
 461       mxtree = mxtree * 2;
 462       MESSAGE( "Redemarrage avec la valeur de mxtree=" << mxtree );
 463       ierr = 0;
 464       goto NEWTREE;
 465     }
 466     else
 467       goto ERREUR;
 468   }
 469
 470   // trianguler les triangles equilateraux feuilles a partir de leurs 3 sommets
 471   // et des points de la frontiere, des points internes imposes interieurs
 472   // ==========================================================================
 473   tetrte( comxmi, aretmx, nbarpi, mxsomm, mnpxyd,
 474            mxqueu, mnqueu, mntree, mosoar, mxsoar, n1soar, mnsoar,
 475            moartr, mxartr, n1artr, mnartr, mnarst,
 476            ierr );
 477
 478   // destruction de la queue et de l'arbre devenus inutiles
 479   delete [] mnqueu;  mnqueu=NULL;
 480   delete [] mntree;  mntree=NULL;
 481
 482   //Temps calcul
 483   deltacpu_( d );
 484   tcpu += d;
 485   MESSAGE( "Temps de la triangulation des TE=" << d << " secondes" );
 486
 487   // ierr =0 si pas d'erreur
 488   //      =1 si le tableau mnsoar est sature
 489   //      =2 si le tableau mnartr est sature
 490   //      =3 si aucun des triangles ne contient l'un des points internes
 491   //      =5 si saturation de la queue de parcours de l'arbre des te
 492   if( ierr != 0 ) goto ERREUR;
 493
 494   //qualites de la triangulation actuelle
 495   qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
 496                 nbt, quamoy, quamin );
 497
 498   // boucle sur les aretes internes (non sur une ligne de la frontiere)
 499   // avec echange des 2 diagonales afin de rendre la triangulation delaunay
 500   // ======================================================================
 501   // formation du chainage 6 des aretes internes a echanger eventuellement
 502   aisoar( mosoar, mxsoar, mnsoar, na );
 503   tedela( mnpxyd, mnarst,
 504            mosoar, mxsoar, n1soar, mnsoar, na,
 505            moartr, mxartr, n1artr, mnartr, n );
 506
 507   MESSAGE( "Nombre d'echanges des diagonales de 2 triangles=" << n );
 508   deltacpu_( d );
 509   tcpu += d;
 510   MESSAGE("Temps de la triangulation Delaunay par echange des diagonales="
 511        << d << " secondes");
 512
 513   //qualites de la triangulation actuelle
 514   qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
 515                 nbt, quamoy, quamin );
 516
 517   // detection des aretes frontalieres initiales perdues
 518   // triangulation frontale pour les restaurer
 519   // ===================================================
 520   mxarcf = mxsomm/5;
 521   if( mn1arcf != NULL ) delete [] mn1arcf;
 522   if( mnarcf  != NULL ) delete [] mnarcf;
 523   if( mnarcf1 != NULL ) delete [] mnarcf1;
 524   if( mnarcf2 != NULL ) delete [] mnarcf2;
 525   mn1arcf = new Z[1+mxarcf];
 526   if( mn1arcf == NULL ) goto ERREUR;
 527   mnarcf  = new Z[3*mxarcf];
 528   if( mnarcf == NULL ) goto ERREUR;
 529   mnarcf1 = new Z[mxarcf];
 530   if( mnarcf1 == NULL ) goto ERREUR;
 531   mnarcf2 = new Z[mxarcf];
 532   if( mnarcf2 == NULL ) goto ERREUR;
 533
 534   terefr( nbarpi, mnpxyd,
 535            mosoar, mxsoar, n1soar, mnsoar,
 536            moartr, mxartr, n1artr, mnartr, mnarst,
 537            mxarcf, mn1arcf, mnarcf, mnarcf1, mnarcf2,
 538            n, ierr );
 539
 540   MESSAGE( "Restauration de " << n << " aretes perdues de la frontiere  ierr=" << ierr );
 541   deltacpu_( d );
 542   tcpu += d;
 543   MESSAGE("Temps de la recuperation des aretes perdues de la frontiere="
 544        << d << " secondes");
 545
 546   if( ierr != 0 ) goto ERREUR;
 547
 548   //qualites de la triangulation actuelle
 549   qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
 550                 nbt, quamoy, quamin );
 551
 552   // fin de la triangulation avec respect des aretes initiales frontalieres
 553
 554   // suppression des triangles externes a la surface
 555   // ===============================================
 556   // recherche du dernier triangle utilise
 557   mn = mxartr * moartr;
 558   for ( ndtri0=mxartr; ndtri0<=1; ndtri0-- )
 559   {
 560     mn -= moartr;
 561     if( mnartr[mn] != 0 ) break;
 562   }
 563
 564   if( mntrsu != NULL ) delete [] mntrsu;
 565   mntrsu = new Z[ndtri0];
 566   if( mntrsu == NULL ) goto ERREUR;
 567
 568   if( mnlftr != NULL ) delete [] mnlftr;
 569   mnlftr = new Z[nblf];
 570   if( mnlftr == NULL ) goto ERREUR;
 571
 572   for (n=0; n<nblf; n++)  //numero de la ligne fermee de 1 a nblf
 573     mnlftr[n] = n+1;
 574
 575   tesuex( nblf,   mnlftr,
 576            ndtri0, nbsomm, mnpxyd, mnslig,
 577            mosoar, mxsoar, mnsoar,
 578            moartr, mxartr, n1artr, mnartr, mnarst,
 579            nbt, mntrsu, ierr );
 580
 581   delete [] mnlftr; mnlftr=NULL;
 582   delete [] mntrsu; mntrsu=NULL;
 583
 584   deltacpu_( d );
 585   tcpu += d;
 586   MESSAGE( "Temps de la suppression des triangles externes=" << d << "ierr=" << ierr );
 587   if( ierr != 0 ) goto ERREUR;
 588
 589   //qualites de la triangulation actuelle
 590   qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
 591                 nbt, quamoy, quamin );
 592
 593   // amelioration de la qualite de la triangulation par
 594   // barycentrage des sommets internes a la triangulation
 595   // suppression des aretes trop longues ou trop courtes
 596   // modification de la topologie des groupes de triangles
 597   // mise en delaunay de la triangulation
 598   // =====================================================
 599   mnarcf3 = new Z[mxarcf];
 600   if( mnarcf3 == NULL )
 601   {
 602     MESSAGE ( "aptrte: MC saturee mnarcf3=" << mnarcf3 );
 603     goto ERREUR;
 604   }
 605   teamqt( nutysu,  aretmx,  airemx,
 606            mnarst,  mosoar,  mxsoar, n1soar, mnsoar,
 607            moartr,  mxartr,  n1artr, mnartr,
 608            mxarcf,  mnarcf2, mnarcf3,
 609            mn1arcf, mnarcf,  mnarcf1,
 610            nbarpi,  nbsomm, mxsomm, mnpxyd, mnslig,
 611            ierr );
 612   if( mnarcf3 != NULL ) {delete [] mnarcf3; mnarcf3=NULL;}
 613   if( mn1arcf != NULL ) {delete [] mn1arcf; mn1arcf=NULL;}
 614   if( mnarcf  != NULL ) {delete [] mnarcf;  mnarcf =NULL;}
 615   if( mnarcf1 != NULL ) {delete [] mnarcf1; mnarcf1=NULL;}
 616   if( mnarcf2 != NULL ) {delete [] mnarcf2; mnarcf2=NULL;}
 617
 618   deltacpu_( d );
 619   tcpu += d;
 620   MESSAGE( "Temps de l'amelioration de la qualite de la triangulation=" << d );
 621   if( ierr == -13 ) ierr=0; //6/10/2006 arret de l'amelioration apres boucle infinie dans caetoi
 622   if( ierr !=   0 ) goto ERREUR;
 623
 624   //qualites de la triangulation finale
 625   qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
 626                 nbt, quamoy, quamin );
 627
 628   // renumerotation des sommets internes: mnarst(i)=numero final du sommet
 629   // ===================================
 630   for (i=0; i<=nbsomm; i++)
 631     mnarst[i] = 0;
 632
 633   for (nt=1; nt<=mxartr; nt++)
 634   {
 635     if( mnartr[nt*moartr-moartr] != 0 )
 636     {
 637       //le numero des 3 sommets du triangle nt
 638       nusotr( nt, mosoar, mnsoar, moartr, mnartr, nosotr );
 639       //les 3 sommets du triangle sont actifs
 640       mnarst[ nosotr[0] ] = 1;
 641       mnarst[ nosotr[1] ] = 1;
 642       mnarst[ nosotr[2] ] = 1;
 643     }
 644   }
 645   nbst = 0;
 646   for (i=1; i<=nbsomm; i++)
 647   {
 648     if( mnarst[i] >0 )
 649       mnarst[i] = ++nbst;
 650   }
 651
 652   // generation du tableau uvst de la surface triangulee
 653   // ---------------------------------------------------
 654   if( uvst != NULL ) delete [] uvst;
 655   uvst = new R2[nbst];
 656   if( uvst == NULL ) goto ERREUR;
 657
 658   nbst=-1;
 659   for (i=0; i<nbsomm; i++ )
 660   {
 661     if( mnarst[i+1]>0 )
 662     {
 663       nbst++;
 664       uvst[nbst].x = mnpxyd[i].x;
 665       uvst[nbst].y = mnpxyd[i].y;
 666
 667       //si le sommet est un point ou appartient a une ligne
 668       //ses coordonnees initiales sont restaurees
 669       n = mnslig[i];
 670       if( n > 0 )
 671       {
 672         if( n >= 1000000 )
 673         {
 674           //sommet d'une ligne
 675           //retour aux coordonnees initiales dans uvslf
 676           l = n / 1000000;
 677           n = n - 1000000 * l + nudslf[l-1] - 1;
 678           uvst[nbst].x = uvslf[n].x;
 679           uvst[nbst].y = uvslf[n].y;
 680         }
 681         else
 682         {
 683           //point utilisateur n interne impose
 684           //retour aux coordonnees initiales dans uvpti
 685           uvst[nbst].x = uvpti[n-1].x;
 686           uvst[nbst].y = uvpti[n-1].y;
 687         }
 688       }
 689     }
 690   }
 691   nbst++;
 692
 693   // generation du tableau 'nsef' de la surface triangulee
 694   // -----------------------------------------------------
 695   // boucle sur les triangles occupes (internes et externes)
 696   if( nust != NULL ) delete [] nust;
 697   nust = new Z[nbsttria*nbt];
 698   if( nust == NULL ) goto ERREUR;
 699   nbt = 0;
 700   for (i=1; i<=mxartr; i++)
 701   {
 702     //le triangle i de mnartr
 703     if( mnartr[i*moartr-moartr] != 0 )
 704     {
 705       //le triangle i est interne => nosotr numero de ses 3 sommets
 706       nusotr( i, mosoar, mnsoar, moartr, mnartr,  nosotr );
 707       nust[nbt++] = mnarst[ nosotr[0] ];
 708       nust[nbt++] = mnarst[ nosotr[1] ];
 709       nust[nbt++] = mnarst[ nosotr[2] ];
 710       nust[nbt++] = 0;
 711     }
 712   }
 713   nbt /= nbsttria;  //le nombre final de triangles de la surface
 714   MESSAGE( "APTRTE: Fin de la triangulation plane avec "<<nbst<<" sommets et "
 715            << nbt << " triangles" );
 716   deltacpu_( d );
 717   tcpu += d;
 718   MESSAGE( "APTRTE: Temps total de la triangulation plane=" << tcpu << " secondes" );
 719
 720   // destruction des tableaux auxiliaires
 721   // ------------------------------------
 722  NETTOYAGE:
 723   if( mnarst != NULL ) delete [] mnarst;
 724   if( mnartr != NULL ) delete [] mnartr;
 725   if( mnslig != NULL ) delete [] mnslig;
 726   if( mnsoar != NULL ) delete [] mnsoar;
 727   if( mnpxyd != NULL ) delete [] mnpxyd;
 728   if( mntree != NULL ) delete [] mntree;
 729   if( mnqueu != NULL ) delete [] mnqueu;
 730   if( mntrsu != NULL ) delete [] mntrsu;
 731   if( mnlftr != NULL ) delete [] mnlftr;
 732   if( mn1arcf != NULL ) delete [] mn1arcf;
 733   if( mnarcf  != NULL ) delete [] mnarcf;
 734   if( mnarcf1 != NULL ) delete [] mnarcf1;
 735   if( mnarcf2 != NULL ) delete [] mnarcf2;
 736   if( mnarcf3 != NULL ) delete [] mnarcf3;
 737   return;
 738
 739  ERREUR:
 740   if( ierr == 51 || ierr == 52 )
 741   {
 742     //saturation des sommets => redepart avec 2 fois plus de sommets
 743     mxsomm = 2 * mxsomm;
 744     ierr   = 0;
 745     goto NEWDEPART;
 746   }
 747   else
 748   {
 749     MESSAGE( "APTRTE: Triangulation NON REALISEE  avec erreur=" << ierr );
 750     if( ierr == 0 ) ierr=1;
 751     goto NETTOYAGE;
 752   }
 753 }
 754 void
 755 #ifdef WIN32
 756 #ifdef F2C_BUILD
 757 #else
 758               __stdcall
 759 #endif
 760 #endif
 761  qualitetrte( R3 *mnpxyd,
 762                    Z & mosoar, Z & mxsoar, Z *mnsoar,
 763                    Z & moartr, Z & mxartr, Z *mnartr,
 764                    Z & nbtria, R & quamoy, R & quamin )
 765 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 766 // but :    calculer la qualite moyenne et minimale de la triangulation
 767 // -----    actuelle definie par les tableaux mnsoar et mnartr
 768 // entrees:
 769 // --------
 770 // mnpxyd : tableau des coordonnees 2d des points
 771 //          par point : x  y  distance_souhaitee
 772 // mosoar : nombre maximal d'entiers par arete et
 773 //          indice dans mnsoar de l'arete suivante dans le hachage
 774 // mxsoar : nombre maximal d'aretes stockables dans le tableau mnsoar
 775 //          attention: mxsoar>3*mxsomm obligatoire!
 776 // mnsoar : numero des 2 sommets , no ligne, 2 triangles de l'arete,
 777 //          chainage des aretes frontalieres, chainage du hachage des aretes
 778 //          hachage des aretes = mnsoar(1)+mnsoar(2)*2
 779 //          avec mxsoar>=3*mxsomm
 780 //          une arete i de mnsoar est vide <=> mnsoar(1,i)=0 et
 781 //          mnsoar(2,arete vide)=l'arete vide qui precede
 782 //          mnsoar(3,arete vide)=l'arete vide qui suit
 783 // moartr : nombre maximal d'entiers par arete du tableau mnartr
 784 // mxartr : nombre maximal de triangles declarables
 785 // mnartr : les 3 aretes des triangles +-arete1, +-arete2, +-arete3
 786 //          arete1 = 0 si triangle vide => arete2 = triangle vide suivant
 787 // sorties:
 788 // --------
 789 // nbtria : nombre de triangles internes au domaine
 790 // quamoy : qualite moyenne  des triangles actuels
 791 // quamin : qualite minimale des triangles actuels
 792 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 793 {
 794   R  d, aire, qualite;
 795   Z  nosotr[3], mn, nbtrianeg, nt, ntqmin;
 796
 797   aire   = 0;
 798   quamoy = 0;
 799   quamin = 2.0;
 800   nbtria = 0;
 801   nbtrianeg = 0;
 802   ntqmin = 0;
 803
 804   mn = -moartr;
 805   for ( nt=1; nt<=mxartr; nt++ )
 806   {
 807     mn += moartr;
 808     if( mnartr[mn]!=0 )
 809     {
 810       //un triangle occupe de plus
 811       nbtria++;
 812
 813       //le numero des 3 sommets du triangle nt
 814       nusotr( nt, mosoar, mnsoar, moartr, mnartr,  nosotr );
 815
 816       //la qualite du triangle ns1 ns2 ns3
 817       qutr2d( mnpxyd[nosotr[0]-1], mnpxyd[nosotr[1]-1], mnpxyd[nosotr[2]-1],
 818                qualite );
 819
 820       //la qualite moyenne
 821       quamoy += qualite;
 822
 823       //la qualite minimale
 824       if( qualite < quamin )
 825       {
 826          quamin = qualite;
 827          ntqmin = nt;
 828       }
 829
 830       //aire signee du triangle nt
 831       d = surtd2( mnpxyd[nosotr[0]-1], mnpxyd[nosotr[1]-1], mnpxyd[nosotr[2]-1] );
 832       if( d<0 )
 833       {
 834         //un triangle d'aire negative de plus
 835         nbtrianeg++;
 836         MESSAGE("ATTENTION: le triangle " << nt << " de sommets:"
 837              << nosotr[0] << " " << nosotr[1] << " " << nosotr[2]
 838              << " a une aire " << d <<"<=0");
 839       }
 840
 841       //aire des triangles actuels
 842       aire += Abs(d);
 843     }
 844   }
 845
 846   //les affichages
 847   quamoy /= nbtria;
 848   MESSAGE("Qualite moyenne=" << quamoy
 849        << "  Qualite minimale=" << quamin
 850        << " des " << nbtria << " triangles de surface plane totale="
 851        << aire);
 852
 853   if( quamin<0.3 )
 854   {
 855     //le numero des 3 sommets du triangle ntqmin de qualite minimale
 856     nusotr(ntqmin, mosoar, mnsoar, moartr, mnartr,  nosotr );
 857     MESSAGE("Triangle de qualite minimale "<<quamin<<" de sommets:"
 858             <<nosotr[0]<<" "<<nosotr[1]<<" "<<nosotr[2]<<" ");
 859     for (int i=0;i<3;i++)
 860       MESSAGE("Sommet "<<nosotr[i]<<": x="<< mnpxyd[nosotr[i]-1].x
 861               <<" y="<< mnpxyd[nosotr[i]-1].y);
 862   }
 863
 864   if( nbtrianeg>0 )
 865     MESSAGE( "ATTENTION: "<< nbtrianeg << " TRIANGLES d'AIRE NEGATIVE" );
 866
 867   MESSAGE(" ");
 868   return;
 869 }