]> SALOME platform Git repositories - modules/smesh.git/blob - src/MEFISTO2/aptrte.cxx
Salome HOME
Join modifications from branch OCC_debug_for_3_2_0b1
[modules/smesh.git] / src / MEFISTO2 / aptrte.cxx
1 //  MEFISTO2: a library to compute 2D triangulation from segmented boundaries
2 //
3 //  Copyright (C) 2006  Laboratoire J.-L. Lions UPMC Paris
4 //
5 //  This library is free software; you can redistribute it and/or
6 //  modify it under the terms of the GNU Lesser General Public
7 //  License as published by the Free Software Foundation; either
8 //  version 2.1 of the License.
9 //
10 //  This library is distributed in the hope that it will be useful,
11 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
12 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
13 //  Lesser General Public License for more details.
14 //
15 //  You should have received a copy of the GNU Lesser General Public
16 //  License along with this library; if not, write to the Free Software
17 //  Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
18 //
19 //  See http://www.ann.jussieu.fr/~perronne or email Perronnet@ann.jussieu.fr
20 //
21 //  File   : aptrte.cxx   le C++ de l'appel du trianguleur plan
22 //  Module : SMESH
23 //  Author : Alain PERRONNET
24 //  Date   : 16 mars 2006
25
26 #include "Rn.h"
27 #include "aptrte.h"
28 #include "utilities.h"
29
30 using namespace std;
31
32 extern "C"
33 {
34   R aretemaxface_;
35   R areteideale_( R3 xyz, R3 direction )
36   {
37     return aretemaxface_;
38   }
39 }
40 //calcul de la longueur ideale de l'arete au sommet xyz (z ici inactif)
41 //dans la direction donnee
42 //a ajuster pour chaque surface plane et selon l'entier notysu (voir plus bas)
43
44
45 static double cpunew, cpuold=0;
46
47 void tempscpu_( double & tempsec )
48 //Retourne le temps CPU utilise en secondes
49 {  
50   tempsec = ( (double) clock() ) / CLOCKS_PER_SEC;
51   //MESSAGE( "temps cpu=" << tempsec );
52 }
53
54
55 void deltacpu_( R & dtcpu )
56 //Retourne le temps CPU utilise en secondes depuis le precedent appel
57 {
58   tempscpu_( cpunew );
59   dtcpu  = R( cpunew - cpuold );
60   cpuold = cpunew;
61   //MESSAGE( "delta temps cpu=" << dtcpu );
62   return;
63 }
64
65
66 void  aptrte( Z   nutysu, R      aretmx,
67               Z   nblf,   Z  *   nudslf,  R2 * uvslf,
68               Z   nbpti,  R2 *   uvpti,
69               Z & nbst,   R2 * & uvst,
70               Z & nbt,    Z  * & nust,
71               Z & ierr )
72 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
73 // but : appel de la triangulation par un arbre-4 recouvrant
74 // ----- de triangles equilateraux
75 //       le contour du domaine plan est defini par des lignes fermees
76 //       la premiere ligne etant l'enveloppe de toutes les autres
77 //       la fonction areteideale(s,d) donne la taille d'arete
78 //       au point s dans la direction (actuellement inactive) d
79 //       des lors toute arete issue d'un sommet s devrait avoir une longueur
80 //       comprise entre 0.65 areteideale_(s,d) et 1.3 areteideale_(s,d)
81 //
82 //Attention:
83 //  Les tableaux uvslf et uvpti sont supposes ne pas avoir de sommets identiques!
84 //  De meme, un sommet d'une ligne fermee ne peut appartenir a une autre ligne fermee
85 //
86 // entrees:
87 // --------
88 // nutysu : numero de traitement de areteideale_(s,d) selon le type de surface
89 //          0 pas d'emploi de la fonction areteideale_() et aretmx est active
90 //          1 il existe une fonction areteideale_(s,d)
91 //            dont seules les 2 premieres composantes de uv sont actives
92 //          ... autres options a definir ...
93 // aretmx : longueur maximale des aretes de la future triangulation
94 // nblf   : nombre de lignes fermees de la surface
95 // nudslf : numero du dernier sommet de chacune des nblf lignes fermees
96 //          nudslf(0)=0 pour permettre la difference sans test
97 //          Attention le dernier sommet de chaque ligne est raccorde au premier
98 //          tous les sommets et les points internes ont des coordonnees
99 //          UV differentes <=> Pas de point double!
100 // uvslf  : uv des nudslf(nblf) sommets des lignes fermees
101 // nbpti  : nombre de points internes futurs sommets de la triangulation
102 // uvpti  : uv des points internes futurs sommets de la triangulation
103 //
104 // sorties:
105 // --------
106 // nbst   : nombre de sommets de la triangulation finale
107 // uvst   : coordonnees uv des nbst sommets de la triangulation
108 // nbt    : nombre de triangles de la triangulation finale
109 // nust   : 4 numeros dans uvst des sommets des nbt triangles
110 //          s1, s2, s3, 0: no dans uvst des 3 sommets et 0 car quadrangle!
111 // ierr   : 0 si pas d'erreur
112 //        > 0 sinon
113 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
114 // auteur : Alain Perronnet  Laboratoire J.-L. LIONS Paris UPMC mars 2006
115 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
116 {
117   Z  nbsttria=4; //Attention: 4 sommets stockes par triangle
118                  //no st1, st2, st3, 0 (non quadrangle)
119
120   R  d, tcpu=0;
121   R3 direction=R3(0,0,0);  //direction pour areteideale() inactive ici!
122   Z  nbarfr=nudslf[nblf];  //nombre total d'aretes des lignes fermees
123   Z  mxtrou = Max( 1024, nblf );  //nombre maximal de trous dans la surface
124
125   R3 *mnpxyd=NULL;
126   Z  *mnsoar=NULL, mosoar=7, mxsoar, n1soar; //le hachage des aretes
127   Z  *mnartr=NULL, moartr=3, mxartr, n1artr; //le no des 3 aretes des triangles
128   Z  *mntree=NULL, motree=9, mxtree; //L'arbre 4 de TE et nombre d'entiers par TE
129   Z  *mnqueu=NULL, mxqueu;
130   Z  *mn1arcf=NULL;
131   Z  *mnarcf=NULL, mxarcf;
132   Z  *mnarcf1=NULL;
133   Z  *mnarcf2=NULL;
134   Z  *mnarcf3=NULL;
135   Z  *mntrsu=NULL;
136   Z  *mnslig=NULL;
137   Z  *mnarst=NULL;
138   Z  *mnlftr=NULL;
139
140   R3 comxmi[2];            //coordonnees UV Min et Maximales
141   R  aremin, aremax;       //longueur minimale et maximale des aretes
142   R  quamoy, quamin;
143
144   Z  noar0, noar, na;
145   Z  i, l, n, ns, ns0, ns1, ns2, nosotr[3], nt;
146   Z  mxsomm, nbsomm, nbarpi, nbarli, ndtri0, mn;
147   Z  moins1=-1;
148   R  dist;
149
150   aretemaxface_ = aretmx;
151
152   // initialisation du temps cpu
153   deltacpu_( d );
154   ierr = 0;
155
156   // quelques reservations de tableaux pour faire les calculs
157   // ========================================================
158   // declaration du tableau des coordonnees des sommets de la frontiere
159   // puis des sommets internes ajoutes
160   // majoration empirique du nombre de sommets de la triangulation
161   i = 4*nbarfr/10;
162   mxsomm = Max( 20000, 64*nbpti+i*i );
163   MESSAGE( "APTRTE: Debut de la triangulation plane avec " );
164   MESSAGE( "nutysu=" << nutysu << "  aretmx=" << aretmx << "  mxsomm=" << mxsomm );
165   MESSAGE( nbarfr << " sommets sur la frontiere et " << nbpti << " points internes");
166
167  NEWDEPART:
168   //mnpxyd( 3, mxsomm ) les coordonnees UV des sommets et la taille d'arete aux sommets
169   if( mnpxyd!=NULL ) delete [] mnpxyd;
170   mnpxyd = new R3[mxsomm];
171   if( mnpxyd==NULL ) goto ERREUR;
172
173   // le tableau mnsoar des aretes des triangles
174   // 1: sommet 1 dans pxyd,
175   // 2: sommet 2 dans pxyd,
176   // 3: numero de 1 a nblf de la ligne qui supporte l'arete
177   // 4: numero dans mnartr du triangle 1 partageant cette arete,
178   // 5: numero dans mnartr du triangle 2 partageant cette arete,
179   // 6: chainage des aretes frontalieres ou internes ou
180   //    des aretes simples des etoiles de triangles,
181   // 7: chainage du hachage des aretes
182   // nombre d'aretes = 3 ( nombre de sommets - 1 + nombre de trous )
183   // pour le hachage des aretes mxsoar doit etre > 3*mxsomm!
184   // h(ns1,ns2) = min( ns1, ns2 )
185   if( mnsoar!=NULL ) delete [] mnsoar;
186   mxsoar = 3 * ( mxsomm + mxtrou );
187   mnsoar = new Z[mosoar*mxsoar];
188   if( mnsoar==NULL ) goto ERREUR;
189   //initialiser le tableau mnsoar pour le hachage des aretes
190   insoar_( mxsomm, mosoar, mxsoar, n1soar, mnsoar );
191
192   // mnarst( mxsomm ) numero mnsoar d'une arete pour chacun des sommets
193   if( mnarst!=NULL ) delete [] mnarst;
194   mnarst = new Z[1+mxsomm];
195   if( mnarst==NULL ) goto ERREUR;
196   n = 1+mxsomm;
197   azeroi_( n, mnarst );
198
199   // mnslig( mxsomm ) no de sommet dans sa ligne pour chaque sommet frontalier
200   //               ou no du point si interne forc'e par l'utilisateur
201   //               ou  0 si interne cree par le module
202   if( mnslig!=NULL ) delete [] mnslig;
203   mnslig = new Z[mxsomm];
204   if( mnslig==NULL ) goto ERREUR;
205   azeroi_( mxsomm, mnslig );
206
207   // initialisation des aretes frontalieres de la triangulation future
208   // renumerotation des sommets des aretes des lignes pour la triangulation
209   // mise a l'echelle des coordonnees des sommets pour obtenir une
210   // meilleure precision lors des calculs + quelques verifications
211   // boucle sur les lignes fermees qui forment la frontiere
212   // ======================================================================
213   noar = 0;
214   aremin = 1e100;
215   aremax = 0;
216
217   for (n=1; n<=nblf; n++)
218   {
219     //l'initialisation de la premiere arete de la ligne n dans la triangulation
220     //-------------------------------------------------------------------------
221     //le sommet ns0 est le numero de l'origine de la ligne
222     ns0 = nudslf[n-1];
223     mnpxyd[ns0].x = uvslf[ns0].x;
224     mnpxyd[ns0].y = uvslf[ns0].y;
225     mnpxyd[ns0].z = areteideale_( mnpxyd[ns0], direction );
226 //     MESSAGE("Sommet " << ns0 << ": " << mnpxyd[ns0].x
227 //       << " " << mnpxyd[ns0].y << " longueur arete=" << mnpxyd[ns0].z);
228
229     //carre de la longueur de l'arete 1 de la ligne fermee n
230     d = pow( uvslf[ns0+1].x - uvslf[ns0].x, 2 ) 
231       + pow( uvslf[ns0+1].y - uvslf[ns0].y, 2 ) ;
232     aremin = Min( aremin, d );
233     aremax = Max( aremax, d );
234
235     //le numero des 2 sommets (ns1,ns2) de la premiere arete de la ligne
236     //initialisation de la 1-ere arete ns1-ns1+1 de cette ligne fermee n
237     //le numero des 2 sommets ns1 ns2 de la 1-ere arete
238     //Attention: les numeros ns debutent a 1 (ils ont >0)
239     //           les tableaux c++ demarrent a zero!
240     //           les tableaux fortran demarrent ou l'on veut!
241     ns0++;
242     ns1 = ns0;
243     ns2 = ns1+1;
244
245      //le numero n de la ligne du sommet et son numero ns1 dans la ligne
246     mnslig[ns0-1] = 1000000 * n + ns1-nudslf[n-1];
247     fasoar_( ns1, ns2, moins1, moins1, n,
248              mosoar, mxsoar, n1soar, mnsoar, mnarst,
249              noar0,  ierr );
250     //pas de test sur ierr car pas de saturation possible a ce niveau
251
252     //le pointeur dans le hachage sur la premiere arete de la ligne fermee n
253     //mndalf[n] = noar0;
254
255     //la nouvelle arete est la suivante de l'arete definie juste avant
256     if( noar > 0 )
257       mnsoar[mosoar * noar - mosoar + 5] = noar0;
258
259     //l'initialisation des aretes suivantes de la ligne dans la triangulation
260     //-----------------------------------------------------------------------
261     nbarli = nudslf[n] - nudslf[n-1];  //nombre d'aretes=sommets de la ligne n
262     for (i=2; i<=nbarli; i++)
263     {
264       ns1 = ns2; //le numero de l'arete et le numero du premier sommet de l'arete
265       if( i < nbarli )
266         //nbs+1 est le 2-eme sommet de l'arete i de la ligne fermee n
267         ns2 = ns1+1;
268       else
269         //le 2-eme sommet de la derniere arete est le premier sommet de la ligne
270         ns2 = ns0;
271
272       //l'arete precedente est dotee de sa suivante:celle cree ensuite
273       //les 2 coordonnees du sommet ns2 de la ligne
274       ns = ns1 - 1;
275       mnpxyd[ns].x = uvslf[ns].x;
276       mnpxyd[ns].y = uvslf[ns].y;
277       mnpxyd[ns].z = areteideale_( mnpxyd[ns], direction );
278 //       MESSAGE("Sommet " << ns << ": " << mnpxyd[ns].x
279 //         << " " << mnpxyd[ns].y << " longueur arete=" << mnpxyd[ns].z);
280
281       //carre de la longueur de l'arete
282       d = pow( uvslf[ns2-1].x - uvslf[ns1-1].x, 2) 
283         + pow( uvslf[ns2-1].y - uvslf[ns1-1].y, 2);
284       aremin = Min( aremin, d );
285       aremax = Max( aremax, d );
286
287       //le numero n de la ligne du sommet et son numero ns1 dans la ligne
288       mnslig[ns] = 1000000 * n + ns1-nudslf[n-1];
289
290       //ajout de l'arete dans la liste
291       fasoar_( ns1, ns2, moins1, moins1, n,
292                mosoar, mxsoar, n1soar, mnsoar,
293                mnarst, noar, ierr );
294       //pas de test sur ierr car pas de saturation possible a ce niveau
295
296       //chainage des aretes frontalieres en position 6 du tableau mnsoar
297       //la nouvelle arete est la suivante de l'arete definie juste avant
298       mnsoar[ mosoar * noar0 - mosoar + 5 ] = noar;
299       noar0 = noar;
300    }
301     //attention: la derniere arete de la ligne fermee enveloppe
302     //           devient en fait la premiere arete de cette ligne
303     //           dans le chainage des aretes de la frontiere!
304   }
305   if( ierr != 0 ) goto ERREUR;
306
307   aremin = sqrt( aremin );  //longueur minimale d'une arete des lignes fermees
308   aremax = sqrt( aremax );  //longueur maximale d'une arete
309
310   aretmx = Min( aretmx, aremax );  //pour homogeneiser
311   MESSAGE("nutysu=" << nutysu << "  aretmx=" << aretmx 
312        << "  arete min=" << aremin << "  arete max=" << aremax);
313
314   //chainage des aretes frontalieres : la derniere arete frontaliere
315   mnsoar[ mosoar * noar - mosoar + 5 ] = 0;
316
317   //tous les sommets et aretes frontaliers sont numerotes de 1 a nbarfr
318   //reservation du tableau des numeros des 3 aretes de chaque triangle
319   //mnartr( moartr, mxartr )
320   //En nombre: Triangles = Aretes Internes + Aretes Frontalieres - Sommets + 1-Trous
321   //          3Triangles = 2 Aretes internes + Aretes frontalieres
322   //       d'ou 3T/2 < AI + AF => T < 3T/2  - Sommets + 1-Trous
323   //nombre de triangles < 2 ( nombre de sommets - 1 + nombre de trous )
324   if( mnartr!=NULL ) delete [] mnartr;
325   mxartr = 2 * ( mxsomm + mxtrou );
326   mnartr = new Z[moartr*mxartr];
327   if( mnartr==NULL ) goto ERREUR;
328
329   //Ajout des points internes
330   ns1 = nudslf[ nblf ];
331   for (i=0; i<nbpti; i++)
332   {
333     //les 2 coordonnees du point i de sommet nbs
334     mnpxyd[ns1].x = uvpti[i].x;
335     mnpxyd[ns1].y = uvpti[i].y;
336     mnpxyd[ns1].z = areteideale_( mnpxyd[ns1], direction );
337     //le numero i du point interne
338     mnslig[ns1] = i+1;
339     ns1++;
340   }
341
342   //nombre de sommets de la frontiere et internes
343   nbarpi = ns1;
344
345   // creation de l'arbre-4 des te (tableau letree)
346   // ajout dans les te des sommets des lignes et des points internes imposes
347   // =======================================================================
348   // premiere estimation de mxtree
349   mxtree = 2 * mxsomm;
350
351  NEWTREE:  //en cas de saturation de l'un des tableaux, on boucle
352   MESSAGE( "Debut triangulation avec mxsomm=" << mxsomm );
353   if( mntree != NULL ) delete [] mntree;
354   nbsomm = nbarpi;
355   mntree = new Z[motree*(1+mxtree)];
356   if( mntree==NULL ) goto ERREUR;
357
358   //initialisation du tableau letree et ajout dans letree des sommets 1 a nbsomm
359   teajte_( mxsomm, nbsomm, mnpxyd, comxmi, aretmx, mxtree, mntree, ierr );
360   comxmi[0].z=0;
361   comxmi[1].z=0;
362
363   if( ierr == 51 )
364   {
365     //saturation de letree => sa taille est augmentee et relance
366     mxtree = mxtree * 2;
367     ierr   = 0;
368     MESSAGE( "Nouvelle valeur de mxtree=" << mxtree );
369     goto NEWTREE;
370   }
371
372   deltacpu_( d );
373   tcpu += d;
374   MESSAGE( "Temps de l'ajout arbre-4 des Triangles Equilateraux=" << d << " secondes" );
375   if( ierr != 0 ) goto ERREUR;
376   //ici le tableau mnpxyd contient les sommets des te et les points frontaliers et internes
377
378   // homogeneisation de l'arbre des te a un saut de taille au plus
379   // prise en compte des tailles d'aretes souhaitees autour des sommets initiaux
380   // ===========================================================================
381   // reservation de la queue pour parcourir les te de l'arbre
382   if( mnqueu != NULL ) delete [] mnqueu;
383   mxqueu = mxtree;
384   mnqueu = new Z[mxqueu];
385   if( mnqueu==NULL) goto ERREUR;
386
387   tehote_( nutysu, nbarpi, mxsomm, nbsomm, mnpxyd,
388            comxmi, aretmx,
389            mntree, mxqueu, mnqueu,
390            ierr );
391
392   deltacpu_( d );
393   tcpu += d;
394   MESSAGE("Temps de l'adaptation et l'homogeneisation de l'arbre-4 des TE="
395        << d << " secondes");
396   if( ierr != 0 )
397   {
398     //destruction du tableau auxiliaire et de l'arbre
399     if( ierr == 51 )
400     {
401       //letree sature
402       mxtree = mxtree * 2;
403       MESSAGE( "Redemarrage avec la valeur de mxtree=" << mxtree );
404       ierr = 0;
405       goto NEWTREE;
406     }
407     else
408       goto ERREUR;
409   }
410
411   // trianguler les triangles equilateraux feuilles a partir de leurs 3 sommets
412   // et des points de la frontiere, des points internes imposes interieurs
413   // ==========================================================================
414   tetrte_( comxmi, aretmx, nbarpi, mxsomm, mnpxyd,
415            mxqueu, mnqueu, mntree, mosoar, mxsoar, n1soar, mnsoar,
416            moartr, mxartr, n1artr, mnartr, mnarst,
417            ierr );
418
419   // destruction de la queue et de l'arbre devenus inutiles
420   delete [] mnqueu;  mnqueu=NULL;
421   delete [] mntree;  mntree=NULL;
422
423   //Temps calcul
424   deltacpu_( d );
425   tcpu += d;
426   MESSAGE( "Temps de la triangulation des TE=" << d << " secondes" );
427
428   // ierr =0 si pas d'erreur
429   //      =1 si le tableau mnsoar est sature
430   //      =2 si le tableau mnartr est sature
431   //      =3 si aucun des triangles ne contient l'un des points internes
432   //      =5 si saturation de la queue de parcours de l'arbre des te
433   if( ierr != 0 ) goto ERREUR;
434
435   //qualites de la triangulation actuelle
436   qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
437                nbt, quamoy, quamin );
438
439   // boucle sur les aretes internes (non sur une ligne de la frontiere)
440   // avec echange des 2 diagonales afin de rendre la triangulation delaunay
441   // ======================================================================
442   // formation du chainage 6 des aretes internes a echanger eventuellement
443   aisoar_( mosoar, mxsoar, mnsoar, na );
444   tedela_( mnpxyd, mnarst,
445            mosoar, mxsoar, n1soar, mnsoar, na,
446            moartr, mxartr, n1artr, mnartr, n );
447
448   MESSAGE( "Nombre d'echanges des diagonales de 2 triangles=" << n );
449   deltacpu_( d );
450   tcpu += d;
451   MESSAGE("Temps de la triangulation Delaunay par echange des diagonales="
452        << d << " secondes");
453
454   //qualites de la triangulation actuelle
455   qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
456                nbt, quamoy, quamin );
457
458   // detection des aretes frontalieres initiales perdues
459   // triangulation frontale pour les restaurer
460   // ===================================================
461   mxarcf = mxsomm/5;
462   if( mn1arcf != NULL ) delete [] mn1arcf;
463   if( mnarcf  != NULL ) delete [] mnarcf;
464   if( mnarcf1 != NULL ) delete [] mnarcf1;
465   if( mnarcf2 != NULL ) delete [] mnarcf2;
466   mn1arcf = new Z[1+mxarcf];
467   if( mn1arcf == NULL ) goto ERREUR;
468   mnarcf  = new Z[3*mxarcf];
469   if( mnarcf == NULL ) goto ERREUR;
470   mnarcf1 = new Z[mxarcf];
471   if( mnarcf1 == NULL ) goto ERREUR;
472   mnarcf2 = new Z[mxarcf];
473   if( mnarcf2 == NULL ) goto ERREUR;
474
475   terefr_( nbarpi, mnpxyd,
476            mosoar, mxsoar, n1soar, mnsoar,
477            moartr, n1artr, mnartr, mnarst,
478            mxarcf, mn1arcf, mnarcf, mnarcf1, mnarcf2,
479            n, ierr );
480
481   MESSAGE( "Restauration de " << n << " aretes perdues de la frontiere" );
482   deltacpu_( d );
483   tcpu += d;
484   MESSAGE("Temps de la recuperation des aretes perdues de la frontiere="
485        << d << " secondes");
486
487   if( ierr != 0 ) goto ERREUR;
488
489   //qualites de la triangulation actuelle
490   qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
491                nbt, quamoy, quamin );
492
493   // fin de la triangulation avec respect des aretes initiales frontalieres
494
495   // suppression des triangles externes a la surface
496   // ===============================================
497   // recherche du dernier triangle utilise
498   mn = mxartr * moartr;
499   for ( ndtri0=mxartr; ndtri0<=1; ndtri0-- )
500   {
501     mn -= moartr;
502     if( mnartr[mn] != 0 ) break;
503   }
504
505   if( mntrsu != NULL ) delete [] mntrsu;
506   mntrsu = new Z[ndtri0];
507   if( mntrsu == NULL ) goto ERREUR;
508
509   if( mnlftr != NULL ) delete [] mnlftr;
510   mnlftr = new Z[nblf];
511   if( mnlftr == NULL ) goto ERREUR;
512
513   for (n=0; n<nblf; n++)  //numero de la ligne fermee de 1 a nblf
514     mnlftr[n] = n+1;
515
516   tesuex_( nblf,   mnlftr,
517            ndtri0, nbsomm, mnpxyd, mnslig,
518            mosoar, mxsoar, mnsoar,
519            moartr, mxartr, n1artr, mnartr, mnarst,
520            nbt, mntrsu, ierr );
521
522   delete [] mnlftr; mnlftr=NULL;
523   delete [] mntrsu; mntrsu=NULL;
524
525   deltacpu_( d );
526   tcpu += d;
527   MESSAGE( "Temps de la suppression des triangles externes=" << d );
528   if( ierr != 0 ) goto ERREUR;
529
530   //qualites de la triangulation actuelle
531   qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
532                nbt, quamoy, quamin );
533
534   // amelioration de la qualite de la triangulation par
535   // barycentrage des sommets internes a la triangulation
536   // suppression des aretes trop longues ou trop courtes
537   // modification de la topologie des groupes de triangles
538   // mise en delaunay de la triangulation
539   // =====================================================
540   mnarcf3 = new Z[mxarcf];
541   if( mnarcf3 == NULL )
542   {
543     cout << "aptrte: MC saturee mnarcf3=" << mnarcf3 << endl;
544     goto ERREUR;
545   }
546   teamqt_( nutysu,
547            mnarst, mosoar, mxsoar, n1soar, mnsoar,
548            moartr, mxartr, n1artr, mnartr,
549            mxarcf, mnarcf2, mnarcf3,
550            mn1arcf, mnarcf, mnarcf1,
551            comxmi, nbarpi, nbsomm, mxsomm, mnpxyd, mnslig,
552            ierr );
553   if( mnarcf3 != NULL ) {delete [] mnarcf3; mnarcf3=NULL;}
554   if( mn1arcf != NULL ) {delete [] mn1arcf; mn1arcf=NULL;}
555   if( mnarcf  != NULL ) {delete [] mnarcf;  mnarcf =NULL;}
556   if( mnarcf1 != NULL ) {delete [] mnarcf1; mnarcf1=NULL;}
557   if( mnarcf2 != NULL ) {delete [] mnarcf2; mnarcf2=NULL;}
558
559   deltacpu_( d );
560   tcpu += d;
561   MESSAGE( "Temps de l'amelioration de la qualite de la triangulation=" << d );
562   if( ierr != 0 ) goto ERREUR;
563
564   //qualites de la triangulation finale
565   qualitetrte( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
566                nbt, quamoy, quamin );
567
568   // renumerotation des sommets internes: mnarst(i)=numero final du sommet
569   // ===================================
570   for (i=0; i<=nbsomm; i++)
571     mnarst[i] = 0;
572
573   for (nt=1; nt<=mxartr; nt++)
574   {
575     if( mnartr[nt*moartr-moartr] != 0 )
576     {
577       //le numero des 3 sommets du triangle nt
578       nusotr_( nt, mosoar, mnsoar, moartr, mnartr, nosotr );
579       //les 3 sommets du triangle sont actifs
580       mnarst[ nosotr[0] ] = 1;
581       mnarst[ nosotr[1] ] = 1;
582       mnarst[ nosotr[2] ] = 1;
583     }
584   }
585   nbst = 0;
586   for (i=1; i<=nbsomm; i++)
587   {
588     if( mnarst[i] >0 )
589       mnarst[i] = ++nbst;
590   }
591
592   // generation du tableau uvst de la surface triangulee
593   // ---------------------------------------------------
594   if( uvst != NULL ) delete [] uvst;
595   uvst = new R2[nbst];
596   if( uvst == NULL ) goto ERREUR;
597
598   nbst=-1;
599   for (i=0; i<nbsomm; i++ )
600   {
601     if( mnarst[i+1]>0 )
602     {
603       nbst++;
604       uvst[nbst].x = mnpxyd[i].x;
605       uvst[nbst].y = mnpxyd[i].y;
606
607       //si le sommet est un point ou appartient a une ligne
608       //ses coordonnees initiales sont restaurees
609       n = mnslig[i];
610       if( n > 0 )
611       {
612         if( n >= 1000000 )
613         {
614           //sommet d'une ligne
615           //retour aux coordonnees initiales dans uvslf
616           l = n / 1000000;
617           n = n - 1000000 * l + nudslf[l-1] - 1;
618           uvst[nbst].x = uvslf[n].x;
619           uvst[nbst].y = uvslf[n].y;
620         }
621         else
622         {
623           //point utilisateur n interne impose
624           //retour aux coordonnees initiales dans uvpti
625           uvst[nbst].x = uvpti[n-1].x;
626           uvst[nbst].y = uvpti[n-1].y;
627         }
628       }
629     }
630   }
631   nbst++;
632
633   // generation du tableau 'nsef' de la surface triangulee
634   // -----------------------------------------------------
635   // boucle sur les triangles occupes (internes et externes)
636   if( nust != NULL ) delete [] nust;
637   nust = new Z[nbsttria*nbt];
638   if( nust == NULL ) goto ERREUR;
639   nbt = 0;
640   for (i=1; i<=mxartr; i++)
641   {
642     //le triangle i de mnartr
643     if( mnartr[i*moartr-moartr] != 0 )
644     {
645       //le triangle i est interne => nosotr numero de ses 3 sommets
646       nusotr_( i, mosoar, mnsoar, moartr, mnartr,  nosotr );
647       nust[nbt++] = mnarst[ nosotr[0] ];
648       nust[nbt++] = mnarst[ nosotr[1] ];
649       nust[nbt++] = mnarst[ nosotr[2] ];
650       nust[nbt++] = 0;
651     }
652   }
653   nbt /= nbsttria;  //le nombre final de triangles de la surface
654   MESSAGE( "APTRTE: Fin de la triangulation plane avec "<<nbst<<" sommets et "
655            << nbt << " triangles=" << nbt);
656   deltacpu_( d );
657   tcpu += d;
658   MESSAGE( "APTRTE: Temps total de la triangulation plane=" << tcpu << " secondes" );
659
660   // destruction des tableaux auxiliaires
661   // ------------------------------------
662  NETTOYAGE:
663   if( mnarst != NULL ) delete [] mnarst;
664   if( mnartr != NULL ) delete [] mnartr;
665   if( mnslig != NULL ) delete [] mnslig;
666   if( mnsoar != NULL ) delete [] mnsoar;
667   if( mnpxyd != NULL ) delete [] mnpxyd;
668   if( mntree != NULL ) delete [] mntree;
669   if( mnqueu != NULL ) delete [] mnqueu;
670   if( mntrsu != NULL ) delete [] mntrsu;
671   if( mnlftr != NULL ) delete [] mnlftr;
672   if( mn1arcf != NULL ) delete [] mn1arcf;
673   if( mnarcf  != NULL ) delete [] mnarcf;
674   if( mnarcf1 != NULL ) delete [] mnarcf1;
675   if( mnarcf2 != NULL ) delete [] mnarcf2;
676   if( mnarcf3 != NULL ) delete [] mnarcf3;
677   return;
678
679  ERREUR:
680   if( ierr == 51 || ierr == 52 )
681   {
682     //saturation des sommets => redepart avec 2 fois plus de sommets
683     mxsomm = 2 * mxsomm;
684     ierr   = 0;
685     goto NEWDEPART;
686   }
687   else
688   {
689     MESSAGE( "APTRTE: Triangulation NON REALISEE  avec erreur=" << ierr );
690     if( ierr == 0 ) ierr=1;
691     goto NETTOYAGE;
692   }
693 }
694
695
696 void qualitetrte( R3 *mnpxyd,
697                   Z & mosoar, Z & mxsoar, Z *mnsoar,
698                   Z & moartr, Z & mxartr, Z *mnartr,
699                   Z & nbtria, R & quamoy, R & quamin )
700 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
701 // but :    calculer la qualite moyenne et minimale de la triangulation
702 // -----    actuelle definie par les tableaux mnsoar et mnartr
703 // entrees:
704 // --------
705 // mnpxyd : tableau des coordonnees 2d des points
706 //          par point : x  y  distance_souhaitee
707 // mosoar : nombre maximal d'entiers par arete et
708 //          indice dans mnsoar de l'arete suivante dans le hachage
709 // mxsoar : nombre maximal d'aretes stockables dans le tableau mnsoar
710 //          attention: mxsoar>3*mxsomm obligatoire!
711 // mnsoar : numero des 2 sommets , no ligne, 2 triangles de l'arete,
712 //          chainage des aretes frontalieres, chainage du hachage des aretes
713 //          hachage des aretes = mnsoar(1)+mnsoar(2)*2
714 //          avec mxsoar>=3*mxsomm
715 //          une arete i de mnsoar est vide <=> mnsoar(1,i)=0 et
716 //          mnsoar(2,arete vide)=l'arete vide qui precede
717 //          mnsoar(3,arete vide)=l'arete vide qui suit
718 // moartr : nombre maximal d'entiers par arete du tableau mnartr
719 // mxartr : nombre maximal de triangles declarables
720 // mnartr : les 3 aretes des triangles +-arete1, +-arete2, +-arete3
721 //          arete1 = 0 si triangle vide => arete2 = triangle vide suivant
722 // sorties:
723 // --------
724 // nbtria : nombre de triangles internes au domaine
725 // quamoy : qualite moyenne  des triangles actuels
726 // quamin : qualite minimale des triangles actuels
727 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
728 {
729   R  d, aire, qualite;
730   Z  nosotr[3], mn, nbtrianeg, nt;
731
732   aire   = 0;
733   quamoy = 0;
734   quamin = 2.0;
735   nbtria = 0;
736   nbtrianeg = 0;
737
738   mn = -moartr;
739   for ( nt=1; nt<=mxartr; nt++ )
740   {
741     mn += moartr;
742     if( mnartr[mn]!=0 )
743     {
744       //un triangle occupe de plus
745       nbtria++;
746
747       //le numero des 3 sommets du triangle nt
748       nusotr_( nt, mosoar, mnsoar, moartr, mnartr,  nosotr );
749
750       //la qualite du triangle ns1 ns2 ns3
751       qutr2d_( mnpxyd[nosotr[0]-1], mnpxyd[nosotr[1]-1], mnpxyd[nosotr[2]-1],
752                qualite );
753
754       //la qualite moyenne
755       quamoy += qualite;
756
757       //la qualite minimale
758       quamin = Min( quamin, qualite );
759
760       //aire signee du triangle nt
761       d = surtd2_( mnpxyd[nosotr[0]-1], mnpxyd[nosotr[1]-1], mnpxyd[nosotr[2]-1] );
762       if( d<0 )
763       {
764         //un triangle d'aire negative de plus
765         nbtrianeg++;
766         MESSAGE("ATTENTION: le triangle " << nt << " de sommets:"
767              << nosotr[0] << " " << nosotr[1] << " " << nosotr[2]
768              << " a une aire " << d <<"<=0");
769       }
770
771       //aire des triangles actuels
772       aire += Abs(d);
773     }
774   }
775
776   //les affichages
777   quamoy /= nbtria;
778   MESSAGE("Qualite moyenne=" << quamoy
779        << "  Qualite minimale=" << quamin
780        << " des " << nbtria << " triangles de surface plane totale="
781        << aire);
782
783   if( nbtrianeg>0 )
784     MESSAGE( "ATTENTION: nombre de triangles d'aire negative=" << nbtrianeg );
785   return;
786 }