Salome HOME
PAL14419 (IMP: a filter predicate to find nodes/elements lying on any
[modules/smesh.git] / src / MEFISTO2 / aptrte.cxx
index 58da98168175e2835bd00c924cc51a64c728ab1f..46429ead71ba0df4cf1c95b986a9ecd3c0b625bc 100755 (executable)
@@ -1,8 +1,34 @@
-using namespace std;
+//  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.
+//
+//  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.
+//
+//  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.ann.jussieu.fr/~perronne or email Perronnet@ann.jussieu.fr
+//
+//  File   : aptrte.cxx   le C++ de l'appel du trianguleur plan
+//  Module : SMESH
+//  Author : Alain PERRONNET
+//  Date   : 16 mars 2006
+
 #include "Rn.h"
 #include "aptrte.h"
 #include "utilities.h"
 
+using namespace std;
+
 extern "C"
 {
   R aretemaxface_;
@@ -37,10 +63,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,
+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 )
 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 // but : appel de la triangulation par un arbre-4 recouvrant
@@ -84,9 +111,12 @@ void  aptrte( Z nutysu, R aretmx,
 // ierr   : 0 si pas d'erreur
 //        > 0 sinon
 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-// auteur : Alain Perronnet  Analyse Numerique Paris UPMC   decembre 2001
+// auteur : Alain Perronnet  Laboratoire J.-L. LIONS Paris UPMC mars 2006
 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 {
+  Z  nbsttria=4; //Attention: 4 sommets stockes par triangle
+                 //no st1, st2, st3, 0 (non quadrangle)
+
   R  d, tcpu=0;
   R3 direction=R3(0,0,0);  //direction pour areteideale() inactive ici!
   Z  nbarfr=nudslf[nblf];  //nombre total d'aretes des lignes fermees
@@ -95,7 +125,7 @@ void  aptrte( Z nutysu, R aretmx,
   R3 *mnpxyd=NULL;
   Z  *mnsoar=NULL, mosoar=7, mxsoar, n1soar; //le hachage des aretes
   Z  *mnartr=NULL, moartr=3, mxartr, n1artr; //le no des 3 aretes des triangles
-  Z  *mntree=NULL, motree=9, mxtree;   //L'arbre 4 de TE et nombre d'entiers par TE
+  Z  *mntree=NULL, motree=9, mxtree; //L'arbre 4 de TE et nombre d'entiers par TE
   Z  *mnqueu=NULL, mxqueu;
   Z  *mn1arcf=NULL;
   Z  *mnarcf=NULL, mxarcf;
@@ -103,7 +133,6 @@ void  aptrte( Z nutysu, R aretmx,
   Z  *mnarcf2=NULL;
   Z  *mnarcf3=NULL;
   Z  *mntrsu=NULL;
-  Z  *mndalf=NULL;
   Z  *mnslig=NULL;
   Z  *mnarst=NULL;
   Z  *mnlftr=NULL;
@@ -116,6 +145,7 @@ void  aptrte( Z nutysu, R aretmx,
   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;
 
@@ -125,19 +155,14 @@ void  aptrte( Z nutysu, R aretmx,
 
   // quelques reservations de tableaux pour faire les calculs
   // ========================================================
-  // le tableau pointeur sur la premiere arete de chaque ligne fermee
-  if( mndalf!=NULL ) delete [] mndalf;
-  mndalf = new Z[1+nblf];
-  if( mndalf==NULL ) goto ERREUR;
-  mndalf[0]=0;
-
   // declaration du tableau des coordonnees des sommets de la frontiere
   // puis des sommets internes ajoutes
   // majoration empirique du nombre de sommets de la triangulation
-  i =  4*nbarfr/10;
+  i = 4*nbarfr/10;
   mxsomm = Max( 20000, 64*nbpti+i*i );
-  MESSAGE( "APTRTE: Depart de la triangulation avec " );
+  MESSAGE( "APTRTE: Debut de la triangulation plane avec " );
   MESSAGE( "nutysu=" << nutysu << "  aretmx=" << aretmx << "  mxsomm=" << mxsomm );
+  MESSAGE( nbarfr << " sommets sur la frontiere et " << nbpti << " points internes");
 
  NEWDEPART:
   //mnpxyd( 3, mxsomm ) les coordonnees UV des sommets et la taille d'arete aux sommets
@@ -198,8 +223,8 @@ void  aptrte( Z nutysu, R aretmx,
     mnpxyd[ns0].x = uvslf[ns0].x;
     mnpxyd[ns0].y = uvslf[ns0].y;
     mnpxyd[ns0].z = areteideale_( mnpxyd[ns0], direction );
-//     cout << "Sommet " << ns0 << ": " << mnpxyd[ns0].x
-//      << " " << mnpxyd[ns0].y << " longueur arete=" << mnpxyd[ns0].z << endl;
+//     MESSAGE("Sommet " << ns0 << ": " << mnpxyd[ns0].x
+//      << " " << 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 ) 
@@ -225,7 +250,7 @@ void  aptrte( Z nutysu, R aretmx,
     //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
-    mndalf[n] = noar0;
+    //mndalf[n] = noar0;
 
     //la nouvelle arete est la suivante de l'arete definie juste avant
     if( noar > 0 )
@@ -250,8 +275,8 @@ void  aptrte( Z nutysu, R aretmx,
       mnpxyd[ns].x = uvslf[ns].x;
       mnpxyd[ns].y = uvslf[ns].y;
       mnpxyd[ns].z = areteideale_( mnpxyd[ns], direction );
-//       cout << "Sommet " << ns << ": " << mnpxyd[ns].x
-//        << " " << mnpxyd[ns].y << " longueur arete=" << mnpxyd[ns].z << endl;
+//       MESSAGE("Sommet " << ns << ": " << mnpxyd[ns].x
+//        << " " << 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) 
@@ -283,8 +308,8 @@ void  aptrte( Z nutysu, R aretmx,
   aremax = sqrt( aremax );  //longueur maximale d'une arete
 
   aretmx = Min( aretmx, aremax );  //pour homogeneiser
-  cout << "nutysu=" << nutysu << "  aretmx=" << aretmx 
-       << "  arete min=" << aremin << "  arete max=" << aremax << endl;
+  MESSAGE("nutysu=" << nutysu << "  aretmx=" << aretmx 
+       << "  arete min=" << aremin << "  arete max=" << aremax);
 
   //chainage des aretes frontalieres : la derniere arete frontaliere
   mnsoar[ mosoar * noar - mosoar + 5 ] = 0;
@@ -366,8 +391,8 @@ void  aptrte( Z nutysu, R aretmx,
 
   deltacpu_( d );
   tcpu += d;
-  cout << "Temps de l'adaptation et l'homogeneisation de l'arbre-4 des TE="
-       << d << " secondes" << endl;
+  MESSAGE("Temps de l'adaptation et l'homogeneisation de l'arbre-4 des TE="
+       << d << " secondes");
   if( ierr != 0 )
   {
     //destruction du tableau auxiliaire et de l'arbre
@@ -423,8 +448,8 @@ void  aptrte( Z nutysu, R aretmx,
   MESSAGE( "Nombre d'echanges des diagonales de 2 triangles=" << n );
   deltacpu_( d );
   tcpu += d;
-  cout << "Temps de la triangulation Delaunay par echange des diagonales="
-       << d << " secondes" << endl;
+  MESSAGE("Temps de la triangulation Delaunay par echange des diagonales="
+       << d << " secondes");
 
   //qualites de la triangulation actuelle
   qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
@@ -456,8 +481,8 @@ void  aptrte( Z nutysu, R aretmx,
   MESSAGE( "Restauration de " << n << " aretes perdues de la frontiere" );
   deltacpu_( d );
   tcpu += d;
-  cout << "Temps de la recuperation des aretes perdues de la frontiere="
-       << d << " secondes" << endl;
+  MESSAGE("Temps de la recuperation des aretes perdues de la frontiere="
+       << d << " secondes");
 
   if( ierr != 0 ) goto ERREUR;
 
@@ -513,8 +538,11 @@ void  aptrte( Z nutysu, R aretmx,
   // mise en delaunay de la triangulation
   // =====================================================
   mnarcf3 = new Z[mxarcf];
-  if( mnarcf3 == NULL ) goto ERREUR;
-
+  if( mnarcf3 == NULL )
+  {
+    cout << "aptrte: MC saturee mnarcf3=" << mnarcf3 << endl;
+    goto ERREUR;
+  }
   teamqt_( nutysu,
           mnarst, mosoar, mxsoar, n1soar, mnsoar,
           moartr, mxartr, n1artr, mnartr,
@@ -522,11 +550,11 @@ void  aptrte( Z nutysu, R aretmx,
           mn1arcf, mnarcf, mnarcf1,
           comxmi, 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;}
   if( mnarcf1 != NULL ) {delete [] mnarcf1; mnarcf1=NULL;}
   if( mnarcf2 != NULL ) {delete [] mnarcf2; mnarcf2=NULL;}
-  if( mnarcf3 != NULL ) {delete [] mnarcf3; mnarcf3=NULL;}
 
   deltacpu_( d );
   tcpu += d;
@@ -606,7 +634,7 @@ void  aptrte( Z nutysu, R aretmx,
   // -----------------------------------------------------
   // boucle sur les triangles occupes (internes et externes)
   if( nust != NULL ) delete [] nust;
-  nust = new Z[4*nbt];
+  nust = new Z[nbsttria*nbt];
   if( nust == NULL ) goto ERREUR;
   nbt = 0;
   for (i=1; i<=mxartr; i++)
@@ -622,13 +650,12 @@ void  aptrte( Z nutysu, R aretmx,
       nust[nbt++] = 0;
     }
   }
-  nbt /= 4;  //le nombre final de triangles de la surface
-  cout << "Nombre de sommets=" << nbst
-       << "  Nombre de triangles=" << nbt << endl;
-
+  nbt /= nbsttria;  //le nombre final de triangles de la surface
+  MESSAGE( "APTRTE: Fin de la triangulation plane avec "<<nbst<<" sommets et "
+          << nbt << " triangles=" << nbt);
   deltacpu_( d );
   tcpu += d;
-  MESSAGE( "Temps total de la triangulation=" << tcpu << " secondes" );
+  MESSAGE( "APTRTE: Temps total de la triangulation plane=" << tcpu << " secondes" );
 
   // destruction des tableaux auxiliaires
   // ------------------------------------
@@ -638,7 +665,6 @@ void  aptrte( Z nutysu, R aretmx,
   if( mnslig != NULL ) delete [] mnslig;
   if( mnsoar != NULL ) delete [] mnsoar;
   if( mnpxyd != NULL ) delete [] mnpxyd;
-  if( mndalf != NULL ) delete [] mndalf;
   if( mntree != NULL ) delete [] mntree;
   if( mnqueu != NULL ) delete [] mnqueu;
   if( mntrsu != NULL ) delete [] mntrsu;
@@ -660,7 +686,7 @@ void  aptrte( Z nutysu, R aretmx,
   }
   else
   {
-    MESSAGE( "Triangulation non realisee " << ierr );
+    MESSAGE( "APTRTE: Triangulation NON REALISEE  avec erreur=" << ierr );
     if( ierr == 0 ) ierr=1;
     goto NETTOYAGE;
   }
@@ -737,9 +763,9 @@ void qualitetrte( R3 *mnpxyd,
       {
        //un triangle d'aire negative de plus
        nbtrianeg++;
-       cout << "ATTENTION: le triangle " << nt << " de sommets:"
+       MESSAGE("ATTENTION: le triangle " << nt << " de sommets:"
             << nosotr[0] << " " << nosotr[1] << " " << nosotr[2]
-            << " a une aire " << d <<"<=0" << endl;
+            << " a une aire " << d <<"<=0");
       }
 
       //aire des triangles actuels
@@ -749,10 +775,10 @@ void qualitetrte( R3 *mnpxyd,
 
   //les affichages
   quamoy /= nbtria;
-  cout << "Qualite moyenne=" << quamoy
+  MESSAGE("Qualite moyenne=" << quamoy
        << "  Qualite minimale=" << quamin
-       << " des " << nbtria << " triangles de surface totale="
-       << aire << endl;
+       << " des " << nbtria << " triangles de surface plane totale="
+       << aire);
 
   if( nbtrianeg>0 )
     MESSAGE( "ATTENTION: nombre de triangles d'aire negative=" << nbtrianeg );