Salome HOME
23080: [CEA 1497] Do not merge a middle node in quadratic with the extreme nodes...
[modules/smesh.git] / src / MEFISTO2 / aptrte.cxx
index 46429ead71ba0df4cf1c95b986a9ecd3c0b625bc..ca2474262911681b847433498632c4a1f640baa1 100755 (executable)
@@ -1,27 +1,28 @@
 //  MEFISTO2: a library to compute 2D triangulation from segmented boundaries
 //
-//  Copyright (C) 2006  Laboratoire J.-L. Lions UPMC Paris
 //
-//  This library is free software; you can redistribute it and/or
-//  modify it under the terms of the GNU Lesser General Public
-//  License as published by the Free Software Foundation; either
-//  version 2.1 of the License.
+// Copyright (C) 2006-2015  CEA/DEN, EDF R&D, OPEN CASCADE
 //
-//  This library is distributed in the hope that it will be useful,
-//  but WITHOUT ANY WARRANTY; without even the implied warranty of
-//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
-//  Lesser General Public License for more details.
+// This library is free software; you can redistribute it and/or
+// modify it under the terms of the GNU Lesser General Public
+// License as published by the Free Software Foundation; either
+// version 2.1 of the License, or (at your option) any later version.
 //
-//  You should have received a copy of the GNU Lesser General Public
-//  License along with this library; if not, write to the Free Software
-//  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+// This library is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+// Lesser General Public License for more details.
 //
-//  See http://www.ann.jussieu.fr/~perronne or email Perronnet@ann.jussieu.fr
+// You should have received a copy of the GNU Lesser General Public
+// License along with this library; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
+//
+// See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
 //
 //  File   : aptrte.cxx   le C++ de l'appel du trianguleur plan
 //  Module : SMESH
 //  Author : Alain PERRONNET
-//  Date   : 16 mars 2006
+//  Date   : 13 novembre 2006
 
 #include "Rn.h"
 #include "aptrte.h"
