Back to index

salome-smesh  6.5.0
aptrte.cxx
Go to the documentation of this file.
00001 //  MEFISTO2: a library to compute 2D triangulation from segmented boundaries
00002 //
00003 //
00004 // Copyright (C) 2006-2012  CEA/DEN, EDF R&D, OPEN CASCADE
00005 //
00006 // This library is free software; you can redistribute it and/or
00007 // modify it under the terms of the GNU Lesser General Public
00008 // License as published by the Free Software Foundation; either
00009 // version 2.1 of the License.
00010 //
00011 // This library is distributed in the hope that it will be useful,
00012 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00013 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014 // Lesser General Public License for more details.
00015 //
00016 // You should have received a copy of the GNU Lesser General Public
00017 // License along with this library; if not, write to the Free Software
00018 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00019 //
00020 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00021 //
00022 //  File   : aptrte.cxx   le C++ de l'appel du trianguleur plan
00023 //  Module : SMESH
00024 //  Author : Alain PERRONNET
00025 //  Date   : 13 novembre 2006
00026 
00027 #include "Rn.h"
00028 #include "aptrte.h"
00029 #include "utilities.h"
00030 
00031 using namespace std;
00032 
00033 extern "C"
00034 {
00035   R aretemaxface_;
00036   MEFISTO2D_EXPORT   
00037     R
00038   #ifdef WIN32
00039   #ifdef F2C_BUILD
00040   #else
00041       __stdcall
00042   #endif
00043   #endif
00044       areteideale()//( R3 xyz, R3 direction )
00045   {
00046     return aretemaxface_;
00047   }
00048 }
00049 //calcul de la longueur ideale de l'arete au sommet xyz (z ici inactif)
00050 //dans la direction donnee
00051 //a ajuster pour chaque surface plane et selon l'entier notysu (voir plus bas)
00052 
00053 
00054 static double cpunew, cpuold=0;
00055 
00056 void
00057 #ifdef WIN32
00058 #ifdef F2C_BUILD
00059 #else
00060               __stdcall
00061 #endif
00062 #endif
00063 tempscpu_( double & tempsec )
00064 //Retourne le temps CPU utilise en secondes
00065 {  
00066   tempsec = ( (double) clock() ) / CLOCKS_PER_SEC;
00067   //MESSAGE( "temps cpu=" << tempsec );
00068 }
00069 
00070 
00071 void
00072 #ifdef WIN32
00073 #ifdef F2C_BUILD
00074 #else
00075               __stdcall
00076 #endif
00077 #endif
00078 deltacpu_( R & dtcpu )
00079 //Retourne le temps CPU utilise en secondes depuis le precedent appel
00080 {
00081   tempscpu_( cpunew );
00082   dtcpu  = R( cpunew - cpuold );
00083   cpuold = cpunew;
00084   //MESSAGE( "delta temps cpu=" << dtcpu );
00085   return;
00086 }
00087 
00088 
00089 void  aptrte( Z   nutysu, R      aretmx,
00090               Z   nblf,   Z  *   nudslf,  R2 * uvslf,
00091               Z   nbpti,  R2 *   uvpti,
00092               Z & nbst,   R2 * & uvst,
00093               Z & nbt,    Z  * & nust,
00094               Z & ierr )
00095 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00096 // but : appel de la triangulation par un arbre-4 recouvrant
00097 // ----- de triangles equilateraux
00098 //       le contour du domaine plan est defini par des lignes fermees
00099 //       la premiere ligne etant l'enveloppe de toutes les autres
00100 //       la fonction areteideale(s,d) donne la taille d'arete
00101 //       au point s dans la direction (actuellement inactive) d
00102 //       des lors toute arete issue d'un sommet s devrait avoir une longueur
00103 //       comprise entre 0.65 areteideale_(s,d) et 1.3 areteideale_(s,d)
00104 //
00105 //Attention:
00106 //  Les tableaux uvslf et uvpti sont supposes ne pas avoir de sommets identiques!
00107 //  De meme, un sommet d'une ligne fermee ne peut appartenir a une autre ligne fermee
00108 //
00109 // entrees:
00110 // --------
00111 // nutysu : numero de traitement de areteideale_(s,d) selon le type de surface
00112 //          0 pas d'emploi de la fonction areteideale_() et aretmx est active
00113 //          1 il existe une fonction areteideale_(s,d)
00114 //            dont seules les 2 premieres composantes de uv sont actives
00115 //          ... autres options a definir ...
00116 // aretmx : longueur maximale des aretes de la future triangulation
00117 // nblf   : nombre de lignes fermees de la surface
00118 // nudslf : numero du dernier sommet de chacune des nblf lignes fermees
00119 //          nudslf(0)=0 pour permettre la difference sans test
00120 //          Attention le dernier sommet de chaque ligne est raccorde au premier
00121 //          tous les sommets et les points internes ont des coordonnees
00122 //          UV differentes <=> Pas de point double!
00123 // uvslf  : uv des nudslf(nblf) sommets des lignes fermees
00124 // nbpti  : nombre de points internes futurs sommets de la triangulation
00125 // uvpti  : uv des points internes futurs sommets de la triangulation
00126 //
00127 // sorties:
00128 // --------
00129 // nbst   : nombre de sommets de la triangulation finale
00130 // uvst   : coordonnees uv des nbst sommets de la triangulation
00131 // nbt    : nombre de triangles de la triangulation finale
00132 // nust   : 4 numeros dans uvst des sommets des nbt triangles
00133 //          s1, s2, s3, 0: no dans uvst des 3 sommets et 0 car quadrangle!
00134 // ierr   : 0 si pas d'erreur
00135 //        > 0 sinon
00136 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00137 // auteur : Alain Perronnet  Laboratoire J.-L. LIONS Paris UPMC mars 2006
00138 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00139 {
00140   Z  nbsttria=4; //Attention: 4 sommets stockes par triangle
00141                  //no st1, st2, st3, 0 (non quadrangle)
00142 
00143   R  d, tcpu=0;
00144 //  R3 direction=R3(0,0,0);  //direction pour areteideale() inactive ici!
00145   Z  nbarfr=nudslf[nblf];  //nombre total d'aretes des lignes fermees
00146   Z  mxtrou = Max( 1024, nblf );  //nombre maximal de trous dans la surface
00147 
00148   R3 *mnpxyd=NULL;
00149   Z  *mnsoar=NULL, mosoar=7, mxsoar, n1soar; //le hachage des aretes
00150   Z  *mnartr=NULL, moartr=3, mxartr, n1artr; //le no des 3 aretes des triangles
00151   Z  *mntree=NULL, motree=9, mxtree; //L'arbre 4 de TE et nombre d'entiers par TE
00152   Z  *mnqueu=NULL, mxqueu;
00153   Z  *mn1arcf=NULL;
00154   Z  *mnarcf=NULL, mxarcf;
00155   Z  *mnarcf1=NULL;
00156   Z  *mnarcf2=NULL;
00157   Z  *mnarcf3=NULL;
00158   Z  *mntrsu=NULL;
00159   Z  *mnslig=NULL;
00160   Z  *mnarst=NULL;
00161   Z  *mnlftr=NULL;
00162 
00163   R3 comxmi[2];            //coordonnees UV Min et Maximales
00164   R  aremin, aremax;       //longueur minimale et maximale des aretes
00165   R  airemx;               //aire maximale souhaitee d'un triangle
00166   R  quamoy, quamin;
00167 
00168   Z  noar0, noar, na;
00169   Z  i, l, n, ns, ns0, ns1, ns2, nosotr[3], nt;
00170   Z  mxsomm, nbsomm, nbarpi, nbarli, ndtri0, mn;
00171   Z  moins1=-1;
00172   Z  nuds = 0;
00173 
00174   // initialisation du temps cpu
00175   deltacpu_( d );
00176   ierr = 0;
00177 
00178   // quelques reservations de tableaux pour faire les calculs
00179   // ========================================================
00180   // declaration du tableau des coordonnees des sommets de la frontiere
00181   // puis des sommets internes ajoutes
00182   // majoration empirique du nombre de sommets de la triangulation
00183   i = 4*nbarfr/10;
00184   mxsomm = Max( 20000, 64*nbpti+i*i );
00185   MESSAGE( "APTRTE: Debut de la triangulation plane avec " );
00186   MESSAGE( "nutysu=" << nutysu << "  aretmx=" << aretmx
00187            << "  mxsomm=" << mxsomm );
00188   MESSAGE( nbarfr << " sommets sur la frontiere et " << nbpti << " points internes");
00189 
00190  NEWDEPART:
00191   //mnpxyd( 3, mxsomm ) les coordonnees UV des sommets et la taille d'arete aux sommets
00192   if( mnpxyd!=NULL ) delete [] mnpxyd;
00193   mnpxyd = new R3[mxsomm];
00194   if( mnpxyd==NULL ) goto ERREUR;
00195 
00196   // le tableau mnsoar des aretes des triangles
00197   // 1: sommet 1 dans pxyd,
00198   // 2: sommet 2 dans pxyd,
00199   // 3: numero de 1 a nblf de la ligne qui supporte l'arete
00200   // 4: numero dans mnartr du triangle 1 partageant cette arete,
00201   // 5: numero dans mnartr du triangle 2 partageant cette arete,
00202   // 6: chainage des aretes frontalieres ou internes ou
00203   //    des aretes simples des etoiles de triangles,
00204   // 7: chainage du hachage des aretes
00205   // nombre d'aretes = 3 ( nombre de sommets - 1 + nombre de trous )
00206   // pour le hachage des aretes mxsoar doit etre > 3*mxsomm!
00207   // h(ns1,ns2) = min( ns1, ns2 )
00208   if( mnsoar!=NULL ) delete [] mnsoar;
00209   mxsoar = 3 * ( mxsomm + mxtrou );
00210   mnsoar = new Z[mosoar*mxsoar];
00211   if( mnsoar==NULL ) goto ERREUR;
00212   //initialiser le tableau mnsoar pour le hachage des aretes
00213   insoar( mxsomm, mosoar, mxsoar, n1soar, mnsoar );
00214 
00215   // mnarst( mxsomm ) numero mnsoar d'une arete pour chacun des sommets
00216   if( mnarst!=NULL ) delete [] mnarst;
00217   mnarst = new Z[1+mxsomm];
00218   if( mnarst==NULL ) goto ERREUR;
00219   n = 1+mxsomm;
00220   azeroi( n, mnarst );
00221 
00222   // mnslig( mxsomm ) no de sommet dans sa ligne pour chaque sommet frontalier
00223   //               ou no du point si interne forc'e par l'utilisateur
00224   //               ou  0 si interne cree par le module
00225   if( mnslig!=NULL ) delete [] mnslig;
00226   mnslig = new Z[mxsomm];
00227   if( mnslig==NULL ) goto ERREUR;
00228   azeroi( mxsomm, mnslig );
00229 
00230   // initialisation des aretes frontalieres de la triangulation future
00231   // renumerotation des sommets des aretes des lignes pour la triangulation
00232   // mise a l'echelle des coordonnees des sommets pour obtenir une
00233   // meilleure precision lors des calculs + quelques verifications
00234   // boucle sur les lignes fermees qui forment la frontiere
00235   // ======================================================================
00236   noar = 0;
00237   aremin = 1e100;
00238   aremax = 0;
00239 
00240   for (n=1; n<=nblf; n++)
00241   {
00242     //l'initialisation de la premiere arete de la ligne n dans la triangulation
00243     //-------------------------------------------------------------------------
00244     //le sommet ns0 est le numero de l'origine de la ligne
00245     ns0 = nudslf[n-1];
00246     mnpxyd[ns0].x = uvslf[ns0].x;
00247     mnpxyd[ns0].y = uvslf[ns0].y;
00248     mnpxyd[ns0].z = areteideale();//( mnpxyd[ns0], direction );
00249 //     MESSAGE("Sommet " << ns0 << ": " << mnpxyd[ns0].x
00250 //       << " " << mnpxyd[ns0].y << " longueur arete=" << mnpxyd[ns0].z);
00251 
00252     //carre de la longueur de l'arete 1 de la ligne fermee n
00253     d = pow( uvslf[ns0+1].x - uvslf[ns0].x, 2 ) 
00254       + pow( uvslf[ns0+1].y - uvslf[ns0].y, 2 ) ;
00255     aremin = Min( aremin, d );
00256     aremax = Max( aremax, d );
00257 
00258     //le numero des 2 sommets (ns1,ns2) de la premiere arete de la ligne
00259     //initialisation de la 1-ere arete ns1-ns1+1 de cette ligne fermee n
00260     //le numero des 2 sommets ns1 ns2 de la 1-ere arete
00261     //Attention: les numeros ns debutent a 1 (ils ont >0)
00262     //           les tableaux c++ demarrent a zero!
00263     //           les tableaux fortran demarrent ou l'on veut!
00264     ns0++;
00265     ns1 = ns0;
00266     ns2 = ns1+1;
00267 
00268      //le numero n de la ligne du sommet et son numero ns1 dans la ligne
00269     mnslig[ns0-1] = 1000000 * n + ns1-nudslf[n-1];
00270     fasoar( ns1, ns2, moins1, moins1, n,
00271              mosoar, mxsoar, n1soar, mnsoar, mnarst,
00272              noar0,  ierr );
00273     //pas de test sur ierr car pas de saturation possible a ce niveau
00274 
00275     //le pointeur dans le hachage sur la premiere arete de la ligne fermee n
00276     //mndalf[n] = noar0;
00277 
00278     //la nouvelle arete est la suivante de l'arete definie juste avant
00279     if( noar > 0 )
00280       mnsoar[mosoar * noar - mosoar + 5] = noar0;
00281 
00282     //l'initialisation des aretes suivantes de la ligne dans la triangulation
00283     //-----------------------------------------------------------------------
00284     nbarli = nudslf[n] - nudslf[n-1];  //nombre d'aretes=sommets de la ligne n
00285     for (i=2; i<=nbarli; i++)
00286     {
00287       ns1 = ns2; //le numero de l'arete et le numero du premier sommet de l'arete
00288       if( i < nbarli )
00289         //nbs+1 est le 2-eme sommet de l'arete i de la ligne fermee n
00290         ns2 = ns1+1;
00291       else
00292         //le 2-eme sommet de la derniere arete est le premier sommet de la ligne
00293         ns2 = ns0;
00294 
00295       //l'arete precedente est dotee de sa suivante:celle cree ensuite
00296       //les 2 coordonnees du sommet ns2 de la ligne
00297       ns = ns1 - 1;
00298 //debut ajout  5/10/2006  ................................................
00299       nuds = Max( nuds, ns );   //le numero du dernier sommet traite
00300 //fin   ajout  5/10/2006  ................................................
00301       mnpxyd[ns].x = uvslf[ns].x;
00302       mnpxyd[ns].y = uvslf[ns].y;
00303       mnpxyd[ns].z = areteideale();//( mnpxyd[ns], direction );
00304 //       MESSAGE("Sommet " << ns << ": " << mnpxyd[ns].x
00305 //         << " " << mnpxyd[ns].y << " longueur arete=" << mnpxyd[ns].z);
00306 
00307       //carre de la longueur de l'arete
00308       d = pow( uvslf[ns2-1].x - uvslf[ns1-1].x, 2) 
00309         + pow( uvslf[ns2-1].y - uvslf[ns1-1].y, 2);
00310       aremin = Min( aremin, d );
00311       aremax = Max( aremax, d );
00312 
00313 //debut ajout du 5/10/2006  .............................................
00314       //la longueur de l'arete ns1-ns2
00315       d = sqrt( d );
00316       //longueur arete = Min ( aretmx, aretes incidentes )
00317       mnpxyd[ns   ].z = Min( mnpxyd[ns   ].z, d );
00318       mnpxyd[ns2-1].z = Min( mnpxyd[ns2-1].z, d );
00319 //fin ajout du 5/10/2006  ...............................................
00320 
00321       //le numero n de la ligne du sommet et son numero ns1 dans la ligne
00322       mnslig[ns] = 1000000 * n + ns1-nudslf[n-1];
00323 
00324       //ajout de l'arete dans la liste
00325       fasoar( ns1, ns2, moins1, moins1, n,
00326                mosoar, mxsoar, n1soar, mnsoar,
00327                mnarst, noar, ierr );
00328       //pas de test sur ierr car pas de saturation possible a ce niveau
00329 
00330       //chainage des aretes frontalieres en position 6 du tableau mnsoar
00331       //la nouvelle arete est la suivante de l'arete definie juste avant
00332       mnsoar[ mosoar * noar0 - mosoar + 5 ] = noar;
00333       noar0 = noar;
00334    }
00335     //attention: la derniere arete de la ligne fermee enveloppe
00336     //           devient en fait la premiere arete de cette ligne
00337     //           dans le chainage des aretes de la frontiere!
00338   }
00339   if( ierr != 0 ) goto ERREUR;
00340 
00341   aremin = sqrt( aremin );  //longueur minimale d'une arete des lignes fermees
00342   aremax = sqrt( aremax );  //longueur maximale d'une arete
00343 
00344 //debut ajout  9/11/2006  ................................................
00345   // devenu un commentaire aretmx = Min( aretmx, aremax ); //pour homogeneiser
00346 
00347   // protection contre une arete max desiree trop grande ou trop petite
00348   if( aretmx > aremax*2.05 ) aretmx = aremax;
00349 
00350   // protection contre une arete max desiree trop petite
00351   if( (aremax-aremin) > (aremin+aremax)*0.05 && aretmx < aremin*0.5 )
00352     aretmx =(aremin+aremax*2)/3.0;
00353 
00354   if( aretmx < aremin  && aremin > 0 )
00355     aretmx = aremin;
00356 
00357   //sauvegarde pour la fonction areteideale_
00358   aretemaxface_ = aretmx;
00359 
00360   //aire maximale souhaitee des triangles
00361   airemx = aretmx * aretmx * sqrt(3.0) / 2.0;  //Aire triangle equilateral
00362 
00363   for(i=0; i<=nuds; i++ )
00364     mnpxyd[i].z = Min( mnpxyd[i].z, aretmx );
00365   //MESSAGE("Numero du dernier sommet frontalier=" << nuds+1);
00366 //fin  ajout 9/11/2006  .................................................
00367 
00368 
00369   MESSAGE("Sur  le  bord: arete min=" << aremin << " arete max=" << aremax );
00370   MESSAGE("Triangulation: arete mx=" << aretmx
00371           << " triangle aire mx=" << airemx );
00372 
00373   //chainage des aretes frontalieres : la derniere arete frontaliere
00374   mnsoar[ mosoar * noar - mosoar + 5 ] = 0;
00375 
00376   //tous les sommets et aretes frontaliers sont numerotes de 1 a nbarfr
00377   //reservation du tableau des numeros des 3 aretes de chaque triangle
00378   //mnartr( moartr, mxartr )
00379   //En nombre: Triangles = Aretes Internes + Aretes Frontalieres - Sommets + 1-Trous
00380   //          3Triangles = 2 Aretes internes + Aretes frontalieres
00381   //       d'ou 3T/2 < AI + AF => T < 3T/2  - Sommets + 1-Trous
00382   //nombre de triangles < 2 ( nombre de sommets - 1 + nombre de trous )
00383   if( mnartr!=NULL ) delete [] mnartr;
00384   mxartr = 2 * ( mxsomm + mxtrou );
00385   mnartr = new Z[moartr*mxartr];
00386   if( mnartr==NULL ) goto ERREUR;
00387 
00388   //Ajout des points internes
00389   ns1 = nudslf[ nblf ];
00390   for (i=0; i<nbpti; i++)
00391   {
00392     //les 2 coordonnees du point i de sommet nbs
00393     mnpxyd[ns1].x = uvpti[i].x;
00394     mnpxyd[ns1].y = uvpti[i].y;
00395     mnpxyd[ns1].z = areteideale();//( mnpxyd[ns1], direction );
00396     //le numero i du point interne
00397     mnslig[ns1] = i+1;
00398     ns1++;
00399   }
00400 
00401   //nombre de sommets de la frontiere et internes
00402   nbarpi = ns1;
00403 
00404   // creation de l'arbre-4 des te (tableau letree)
00405   // ajout dans les te des sommets des lignes et des points internes imposes
00406   // =======================================================================
00407   // premiere estimation de mxtree
00408   mxtree = 2 * mxsomm;
00409 
00410  NEWTREE:  //en cas de saturation de l'un des tableaux, on boucle
00411   MESSAGE( "Debut triangulation avec mxsomm=" << mxsomm );
00412   if( mntree != NULL ) delete [] mntree;
00413   nbsomm = nbarpi;
00414   mntree = new Z[motree*(1+mxtree)];
00415   if( mntree==NULL ) goto ERREUR;
00416 
00417   //initialisation du tableau letree et ajout dans letree des sommets 1 a nbsomm
00418   comxmi[0].x = comxmi[1].x = uvslf[0].x;
00419   comxmi[0].y = comxmi[1].y = uvslf[0].y;
00420   teajte( mxsomm, nbsomm, mnpxyd, comxmi, aretmx, mxtree, mntree, ierr );
00421   comxmi[0].z=0;
00422   comxmi[1].z=0;
00423 
00424   if( ierr == 51 )
00425   {
00426     //saturation de letree => sa taille est augmentee et relance
00427     mxtree = mxtree * 2;
00428     ierr   = 0;
00429     MESSAGE( "Nouvelle valeur de mxtree=" << mxtree );
00430     goto NEWTREE;
00431   }
00432 
00433   deltacpu_( d );
00434   tcpu += d;
00435   MESSAGE( "Temps de l'ajout arbre-4 des Triangles Equilateraux=" << d << " secondes" );
00436   if( ierr != 0 ) goto ERREUR;
00437   //ici le tableau mnpxyd contient les sommets des te et les points frontaliers et internes
00438 
00439   // homogeneisation de l'arbre des te a un saut de taille au plus
00440   // prise en compte des tailles d'aretes souhaitees autour des sommets initiaux
00441   // ===========================================================================
00442   // reservation de la queue pour parcourir les te de l'arbre
00443   if( mnqueu != NULL ) delete [] mnqueu;
00444   mxqueu = mxtree;
00445   mnqueu = new Z[mxqueu];
00446   if( mnqueu==NULL) goto ERREUR;
00447 
00448   tehote( nutysu, nbarpi, mxsomm, nbsomm, mnpxyd,
00449            comxmi, aretmx,
00450            mntree, mxqueu, mnqueu,
00451            ierr );
00452 
00453   deltacpu_( d );
00454   tcpu += d;
00455   MESSAGE("Temps de l'adaptation et l'homogeneisation de l'arbre-4 des TE="
00456        << d << " secondes");
00457   if( ierr != 0 )
00458   {
00459     //destruction du tableau auxiliaire et de l'arbre
00460     if( ierr == 51 )
00461     {
00462       //letree sature
00463       mxtree = mxtree * 2;
00464       MESSAGE( "Redemarrage avec la valeur de mxtree=" << mxtree );
00465       ierr = 0;
00466       goto NEWTREE;
00467     }
00468     else
00469       goto ERREUR;
00470   }
00471 
00472   // trianguler les triangles equilateraux feuilles a partir de leurs 3 sommets
00473   // et des points de la frontiere, des points internes imposes interieurs
00474   // ==========================================================================
00475   tetrte( comxmi, aretmx, nbarpi, mxsomm, mnpxyd,
00476            mxqueu, mnqueu, mntree, mosoar, mxsoar, n1soar, mnsoar,
00477            moartr, mxartr, n1artr, mnartr, mnarst,
00478            ierr );
00479 
00480   // destruction de la queue et de l'arbre devenus inutiles
00481   delete [] mnqueu;  mnqueu=NULL;
00482   delete [] mntree;  mntree=NULL;
00483 
00484   //Temps calcul
00485   deltacpu_( d );
00486   tcpu += d;
00487   MESSAGE( "Temps de la triangulation des TE=" << d << " secondes" );
00488 
00489   // ierr =0 si pas d'erreur
00490   //      =1 si le tableau mnsoar est sature
00491   //      =2 si le tableau mnartr est sature
00492   //      =3 si aucun des triangles ne contient l'un des points internes
00493   //      =5 si saturation de la queue de parcours de l'arbre des te
00494   if( ierr != 0 ) goto ERREUR;
00495 
00496   //qualites de la triangulation actuelle
00497   qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
00498                 nbt, quamoy, quamin );
00499 
00500   // boucle sur les aretes internes (non sur une ligne de la frontiere)
00501   // avec echange des 2 diagonales afin de rendre la triangulation delaunay
00502   // ======================================================================
00503   // formation du chainage 6 des aretes internes a echanger eventuellement
00504   aisoar( mosoar, mxsoar, mnsoar, na );
00505   tedela( mnpxyd, mnarst,
00506            mosoar, mxsoar, n1soar, mnsoar, na,
00507            moartr, mxartr, n1artr, mnartr, n );
00508 
00509   MESSAGE( "Nombre d'echanges des diagonales de 2 triangles=" << n );
00510   deltacpu_( d );
00511   tcpu += d;
00512   MESSAGE("Temps de la triangulation Delaunay par echange des diagonales="
00513        << d << " secondes");
00514 
00515   //qualites de la triangulation actuelle
00516   qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
00517                 nbt, quamoy, quamin );
00518 
00519   // detection des aretes frontalieres initiales perdues
00520   // triangulation frontale pour les restaurer
00521   // ===================================================
00522   mxarcf = mxsomm/5;
00523   if( mn1arcf != NULL ) delete [] mn1arcf;
00524   if( mnarcf  != NULL ) delete [] mnarcf;
00525   if( mnarcf1 != NULL ) delete [] mnarcf1;
00526   if( mnarcf2 != NULL ) delete [] mnarcf2;
00527   mn1arcf = new Z[1+mxarcf];
00528   if( mn1arcf == NULL ) goto ERREUR;
00529   mnarcf  = new Z[3*mxarcf];
00530   if( mnarcf == NULL ) goto ERREUR;
00531   mnarcf1 = new Z[mxarcf];
00532   if( mnarcf1 == NULL ) goto ERREUR;
00533   mnarcf2 = new Z[mxarcf];
00534   if( mnarcf2 == NULL ) goto ERREUR;
00535 
00536   terefr( nbarpi, mnpxyd,
00537            mosoar, mxsoar, n1soar, mnsoar,
00538            moartr, mxartr, n1artr, mnartr, mnarst,
00539            mxarcf, mn1arcf, mnarcf, mnarcf1, mnarcf2,
00540            n, ierr );
00541 
00542   MESSAGE( "Restauration de " << n << " aretes perdues de la frontiere  ierr=" << ierr );
00543   deltacpu_( d );
00544   tcpu += d;
00545   MESSAGE("Temps de la recuperation des aretes perdues de la frontiere="
00546        << d << " secondes");
00547 
00548   if( ierr != 0 ) goto ERREUR;
00549 
00550   //qualites de la triangulation actuelle
00551   qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
00552                 nbt, quamoy, quamin );
00553 
00554   // fin de la triangulation avec respect des aretes initiales frontalieres
00555 
00556   // suppression des triangles externes a la surface
00557   // ===============================================
00558   // recherche du dernier triangle utilise
00559   mn = mxartr * moartr;
00560   for ( ndtri0=mxartr; ndtri0<=1; ndtri0-- )
00561   {
00562     mn -= moartr;
00563     if( mnartr[mn] != 0 ) break;
00564   }
00565 
00566   if( mntrsu != NULL ) delete [] mntrsu;
00567   mntrsu = new Z[ndtri0];
00568   if( mntrsu == NULL ) goto ERREUR;
00569 
00570   if( mnlftr != NULL ) delete [] mnlftr;
00571   mnlftr = new Z[nblf];
00572   if( mnlftr == NULL ) goto ERREUR;
00573 
00574   for (n=0; n<nblf; n++)  //numero de la ligne fermee de 1 a nblf
00575     mnlftr[n] = n+1;
00576 
00577   tesuex( nblf,   mnlftr,
00578            ndtri0, nbsomm, mnpxyd, mnslig,
00579            mosoar, mxsoar, mnsoar,
00580            moartr, mxartr, n1artr, mnartr, mnarst,
00581            nbt, mntrsu, ierr );
00582 
00583   delete [] mnlftr; mnlftr=NULL;
00584   delete [] mntrsu; mntrsu=NULL;
00585 
00586   deltacpu_( d );
00587   tcpu += d;
00588   MESSAGE( "Temps de la suppression des triangles externes=" << d << "ierr=" << ierr );
00589   if( ierr != 0 ) goto ERREUR;
00590 
00591   //qualites de la triangulation actuelle
00592   qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
00593                 nbt, quamoy, quamin );
00594 
00595   // amelioration de la qualite de la triangulation par
00596   // barycentrage des sommets internes a la triangulation
00597   // suppression des aretes trop longues ou trop courtes
00598   // modification de la topologie des groupes de triangles
00599   // mise en delaunay de la triangulation
00600   // =====================================================
00601   mnarcf3 = new Z[mxarcf];
00602   if( mnarcf3 == NULL )
00603   {
00604     MESSAGE ( "aptrte: MC saturee mnarcf3=" << mnarcf3 );
00605     goto ERREUR;
00606   }
00607   teamqt( nutysu,  aretmx,  airemx,
00608            mnarst,  mosoar,  mxsoar, n1soar, mnsoar,
00609            moartr,  mxartr,  n1artr, mnartr,
00610            mxarcf,  mnarcf2, mnarcf3,
00611            mn1arcf, mnarcf,  mnarcf1,
00612            nbarpi,  nbsomm, mxsomm, mnpxyd, mnslig,
00613            ierr );
00614   if( mnarcf3 != NULL ) {delete [] mnarcf3; mnarcf3=NULL;}
00615   if( mn1arcf != NULL ) {delete [] mn1arcf; mn1arcf=NULL;}
00616   if( mnarcf  != NULL ) {delete [] mnarcf;  mnarcf =NULL;}
00617   if( mnarcf1 != NULL ) {delete [] mnarcf1; mnarcf1=NULL;}
00618   if( mnarcf2 != NULL ) {delete [] mnarcf2; mnarcf2=NULL;}
00619 
00620   deltacpu_( d );
00621   tcpu += d;
00622   MESSAGE( "Temps de l'amelioration de la qualite de la triangulation=" << d );
00623   if( ierr == -13 ) ierr=0; //6/10/2006 arret de l'amelioration apres boucle infinie dans caetoi
00624   if( ierr !=   0 ) goto ERREUR;
00625 
00626   //qualites de la triangulation finale
00627   qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
00628                 nbt, quamoy, quamin );
00629 
00630   // renumerotation des sommets internes: mnarst(i)=numero final du sommet
00631   // ===================================
00632   for (i=0; i<=nbsomm; i++)
00633     mnarst[i] = 0;
00634 
00635   for (nt=1; nt<=mxartr; nt++)
00636   {
00637     if( mnartr[nt*moartr-moartr] != 0 )
00638     {
00639       //le numero des 3 sommets du triangle nt
00640       nusotr( nt, mosoar, mnsoar, moartr, mnartr, nosotr );
00641       //les 3 sommets du triangle sont actifs
00642       mnarst[ nosotr[0] ] = 1;
00643       mnarst[ nosotr[1] ] = 1;
00644       mnarst[ nosotr[2] ] = 1;
00645     }
00646   }
00647   nbst = 0;
00648   for (i=1; i<=nbsomm; i++)
00649   {
00650     if( mnarst[i] >0 )
00651       mnarst[i] = ++nbst;
00652   }
00653 
00654   // generation du tableau uvst de la surface triangulee
00655   // ---------------------------------------------------
00656   if( uvst != NULL ) delete [] uvst;
00657   uvst = new R2[nbst];
00658   if( uvst == NULL ) goto ERREUR;
00659 
00660   nbst=-1;
00661   for (i=0; i<nbsomm; i++ )
00662   {
00663     if( mnarst[i+1]>0 )
00664     {
00665       nbst++;
00666       uvst[nbst].x = mnpxyd[i].x;
00667       uvst[nbst].y = mnpxyd[i].y;
00668 
00669       //si le sommet est un point ou appartient a une ligne
00670       //ses coordonnees initiales sont restaurees
00671       n = mnslig[i];
00672       if( n > 0 )
00673       {
00674         if( n >= 1000000 )
00675         {
00676           //sommet d'une ligne
00677           //retour aux coordonnees initiales dans uvslf
00678           l = n / 1000000;
00679           n = n - 1000000 * l + nudslf[l-1] - 1;
00680           uvst[nbst].x = uvslf[n].x;
00681           uvst[nbst].y = uvslf[n].y;
00682         }
00683         else
00684         {
00685           //point utilisateur n interne impose
00686           //retour aux coordonnees initiales dans uvpti
00687           uvst[nbst].x = uvpti[n-1].x;
00688           uvst[nbst].y = uvpti[n-1].y;
00689         }
00690       }
00691     }
00692   }
00693   nbst++;
00694 
00695   // generation du tableau 'nsef' de la surface triangulee
00696   // -----------------------------------------------------
00697   // boucle sur les triangles occupes (internes et externes)
00698   if( nust != NULL ) delete [] nust;
00699   nust = new Z[nbsttria*nbt];
00700   if( nust == NULL ) goto ERREUR;
00701   nbt = 0;
00702   for (i=1; i<=mxartr; i++)
00703   {
00704     //le triangle i de mnartr
00705     if( mnartr[i*moartr-moartr] != 0 )
00706     {
00707       //le triangle i est interne => nosotr numero de ses 3 sommets
00708       nusotr( i, mosoar, mnsoar, moartr, mnartr,  nosotr );
00709       nust[nbt++] = mnarst[ nosotr[0] ];
00710       nust[nbt++] = mnarst[ nosotr[1] ];
00711       nust[nbt++] = mnarst[ nosotr[2] ];
00712       nust[nbt++] = 0;
00713     }
00714   }
00715   nbt /= nbsttria;  //le nombre final de triangles de la surface
00716   MESSAGE( "APTRTE: Fin de la triangulation plane avec "<<nbst<<" sommets et "
00717            << nbt << " triangles" );
00718   deltacpu_( d );
00719   tcpu += d;
00720   MESSAGE( "APTRTE: Temps total de la triangulation plane=" << tcpu << " secondes" );
00721 
00722   // destruction des tableaux auxiliaires
00723   // ------------------------------------
00724  NETTOYAGE:
00725   if( mnarst != NULL ) delete [] mnarst;
00726   if( mnartr != NULL ) delete [] mnartr;
00727   if( mnslig != NULL ) delete [] mnslig;
00728   if( mnsoar != NULL ) delete [] mnsoar;
00729   if( mnpxyd != NULL ) delete [] mnpxyd;
00730   if( mntree != NULL ) delete [] mntree;
00731   if( mnqueu != NULL ) delete [] mnqueu;
00732   if( mntrsu != NULL ) delete [] mntrsu;
00733   if( mnlftr != NULL ) delete [] mnlftr;
00734   if( mn1arcf != NULL ) delete [] mn1arcf;
00735   if( mnarcf  != NULL ) delete [] mnarcf;
00736   if( mnarcf1 != NULL ) delete [] mnarcf1;
00737   if( mnarcf2 != NULL ) delete [] mnarcf2;
00738   if( mnarcf3 != NULL ) delete [] mnarcf3;
00739   return;
00740 
00741  ERREUR:
00742   if( ierr == 51 || ierr == 52 )
00743   {
00744     //saturation des sommets => redepart avec 2 fois plus de sommets
00745     mxsomm = 2 * mxsomm;
00746     ierr   = 0;
00747     goto NEWDEPART;
00748   }
00749   else
00750   {
00751     MESSAGE( "APTRTE: Triangulation NON REALISEE  avec erreur=" << ierr );
00752     if( ierr == 0 ) ierr=1;
00753     goto NETTOYAGE;
00754   }
00755 }
00756 void
00757 #ifdef WIN32
00758 #ifdef F2C_BUILD
00759 #else
00760               __stdcall
00761 #endif
00762 #endif
00763  qualitetrte( R3 *mnpxyd,
00764                    Z & mosoar, Z & mxsoar, Z *mnsoar,
00765                    Z & moartr, Z & mxartr, Z *mnartr,
00766                    Z & nbtria, R & quamoy, R & quamin )
00767 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00768 // but :    calculer la qualite moyenne et minimale de la triangulation
00769 // -----    actuelle definie par les tableaux mnsoar et mnartr
00770 // entrees:
00771 // --------
00772 // mnpxyd : tableau des coordonnees 2d des points
00773 //          par point : x  y  distance_souhaitee
00774 // mosoar : nombre maximal d'entiers par arete et
00775 //          indice dans mnsoar de l'arete suivante dans le hachage
00776 // mxsoar : nombre maximal d'aretes stockables dans le tableau mnsoar
00777 //          attention: mxsoar>3*mxsomm obligatoire!
00778 // mnsoar : numero des 2 sommets , no ligne, 2 triangles de l'arete,
00779 //          chainage des aretes frontalieres, chainage du hachage des aretes
00780 //          hachage des aretes = mnsoar(1)+mnsoar(2)*2
00781 //          avec mxsoar>=3*mxsomm
00782 //          une arete i de mnsoar est vide <=> mnsoar(1,i)=0 et
00783 //          mnsoar(2,arete vide)=l'arete vide qui precede
00784 //          mnsoar(3,arete vide)=l'arete vide qui suit
00785 // moartr : nombre maximal d'entiers par arete du tableau mnartr
00786 // mxartr : nombre maximal de triangles declarables
00787 // mnartr : les 3 aretes des triangles +-arete1, +-arete2, +-arete3
00788 //          arete1 = 0 si triangle vide => arete2 = triangle vide suivant
00789 // sorties:
00790 // --------
00791 // nbtria : nombre de triangles internes au domaine
00792 // quamoy : qualite moyenne  des triangles actuels
00793 // quamin : qualite minimale des triangles actuels
00794 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
00795 {
00796   R  d, aire, qualite;
00797   Z  nosotr[3], mn, nbtrianeg, nt, ntqmin;
00798 
00799   aire   = 0;
00800   quamoy = 0;
00801   quamin = 2.0;
00802   nbtria = 0;
00803   nbtrianeg = 0;
00804   ntqmin = 0;
00805 
00806   mn = -moartr;
00807   for ( nt=1; nt<=mxartr; nt++ )
00808   {
00809     mn += moartr;
00810     if( mnartr[mn]!=0 )
00811     {
00812       //un triangle occupe de plus
00813       nbtria++;
00814 
00815       //le numero des 3 sommets du triangle nt
00816       nusotr( nt, mosoar, mnsoar, moartr, mnartr,  nosotr );
00817 
00818       //la qualite du triangle ns1 ns2 ns3
00819       qutr2d( mnpxyd[nosotr[0]-1], mnpxyd[nosotr[1]-1], mnpxyd[nosotr[2]-1],
00820                qualite );
00821 
00822       //la qualite moyenne
00823       quamoy += qualite;
00824 
00825       //la qualite minimale
00826       if( qualite < quamin )
00827       {
00828          quamin = qualite;
00829          ntqmin = nt;
00830       }
00831 
00832       //aire signee du triangle nt
00833       d = surtd2( mnpxyd[nosotr[0]-1], mnpxyd[nosotr[1]-1], mnpxyd[nosotr[2]-1] );
00834       if( d<0 )
00835       {
00836         //un triangle d'aire negative de plus
00837         nbtrianeg++;
00838         MESSAGE("ATTENTION: le triangle " << nt << " de sommets:"
00839              << nosotr[0] << " " << nosotr[1] << " " << nosotr[2]
00840              << " a une aire " << d <<"<=0");
00841       }
00842 
00843       //aire des triangles actuels
00844       aire += Abs(d);
00845     }
00846   }
00847 
00848   //les affichages
00849   quamoy /= nbtria;
00850   MESSAGE("Qualite moyenne=" << quamoy
00851        << "  Qualite minimale=" << quamin
00852        << " des " << nbtria << " triangles de surface plane totale="
00853        << aire);
00854 
00855   if( quamin<0.3 )
00856   {
00857     //le numero des 3 sommets du triangle ntqmin de qualite minimale
00858     nusotr(ntqmin, mosoar, mnsoar, moartr, mnartr,  nosotr );
00859     MESSAGE("Triangle de qualite minimale "<<quamin<<" de sommets:"
00860             <<nosotr[0]<<" "<<nosotr[1]<<" "<<nosotr[2]<<" ");
00861     for (int i=0;i<3;i++)
00862       MESSAGE("Sommet "<<nosotr[i]<<": x="<< mnpxyd[nosotr[i]-1].x
00863               <<" y="<< mnpxyd[nosotr[i]-1].y);
00864   }
00865 
00866   if( nbtrianeg>0 )
00867     MESSAGE( "ATTENTION: "<< nbtrianeg << " TRIANGLES d'AIRE NEGATIVE" );
00868 
00869   MESSAGE(" ");
00870   return;
00871 }