@@ -32,7 +33,15 @@ using namespace std;
 extern "C"
 {
   R aretemaxface_;
-  R areteideale_( R3 xyz, R3 direction )
+  MEFISTO2D_EXPORT   
+    R
+  #ifdef WIN32
+  #ifdef F2C_BUILD
+  #else
+      __stdcall
+  #endif
+  #endif
+      areteideale()//( R3 xyz, R3 direction )
   {
     return aretemaxface_;
   }
@@ -44,7 +53,14 @@ extern "C"
 
 static double cpunew, cpuold=0;
 
-void tempscpu_( double & tempsec )
+void
+#ifdef WIN32
+#ifdef F2C_BUILD
+#else
+              __stdcall
+#endif
+#endif
+tempscpu_( double & tempsec )
 //Retourne le temps CPU utilise en secondes
 {  
   tempsec = ( (double) clock() ) / CLOCKS_PER_SEC;
@@ -52,7 +68,14 @@ void tempscpu_( double & tempsec )
 }
 
 
-void deltacpu_( R & dtcpu )
+void
+#ifdef WIN32
+#ifdef F2C_BUILD
+#else
+              __stdcall
+#endif
+#endif
+deltacpu_( R & dtcpu )
 //Retourne le temps CPU utilise en secondes depuis le precedent appel
 {
   tempscpu_( cpunew );
@@ -64,11 +87,11 @@ void deltacpu_( R & dtcpu )
 
 
 void  aptrte( Z   nutysu, R      aretmx,
-             Z   nblf,   Z  *   nudslf,  R2 * uvslf,
-             Z   nbpti,  R2 *   uvpti,
-             Z & nbst,   R2 * & uvst,
-             Z & nbt,    Z  * & nust,
-             Z & ierr )
+              Z   nblf,   Z  *   nudslf,  R2 * uvslf,
+              Z   nbpti,  R2 *   uvpti,
+              Z & nbst,   R2 * & uvst,
+              Z & nbt,    Z  * & nust,
+              Z & ierr )
 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 // but : appel de la triangulation par un arbre-4 recouvrant
 // ----- de triangles equilateraux
@@ -118,7 +141,7 @@ void  aptrte( Z   nutysu, R      aretmx,
                  //no st1, st2, st3, 0 (non quadrangle)
 
   R  d, tcpu=0;
-  R3 direction=R3(0,0,0);  //direction pour areteideale() inactive ici!
+//  R3 direction=R3(0,0,0);  //direction pour areteideale() inactive ici!
   Z  nbarfr=nudslf[nblf];  //nombre total d'aretes des lignes fermees
   Z  mxtrou = Max( 1024, nblf );  //nombre maximal de trous dans la surface
 
@@ -139,15 +162,14 @@ void  aptrte( Z   nutysu, R      aretmx,
 
   R3 comxmi[2];            //coordonnees UV Min et Maximales
   R  aremin, aremax;       //longueur minimale et maximale des aretes
+  R  airemx;               //aire maximale souhaitee d'un triangle
   R  quamoy, quamin;
 
   Z  noar0, noar, na;
   Z  i, l, n, ns, ns0, ns1, ns2, nosotr[3], nt;
   Z  mxsomm, nbsomm, nbarpi, nbarli, ndtri0, mn;
   Z  moins1=-1;
-  R  dist;
-
-  aretemaxface_ = aretmx;
+  Z  nuds = 0;
 
   // initialisation du temps cpu
   deltacpu_( d );
@@ -161,7 +183,8 @@ void  aptrte( Z   nutysu, R      aretmx,
   i = 4*nbarfr/10;
   mxsomm = Max( 20000, 64*nbpti+i*i );
   MESSAGE( "APTRTE: Debut de la triangulation plane avec " );
-  MESSAGE( "nutysu=" << nutysu << "  aretmx=" << aretmx << "  mxsomm=" << mxsomm );
+  MESSAGE( "nutysu=" << nutysu << "  aretmx=" << aretmx
+           << "  mxsomm=" << mxsomm );
   MESSAGE( nbarfr << " sommets sur la frontiere et " << nbpti << " points internes");
 
  NEWDEPART:
@@ -187,14 +210,14 @@ void  aptrte( Z   nutysu, R      aretmx,
   mnsoar = new Z[mosoar*mxsoar];
   if( mnsoar==NULL ) goto ERREUR;
   //initialiser le tableau mnsoar pour le hachage des aretes
-  insoar_( mxsomm, mosoar, mxsoar, n1soar, mnsoar );
+  insoar( mxsomm, mosoar, mxsoar, n1soar, mnsoar );
 
   // mnarst( mxsomm ) numero mnsoar d'une arete pour chacun des sommets
   if( mnarst!=NULL ) delete [] mnarst;
   mnarst = new Z[1+mxsomm];
   if( mnarst==NULL ) goto ERREUR;
   n = 1+mxsomm;
-  azeroi_( n, mnarst );
+  azeroi( n, mnarst );
 
   // mnslig( mxsomm ) no de sommet dans sa ligne pour chaque sommet frontalier
   //               ou no du point si interne forc'e par l'utilisateur
@@ -202,7 +225,7 @@ void  aptrte( Z   nutysu, R      aretmx,
   if( mnslig!=NULL ) delete [] mnslig;
   mnslig = new Z[mxsomm];
   if( mnslig==NULL ) goto ERREUR;
-  azeroi_( mxsomm, mnslig );
+  azeroi( mxsomm, mnslig );
 
   // initialisation des aretes frontalieres de la triangulation future
   // renumerotation des sommets des aretes des lignes pour la triangulation
@@ -222,9 +245,9 @@ void  aptrte( Z   nutysu, R      aretmx,
     ns0 = nudslf[n-1];
     mnpxyd[ns0].x = uvslf[ns0].x;
     mnpxyd[ns0].y = uvslf[ns0].y;
-    mnpxyd[ns0].z = areteideale_( mnpxyd[ns0], direction );
+    mnpxyd[ns0].z = areteideale();//( mnpxyd[ns0], direction );
 //     MESSAGE("Sommet " << ns0 << ": " << mnpxyd[ns0].x
-//      << " " << mnpxyd[ns0].y << " longueur arete=" << mnpxyd[ns0].z);
+//       << " " << mnpxyd[ns0].y << " longueur arete=" << mnpxyd[ns0].z);
 
     //carre de la longueur de l'arete 1 de la ligne fermee n
     d = pow( uvslf[ns0+1].x - uvslf[ns0].x, 2 ) 
@@ -244,9 +267,9 @@ void  aptrte( Z   nutysu, R      aretmx,
 
      //le numero n de la ligne du sommet et son numero ns1 dans la ligne
     mnslig[ns0-1] = 1000000 * n + ns1-nudslf[n-1];
-    fasoar_( ns1, ns2, moins1, moins1, n,
-            mosoar, mxsoar, n1soar, mnsoar, mnarst,
-            noar0,  ierr );
+    fasoar( ns1, ns2, moins1, moins1, n,
+             mosoar, mxsoar, n1soar, mnsoar, mnarst,
+             noar0,  ierr );
     //pas de test sur ierr car pas de saturation possible a ce niveau
 
     //le pointeur dans le hachage sur la premiere arete de la ligne fermee n
@@ -263,20 +286,23 @@ void  aptrte( Z   nutysu, R      aretmx,
     {
       ns1 = ns2; //le numero de l'arete et le numero du premier sommet de l'arete
       if( i < nbarli )
-       //nbs+1 est le 2-eme sommet de l'arete i de la ligne fermee n
-       ns2 = ns1+1;
+        //nbs+1 est le 2-eme sommet de l'arete i de la ligne fermee n
+        ns2 = ns1+1;
       else
-       //le 2-eme sommet de la derniere arete est le premier sommet de la ligne
-       ns2 = ns0;
+        //le 2-eme sommet de la derniere arete est le premier sommet de la ligne
+        ns2 = ns0;
 
       //l'arete precedente est dotee de sa suivante:celle cree ensuite
       //les 2 coordonnees du sommet ns2 de la ligne
       ns = ns1 - 1;
+//debut ajout  5/10/2006  ................................................
+      nuds = Max( nuds, ns );   //le numero du dernier sommet traite
+//fin   ajout  5/10/2006  ................................................
       mnpxyd[ns].x = uvslf[ns].x;
       mnpxyd[ns].y = uvslf[ns].y;
-      mnpxyd[ns].z = areteideale_( mnpxyd[ns], direction );
+      mnpxyd[ns].z = areteideale();//( mnpxyd[ns], direction );
 //       MESSAGE("Sommet " << ns << ": " << mnpxyd[ns].x
-//        << " " << mnpxyd[ns].y << " longueur arete=" << mnpxyd[ns].z);
+//         << " " << mnpxyd[ns].y << " longueur arete=" << mnpxyd[ns].z);
 
       //carre de la longueur de l'arete
       d = pow( uvslf[ns2-1].x - uvslf[ns1-1].x, 2) 
@@ -284,13 +310,21 @@ void  aptrte( Z   nutysu, R      aretmx,
       aremin = Min( aremin, d );
       aremax = Max( aremax, d );
 
+//debut ajout du 5/10/2006  .............................................
+      //la longueur de l'arete ns1-ns2
+      d = sqrt( d );
+      //longueur arete = Min ( aretmx, aretes incidentes )
+      mnpxyd[ns   ].z = Min( mnpxyd[ns   ].z, d );
+      mnpxyd[ns2-1].z = Min( mnpxyd[ns2-1].z, d );
+//fin ajout du 5/10/2006  ...............................................
+
       //le numero n de la ligne du sommet et son numero ns1 dans la ligne
       mnslig[ns] = 1000000 * n + ns1-nudslf[n-1];
 
       //ajout de l'arete dans la liste
-      fasoar_( ns1, ns2, moins1, moins1, n,
-              mosoar, mxsoar, n1soar, mnsoar,
-              mnarst, noar, ierr );
+      fasoar( ns1, ns2, moins1, moins1, n,
+               mosoar, mxsoar, n1soar, mnsoar,
+               mnarst, noar, ierr );
       //pas de test sur ierr car pas de saturation possible a ce niveau
 
       //chainage des aretes frontalieres en position 6 du tableau mnsoar
@@ -307,9 +341,34 @@ void  aptrte( Z   nutysu, R      aretmx,
   aremin = sqrt( aremin );  //longueur minimale d'une arete des lignes fermees
   aremax = sqrt( aremax );  //longueur maximale d'une arete
 
-  aretmx = Min( aretmx, aremax );  //pour homogeneiser
-  MESSAGE("nutysu=" << nutysu << "  aretmx=" << aretmx 
-       << "  arete min=" << aremin << "  arete max=" << aremax);
+//debut ajout  9/11/2006  ................................................
+  // devenu un commentaire aretmx = Min( aretmx, aremax ); //pour homogeneiser
+
+  // protection contre une arete max desiree trop grande ou trop petite
+  if( aretmx > aremax*2.05 ) aretmx = aremax;
+
+  // protection contre une arete max desiree trop petite
+  if( (aremax-aremin) > (aremin+aremax)*0.05 && aretmx < aremin*0.5 )
+    aretmx =(aremin+aremax*2)/3.0;
+
+  if( aretmx < aremin  && aremin > 0 )
+    aretmx = aremin;
+
+  //sauvegarde pour la fonction areteideale_
+  aretemaxface_ = aretmx;
+
+  //aire maximale souhaitee des triangles
+  airemx = aretmx * aretmx * sqrt(3.0) / 2.0;  //Aire triangle equilateral
+
+  for(i=0; i<=nuds; i++ )
+    mnpxyd[i].z = Min( mnpxyd[i].z, aretmx );
+  //MESSAGE("Numero du dernier sommet frontalier=" << nuds+1);
+//fin  ajout 9/11/2006  .................................................
+
+
+  MESSAGE("Sur  le  bord: arete min=" << aremin << " arete max=" << aremax );
+  MESSAGE("Triangulation: arete mx=" << aretmx
+          << " triangle aire mx=" << airemx );
 
   //chainage des aretes frontalieres : la derniere arete frontaliere
   mnsoar[ mosoar * noar - mosoar + 5 ] = 0;
@@ -333,7 +392,7 @@ void  aptrte( Z   nutysu, R      aretmx,
     //les 2 coordonnees du point i de sommet nbs
     mnpxyd[ns1].x = uvpti[i].x;
     mnpxyd[ns1].y = uvpti[i].y;
-    mnpxyd[ns1].z = areteideale_( mnpxyd[ns1], direction );
+    mnpxyd[ns1].z = areteideale();//( mnpxyd[ns1], direction );
     //le numero i du point interne
     mnslig[ns1] = i+1;
     ns1++;
@@ -356,7 +415,9 @@ void  aptrte( Z   nutysu, R      aretmx,
   if( mntree==NULL ) goto ERREUR;
 
   //initialisation du tableau letree et ajout dans letree des sommets 1 a nbsomm
-  teajte_( mxsomm, nbsomm, mnpxyd, comxmi, aretmx, mxtree, mntree, ierr );
+  comxmi[0].x = comxmi[1].x = uvslf[0].x;
+  comxmi[0].y = comxmi[1].y = uvslf[0].y;
+  teajte( mxsomm, nbsomm, mnpxyd, comxmi, aretmx, mxtree, mntree, ierr );
   comxmi[0].z=0;
   comxmi[1].z=0;
 
@@ -384,10 +445,10 @@ void  aptrte( Z   nutysu, R      aretmx,
   mnqueu = new Z[mxqueu];
   if( mnqueu==NULL) goto ERREUR;
 
-  tehote_( nutysu, nbarpi, mxsomm, nbsomm, mnpxyd,
-          comxmi, aretmx,
-          mntree, mxqueu, mnqueu,
-          ierr );
+  tehote( nutysu, nbarpi, mxsomm, nbsomm, mnpxyd,
+           comxmi, aretmx,
+           mntree, mxqueu, mnqueu,
+           ierr );
 
   deltacpu_( d );
   tcpu += d;
@@ -411,10 +472,10 @@ void  aptrte( Z   nutysu, R      aretmx,
   // trianguler les triangles equilateraux feuilles a partir de leurs 3 sommets
   // et des points de la frontiere, des points internes imposes interieurs
   // ==========================================================================
-  tetrte_( comxmi, aretmx, nbarpi, mxsomm, mnpxyd,
-          mxqueu, mnqueu, mntree, mosoar, mxsoar, n1soar, mnsoar,
-          moartr, mxartr, n1artr, mnartr, mnarst,
-          ierr );
+  tetrte( comxmi, aretmx, nbarpi, mxsomm, mnpxyd,
+           mxqueu, mnqueu, mntree, mosoar, mxsoar, n1soar, mnsoar,
+           moartr, mxartr, n1artr, mnartr, mnarst,
+           ierr );
 
   // destruction de la queue et de l'arbre devenus inutiles
   delete [] mnqueu;  mnqueu=NULL;
@@ -434,16 +495,16 @@ void  aptrte( Z   nutysu, R      aretmx,
 
   //qualites de la triangulation actuelle
   qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
-              nbt, quamoy, quamin );
+                nbt, quamoy, quamin );
 
   // boucle sur les aretes internes (non sur une ligne de la frontiere)
   // avec echange des 2 diagonales afin de rendre la triangulation delaunay
   // ======================================================================
   // formation du chainage 6 des aretes internes a echanger eventuellement
-  aisoar_( mosoar, mxsoar, mnsoar, na );
-  tedela_( mnpxyd, mnarst,
-          mosoar, mxsoar, n1soar, mnsoar, na,
-          moartr, mxartr, n1artr, mnartr, n );
+  aisoar( mosoar, mxsoar, mnsoar, na );
+  tedela( mnpxyd, mnarst,
+           mosoar, mxsoar, n1soar, mnsoar, na,
+           moartr, mxartr, n1artr, mnartr, n );
 
   MESSAGE( "Nombre d'echanges des diagonales de 2 triangles=" << n );
   deltacpu_( d );
@@ -453,7 +514,7 @@ void  aptrte( Z   nutysu, R      aretmx,
 
   //qualites de la triangulation actuelle
   qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
-              nbt, quamoy, quamin );
+                nbt, quamoy, quamin );
 
   // detection des aretes frontalieres initiales perdues
   // triangulation frontale pour les restaurer
@@ -472,13 +533,13 @@ void  aptrte( Z   nutysu, R      aretmx,
   mnarcf2 = new Z[mxarcf];
   if( mnarcf2 == NULL ) goto ERREUR;
 
-  terefr_( nbarpi, mnpxyd,
-          mosoar, mxsoar, n1soar, mnsoar,
-          moartr, n1artr, mnartr, mnarst,
-          mxarcf, mn1arcf, mnarcf, mnarcf1, mnarcf2,
-          n, ierr );
+  terefr( nbarpi, mnpxyd,
+           mosoar, mxsoar, n1soar, mnsoar,
+           moartr, mxartr, n1artr, mnartr, mnarst,
+           mxarcf, mn1arcf, mnarcf, mnarcf1, mnarcf2,
+           n, ierr );
 
-  MESSAGE( "Restauration de " << n << " aretes perdues de la frontiere" );
+  MESSAGE( "Restauration de " << n << " aretes perdues de la frontiere  ierr=" << ierr );
   deltacpu_( d );
   tcpu += d;
   MESSAGE("Temps de la recuperation des aretes perdues de la frontiere="
@@ -488,7 +549,7 @@ void  aptrte( Z   nutysu, R      aretmx,
 
   //qualites de la triangulation actuelle
   qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
-              nbt, quamoy, quamin );
+                nbt, quamoy, quamin );
 
   // fin de la triangulation avec respect des aretes initiales frontalieres
 
@@ -513,23 +574,23 @@ void  aptrte( Z   nutysu, R      aretmx,
   for (n=0; n<nblf; n++)  //numero de la ligne fermee de 1 a nblf
     mnlftr[n] = n+1;
 
-  tesuex_( nblf,   mnlftr,
-          ndtri0, nbsomm, mnpxyd, mnslig,
-          mosoar, mxsoar, mnsoar,
-          moartr, mxartr, n1artr, mnartr, mnarst,
-          nbt, mntrsu, ierr );
+  tesuex( nblf,   mnlftr,
+           ndtri0, nbsomm, mnpxyd, mnslig,
+           mosoar, mxsoar, mnsoar,
+           moartr, mxartr, n1artr, mnartr, mnarst,
+           nbt, mntrsu, ierr );
 
   delete [] mnlftr; mnlftr=NULL;
   delete [] mntrsu; mntrsu=NULL;
 
   deltacpu_( d );
   tcpu += d;
-  MESSAGE( "Temps de la suppression des triangles externes=" << d );
+  MESSAGE( "Temps de la suppression des triangles externes=" << d << "ierr=" << ierr );
   if( ierr != 0 ) goto ERREUR;
 
   //qualites de la triangulation actuelle
   qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
-              nbt, quamoy, quamin );
+                nbt, quamoy, quamin );
 
   // amelioration de la qualite de la triangulation par
   // barycentrage des sommets internes a la triangulation
@@ -540,16 +601,16 @@ void  aptrte( Z   nutysu, R      aretmx,
   mnarcf3 = new Z[mxarcf];
   if( mnarcf3 == NULL )
   {
-    cout << "aptrte: MC saturee mnarcf3=" << mnarcf3 << endl;
+    MESSAGE ( "aptrte: MC saturee mnarcf3=" << mnarcf3 );
     goto ERREUR;
   }
-  teamqt_( nutysu,
-          mnarst, mosoar, mxsoar, n1soar, mnsoar,
-          moartr, mxartr, n1artr, mnartr,
-          mxarcf, mnarcf2, mnarcf3,
-          mn1arcf, mnarcf, mnarcf1,
-          comxmi, nbarpi, nbsomm, mxsomm, mnpxyd, mnslig,
-          ierr );
+  teamqt( nutysu,  aretmx,  airemx,
+           mnarst,  mosoar,  mxsoar, n1soar, mnsoar,
+           moartr,  mxartr,  n1artr, mnartr,
+           mxarcf,  mnarcf2, mnarcf3,
+           mn1arcf, mnarcf,  mnarcf1,
+           nbarpi,  nbsomm, mxsomm, mnpxyd, mnslig,
+           ierr );
   if( mnarcf3 != NULL ) {delete [] mnarcf3; mnarcf3=NULL;}
   if( mn1arcf != NULL ) {delete [] mn1arcf; mn1arcf=NULL;}
   if( mnarcf  != NULL ) {delete [] mnarcf;  mnarcf =NULL;}
@@ -559,11 +620,12 @@ void  aptrte( Z   nutysu, R      aretmx,
   deltacpu_( d );
   tcpu += d;
   MESSAGE( "Temps de l'amelioration de la qualite de la triangulation=" << d );
-  if( ierr != 0 ) goto ERREUR;
+  if( ierr == -13 ) ierr=0; //6/10/2006 arret de l'amelioration apres boucle infinie dans caetoi
+  if( ierr !=   0 ) goto ERREUR;
 
   //qualites de la triangulation finale
   qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
-              nbt, quamoy, quamin );
+                nbt, quamoy, quamin );
 
   // renumerotation des sommets internes: mnarst(i)=numero final du sommet
   // ===================================
@@ -575,7 +637,7 @@ void  aptrte( Z   nutysu, R      aretmx,
     if( mnartr[nt*moartr-moartr] != 0 )
     {
       //le numero des 3 sommets du triangle nt
-      nusotr_( nt, mosoar, mnsoar, moartr, mnartr, nosotr );
+      nusotr( nt, mosoar, mnsoar, moartr, mnartr, nosotr );
       //les 3 sommets du triangle sont actifs
       mnarst[ nosotr[0] ] = 1;
       mnarst[ nosotr[1] ] = 1;
@@ -609,22 +671,22 @@ void  aptrte( Z   nutysu, R      aretmx,
       n = mnslig[i];
       if( n > 0 )
       {
-       if( n >= 1000000 )
-       {
-         //sommet d'une ligne
-         //retour aux coordonnees initiales dans uvslf
-         l = n / 1000000;
-         n = n - 1000000 * l + nudslf[l-1] - 1;
-         uvst[nbst].x = uvslf[n].x;
-         uvst[nbst].y = uvslf[n].y;
-       }
-       else
-       {
-         //point utilisateur n interne impose
-         //retour aux coordonnees initiales dans uvpti
-         uvst[nbst].x = uvpti[n-1].x;
-         uvst[nbst].y = uvpti[n-1].y;
-       }
+        if( n >= 1000000 )
+        {
+          //sommet d'une ligne
+          //retour aux coordonnees initiales dans uvslf
+          l = n / 1000000;
+          n = n - 1000000 * l + nudslf[l-1] - 1;
+          uvst[nbst].x = uvslf[n].x;
+          uvst[nbst].y = uvslf[n].y;
+        }
+        else
+        {
+          //point utilisateur n interne impose
+          //retour aux coordonnees initiales dans uvpti
+          uvst[nbst].x = uvpti[n-1].x;
+          uvst[nbst].y = uvpti[n-1].y;
+        }
       }
     }
   }
@@ -643,7 +705,7 @@ void  aptrte( Z   nutysu, R      aretmx,
     if( mnartr[i*moartr-moartr] != 0 )
     {
       //le triangle i est interne => nosotr numero de ses 3 sommets
-      nusotr_( i, mosoar, mnsoar, moartr, mnartr,  nosotr );
+      nusotr( i, mosoar, mnsoar, moartr, mnartr,  nosotr );
       nust[nbt++] = mnarst[ nosotr[0] ];
       nust[nbt++] = mnarst[ nosotr[1] ];
       nust[nbt++] = mnarst[ nosotr[2] ];
@@ -652,7 +714,7 @@ void  aptrte( Z   nutysu, R      aretmx,
   }
   nbt /= nbsttria;  //le nombre final de triangles de la surface
   MESSAGE( "APTRTE: Fin de la triangulation plane avec "<<nbst<<" sommets et "
-          << nbt << " triangles=" << nbt);
+           << nbt << " triangles" );
   deltacpu_( d );
   tcpu += d;
   MESSAGE( "APTRTE: Temps total de la triangulation plane=" << tcpu << " secondes" );
@@ -691,12 +753,17 @@ void  aptrte( Z   nutysu, R      aretmx,
     goto NETTOYAGE;
   }
 }
-
-
-void qualitetrte( R3 *mnpxyd,
-                 Z & mosoar, Z & mxsoar, Z *mnsoar,
-                 Z & moartr, Z & mxartr, Z *mnartr,
-                 Z & nbtria, R & quamoy, R & quamin )
+void
+#ifdef WIN32
+#ifdef F2C_BUILD
+#else
+              __stdcall
+#endif
+#endif
+ qualitetrte( R3 *mnpxyd,
+                   Z & mosoar, Z & mxsoar, Z *mnsoar,
+                   Z & moartr, Z & mxartr, Z *mnartr,
+                   Z & nbtria, R & quamoy, R & quamin )
 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 // but :    calculer la qualite moyenne et minimale de la triangulation
 // -----    actuelle definie par les tableaux mnsoar et mnartr
@@ -727,13 +794,14 @@ void qualitetrte( R3 *mnpxyd,
 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 {
   R  d, aire, qualite;
-  Z  nosotr[3], mn, nbtrianeg, nt;
+  Z  nosotr[3], mn, nbtrianeg, nt, ntqmin;
 
   aire   = 0;
   quamoy = 0;
   quamin = 2.0;
   nbtria = 0;
   nbtrianeg = 0;
+  ntqmin = 0;
 
   mn = -moartr;
   for ( nt=1; nt<=mxartr; nt++ )
@@ -745,27 +813,31 @@ void qualitetrte( R3 *mnpxyd,
       nbtria++;
 
       //le numero des 3 sommets du triangle nt
-      nusotr_( nt, mosoar, mnsoar, moartr, mnartr,  nosotr );
+      nusotr( nt, mosoar, mnsoar, moartr, mnartr,  nosotr );
 
       //la qualite du triangle ns1 ns2 ns3
-      qutr2d_( mnpxyd[nosotr[0]-1], mnpxyd[nosotr[1]-1], mnpxyd[nosotr[2]-1],
-              qualite );
+      qutr2d( mnpxyd[nosotr[0]-1], mnpxyd[nosotr[1]-1], mnpxyd[nosotr[2]-1],
+               qualite );
 
       //la qualite moyenne
       quamoy += qualite;
 
       //la qualite minimale
-      quamin = Min( quamin, qualite );
+      if( qualite < quamin )
+      {
+         quamin = qualite;
+         ntqmin = nt;
+      }
 
       //aire signee du triangle nt
-      d = surtd2_( mnpxyd[nosotr[0]-1], mnpxyd[nosotr[1]-1], mnpxyd[nosotr[2]-1] );
+      d = surtd2( mnpxyd[nosotr[0]-1], mnpxyd[nosotr[1]-1], mnpxyd[nosotr[2]-1] );
       if( d<0 )
       {
-       //un triangle d'aire negative de plus
-       nbtrianeg++;
-       MESSAGE("ATTENTION: le triangle " << nt << " de sommets:"
-            << nosotr[0] << " " << nosotr[1] << " " << nosotr[2]
-            << " a une aire " << d <<"<=0");
+        //un triangle d'aire negative de plus
+        nbtrianeg++;
+        MESSAGE("ATTENTION: le triangle " << nt << " de sommets:"
+             << nosotr[0] << " " << nosotr[1] << " " << nosotr[2]
+             << " a une aire " << d <<"<=0");
       }
 
       //aire des triangles actuels
@@ -780,7 +852,20 @@ void qualitetrte( R3 *mnpxyd,
        << " des " << nbtria << " triangles de surface plane totale="
        << aire);
 
+  if( quamin<0.3 )
+  {
+    //le numero des 3 sommets du triangle ntqmin de qualite minimale
+    nusotr(ntqmin, mosoar, mnsoar, moartr, mnartr,  nosotr );
+    MESSAGE("Triangle de qualite minimale "<<quamin<<" de sommets:"
+            <<nosotr[0]<<" "<<nosotr[1]<<" "<<nosotr[2]<<" ");
+    for (int i=0;i<3;i++)
+      MESSAGE("Sommet "<<nosotr[i]<<": x="<< mnpxyd[nosotr[i]-1].x
+              <<" y="<< mnpxyd[nosotr[i]-1].y);
+  }
+
   if( nbtrianeg>0 )
-    MESSAGE( "ATTENTION: nombre de triangles d'aire negative=" << nbtrianeg );
+    MESSAGE( "ATTENTION: "<< nbtrianeg << " TRIANGLES d'AIRE NEGATIVE" );
+
+  MESSAGE(" ");
   return;
 }