1 // MEFISTO2: a library to compute 2D triangulation from segmented boundaries
4 // Copyright (C) 2006 Laboratoire J.-L. Lions UPMC Paris
6 // This library is free software; you can redistribute it and/or
7 // modify it under the terms of the GNU Lesser General Public
8 // License as published by the Free Software Foundation; either
9 // version 2.1 of the License.
11 // This library is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 // Lesser General Public License for more details.
16 // You should have received a copy of the GNU Lesser General Public
17 // License along with this library; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
20 // See http://www.ann.jussieu.fr/~perronne or email Perronnet@ann.jussieu.fr
22 // File : aptrte.cxx le C++ de l'appel du trianguleur plan
24 // Author : Alain PERRONNET
25 // Date : 13 novembre 2006
29 #include "utilities.h"
41 areteideale()//( R3 xyz, R3 direction )
46 //calcul de la longueur ideale de l'arete au sommet xyz (z ici inactif)
47 //dans la direction donnee
48 //a ajuster pour chaque surface plane et selon l'entier notysu (voir plus bas)
51 static double cpunew, cpuold=0;
53 void tempscpu_( double & tempsec )
54 //Retourne le temps CPU utilise en secondes
56 tempsec = ( (double) clock() ) / CLOCKS_PER_SEC;
57 //MESSAGE( "temps cpu=" << tempsec );
61 void deltacpu_( R & dtcpu )
62 //Retourne le temps CPU utilise en secondes depuis le precedent appel
65 dtcpu = R( cpunew - cpuold );
67 //MESSAGE( "delta temps cpu=" << dtcpu );
72 void aptrte( Z nutysu, R aretmx,
73 Z nblf, Z * nudslf, R2 * uvslf,
75 Z & nbst, R2 * & uvst,
78 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
79 // but : appel de la triangulation par un arbre-4 recouvrant
80 // ----- de triangles equilateraux
81 // le contour du domaine plan est defini par des lignes fermees
82 // la premiere ligne etant l'enveloppe de toutes les autres
83 // la fonction areteideale(s,d) donne la taille d'arete
84 // au point s dans la direction (actuellement inactive) d
85 // des lors toute arete issue d'un sommet s devrait avoir une longueur
86 // comprise entre 0.65 areteideale_(s,d) et 1.3 areteideale_(s,d)
89 // Les tableaux uvslf et uvpti sont supposes ne pas avoir de sommets identiques!
90 // De meme, un sommet d'une ligne fermee ne peut appartenir a une autre ligne fermee
94 // nutysu : numero de traitement de areteideale_(s,d) selon le type de surface
95 // 0 pas d'emploi de la fonction areteideale_() et aretmx est active
96 // 1 il existe une fonction areteideale_(s,d)
97 // dont seules les 2 premieres composantes de uv sont actives
98 // ... autres options a definir ...
99 // aretmx : longueur maximale des aretes de la future triangulation
100 // nblf : nombre de lignes fermees de la surface
101 // nudslf : numero du dernier sommet de chacune des nblf lignes fermees
102 // nudslf(0)=0 pour permettre la difference sans test
103 // Attention le dernier sommet de chaque ligne est raccorde au premier
104 // tous les sommets et les points internes ont des coordonnees
105 // UV differentes <=> Pas de point double!
106 // uvslf : uv des nudslf(nblf) sommets des lignes fermees
107 // nbpti : nombre de points internes futurs sommets de la triangulation
108 // uvpti : uv des points internes futurs sommets de la triangulation
112 // nbst : nombre de sommets de la triangulation finale
113 // uvst : coordonnees uv des nbst sommets de la triangulation
114 // nbt : nombre de triangles de la triangulation finale
115 // nust : 4 numeros dans uvst des sommets des nbt triangles
116 // s1, s2, s3, 0: no dans uvst des 3 sommets et 0 car quadrangle!
117 // ierr : 0 si pas d'erreur
119 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
120 // auteur : Alain Perronnet Laboratoire J.-L. LIONS Paris UPMC mars 2006
121 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
123 Z nbsttria=4; //Attention: 4 sommets stockes par triangle
124 //no st1, st2, st3, 0 (non quadrangle)
127 R3 direction=R3(0,0,0); //direction pour areteideale() inactive ici!
128 Z nbarfr=nudslf[nblf]; //nombre total d'aretes des lignes fermees
129 Z mxtrou = Max( 1024, nblf ); //nombre maximal de trous dans la surface
132 Z *mnsoar=NULL, mosoar=7, mxsoar, n1soar; //le hachage des aretes
133 Z *mnartr=NULL, moartr=3, mxartr, n1artr; //le no des 3 aretes des triangles
134 Z *mntree=NULL, motree=9, mxtree; //L'arbre 4 de TE et nombre d'entiers par TE
135 Z *mnqueu=NULL, mxqueu;
137 Z *mnarcf=NULL, mxarcf;
146 R3 comxmi[2]; //coordonnees UV Min et Maximales
147 R aremin, aremax; //longueur minimale et maximale des aretes
148 R airemx; //aire maximale souhaitee d'un triangle
152 Z i, l, n, ns, ns0, ns1, ns2, nosotr[3], nt;
153 Z mxsomm, nbsomm, nbarpi, nbarli, ndtri0, mn;
157 // initialisation du temps cpu
161 // quelques reservations de tableaux pour faire les calculs
162 // ========================================================
163 // declaration du tableau des coordonnees des sommets de la frontiere
164 // puis des sommets internes ajoutes
165 // majoration empirique du nombre de sommets de la triangulation
167 mxsomm = Max( 20000, 64*nbpti+i*i );
168 MESSAGE( "APTRTE: Debut de la triangulation plane avec " );
169 MESSAGE( "nutysu=" << nutysu << " aretmx=" << aretmx
170 << " mxsomm=" << mxsomm );
171 MESSAGE( nbarfr << " sommets sur la frontiere et " << nbpti << " points internes");
174 //mnpxyd( 3, mxsomm ) les coordonnees UV des sommets et la taille d'arete aux sommets
175 if( mnpxyd!=NULL ) delete [] mnpxyd;
176 mnpxyd = new R3[mxsomm];
177 if( mnpxyd==NULL ) goto ERREUR;
179 // le tableau mnsoar des aretes des triangles
180 // 1: sommet 1 dans pxyd,
181 // 2: sommet 2 dans pxyd,
182 // 3: numero de 1 a nblf de la ligne qui supporte l'arete
183 // 4: numero dans mnartr du triangle 1 partageant cette arete,
184 // 5: numero dans mnartr du triangle 2 partageant cette arete,
185 // 6: chainage des aretes frontalieres ou internes ou
186 // des aretes simples des etoiles de triangles,
187 // 7: chainage du hachage des aretes
188 // nombre d'aretes = 3 ( nombre de sommets - 1 + nombre de trous )
189 // pour le hachage des aretes mxsoar doit etre > 3*mxsomm!
190 // h(ns1,ns2) = min( ns1, ns2 )
191 if( mnsoar!=NULL ) delete [] mnsoar;
192 mxsoar = 3 * ( mxsomm + mxtrou );
193 mnsoar = new Z[mosoar*mxsoar];
194 if( mnsoar==NULL ) goto ERREUR;
195 //initialiser le tableau mnsoar pour le hachage des aretes
196 insoar( mxsomm, mosoar, mxsoar, n1soar, mnsoar );
198 // mnarst( mxsomm ) numero mnsoar d'une arete pour chacun des sommets
199 if( mnarst!=NULL ) delete [] mnarst;
200 mnarst = new Z[1+mxsomm];
201 if( mnarst==NULL ) goto ERREUR;
205 // mnslig( mxsomm ) no de sommet dans sa ligne pour chaque sommet frontalier
206 // ou no du point si interne forc'e par l'utilisateur
207 // ou 0 si interne cree par le module
208 if( mnslig!=NULL ) delete [] mnslig;
209 mnslig = new Z[mxsomm];
210 if( mnslig==NULL ) goto ERREUR;
211 azeroi( mxsomm, mnslig );
213 // initialisation des aretes frontalieres de la triangulation future
214 // renumerotation des sommets des aretes des lignes pour la triangulation
215 // mise a l'echelle des coordonnees des sommets pour obtenir une
216 // meilleure precision lors des calculs + quelques verifications
217 // boucle sur les lignes fermees qui forment la frontiere
218 // ======================================================================
223 for (n=1; n<=nblf; n++)
225 //l'initialisation de la premiere arete de la ligne n dans la triangulation
226 //-------------------------------------------------------------------------
227 //le sommet ns0 est le numero de l'origine de la ligne
229 mnpxyd[ns0].x = uvslf[ns0].x;
230 mnpxyd[ns0].y = uvslf[ns0].y;
231 mnpxyd[ns0].z = areteideale();//( mnpxyd[ns0], direction );
232 // MESSAGE("Sommet " << ns0 << ": " << mnpxyd[ns0].x
233 // << " " << mnpxyd[ns0].y << " longueur arete=" << mnpxyd[ns0].z);
235 //carre de la longueur de l'arete 1 de la ligne fermee n
236 d = pow( uvslf[ns0+1].x - uvslf[ns0].x, 2 )
237 + pow( uvslf[ns0+1].y - uvslf[ns0].y, 2 ) ;
238 aremin = Min( aremin, d );
239 aremax = Max( aremax, d );
241 //le numero des 2 sommets (ns1,ns2) de la premiere arete de la ligne
242 //initialisation de la 1-ere arete ns1-ns1+1 de cette ligne fermee n
243 //le numero des 2 sommets ns1 ns2 de la 1-ere arete
244 //Attention: les numeros ns debutent a 1 (ils ont >0)
245 // les tableaux c++ demarrent a zero!
246 // les tableaux fortran demarrent ou l'on veut!
251 //le numero n de la ligne du sommet et son numero ns1 dans la ligne
252 mnslig[ns0-1] = 1000000 * n + ns1-nudslf[n-1];
253 fasoar( ns1, ns2, moins1, moins1, n,
254 mosoar, mxsoar, n1soar, mnsoar, mnarst,
256 //pas de test sur ierr car pas de saturation possible a ce niveau
258 //le pointeur dans le hachage sur la premiere arete de la ligne fermee n
261 //la nouvelle arete est la suivante de l'arete definie juste avant
263 mnsoar[mosoar * noar - mosoar + 5] = noar0;
265 //l'initialisation des aretes suivantes de la ligne dans la triangulation
266 //-----------------------------------------------------------------------
267 nbarli = nudslf[n] - nudslf[n-1]; //nombre d'aretes=sommets de la ligne n
268 for (i=2; i<=nbarli; i++)
270 ns1 = ns2; //le numero de l'arete et le numero du premier sommet de l'arete
272 //nbs+1 est le 2-eme sommet de l'arete i de la ligne fermee n
275 //le 2-eme sommet de la derniere arete est le premier sommet de la ligne
278 //l'arete precedente est dotee de sa suivante:celle cree ensuite
279 //les 2 coordonnees du sommet ns2 de la ligne
281 //debut ajout 5/10/2006 ................................................
282 nuds = Max( nuds, ns ); //le numero du dernier sommet traite
283 //fin ajout 5/10/2006 ................................................
284 mnpxyd[ns].x = uvslf[ns].x;
285 mnpxyd[ns].y = uvslf[ns].y;
286 mnpxyd[ns].z = areteideale();//( mnpxyd[ns], direction );
287 // MESSAGE("Sommet " << ns << ": " << mnpxyd[ns].x
288 // << " " << mnpxyd[ns].y << " longueur arete=" << mnpxyd[ns].z);
290 //carre de la longueur de l'arete
291 d = pow( uvslf[ns2-1].x - uvslf[ns1-1].x, 2)
292 + pow( uvslf[ns2-1].y - uvslf[ns1-1].y, 2);
293 aremin = Min( aremin, d );
294 aremax = Max( aremax, d );
296 //debut ajout du 5/10/2006 .............................................
297 //la longueur de l'arete ns1-ns2
299 //longueur arete = Min ( aretmx, aretes incidentes )
300 mnpxyd[ns ].z = Min( mnpxyd[ns ].z, d );
301 mnpxyd[ns2-1].z = Min( mnpxyd[ns2-1].z, d );
302 //fin ajout du 5/10/2006 ...............................................
304 //le numero n de la ligne du sommet et son numero ns1 dans la ligne
305 mnslig[ns] = 1000000 * n + ns1-nudslf[n-1];
307 //ajout de l'arete dans la liste
308 fasoar( ns1, ns2, moins1, moins1, n,
309 mosoar, mxsoar, n1soar, mnsoar,
310 mnarst, noar, ierr );
311 //pas de test sur ierr car pas de saturation possible a ce niveau
313 //chainage des aretes frontalieres en position 6 du tableau mnsoar
314 //la nouvelle arete est la suivante de l'arete definie juste avant
315 mnsoar[ mosoar * noar0 - mosoar + 5 ] = noar;
318 //attention: la derniere arete de la ligne fermee enveloppe
319 // devient en fait la premiere arete de cette ligne
320 // dans le chainage des aretes de la frontiere!
322 if( ierr != 0 ) goto ERREUR;
324 aremin = sqrt( aremin ); //longueur minimale d'une arete des lignes fermees
325 aremax = sqrt( aremax ); //longueur maximale d'une arete
327 //debut ajout 9/11/2006 ................................................
328 // devenu un commentaire aretmx = Min( aretmx, aremax ); //pour homogeneiser
330 // protection contre une arete max desiree trop grande ou trop petite
331 if( aretmx > aremax*2.05 ) aretmx = aremax;
333 // protection contre une arete max desiree trop petite
334 if( (aremax-aremin) > (aremin+aremax)*0.05 && aretmx < aremin*0.5 )
335 aretmx =(aremin+aremax*2)/3.0;
337 if( aretmx < aremin && aremin > 0 )
340 //sauvegarde pour la fonction areteideale_
341 aretemaxface_ = aretmx;
343 //aire maximale souhaitee des triangles
344 airemx = aretmx * aretmx * sqrt(3.0) / 2.0; //Aire triangle equilateral
346 for(i=0; i<=nuds; i++ )
347 mnpxyd[i].z = Min( mnpxyd[i].z, aretmx );
348 //MESSAGE("Numero du dernier sommet frontalier=" << nuds+1);
349 //fin ajout 9/11/2006 .................................................
352 MESSAGE("Sur le bord: arete min=" << aremin << " arete max=" << aremax );
353 MESSAGE("Triangulation: arete mx=" << aretmx
354 << " triangle aire mx=" << airemx );
356 //chainage des aretes frontalieres : la derniere arete frontaliere
357 mnsoar[ mosoar * noar - mosoar + 5 ] = 0;
359 //tous les sommets et aretes frontaliers sont numerotes de 1 a nbarfr
360 //reservation du tableau des numeros des 3 aretes de chaque triangle
361 //mnartr( moartr, mxartr )
362 //En nombre: Triangles = Aretes Internes + Aretes Frontalieres - Sommets + 1-Trous
363 // 3Triangles = 2 Aretes internes + Aretes frontalieres
364 // d'ou 3T/2 < AI + AF => T < 3T/2 - Sommets + 1-Trous
365 //nombre de triangles < 2 ( nombre de sommets - 1 + nombre de trous )
366 if( mnartr!=NULL ) delete [] mnartr;
367 mxartr = 2 * ( mxsomm + mxtrou );
368 mnartr = new Z[moartr*mxartr];
369 if( mnartr==NULL ) goto ERREUR;
371 //Ajout des points internes
372 ns1 = nudslf[ nblf ];
373 for (i=0; i<nbpti; i++)
375 //les 2 coordonnees du point i de sommet nbs
376 mnpxyd[ns1].x = uvpti[i].x;
377 mnpxyd[ns1].y = uvpti[i].y;
378 mnpxyd[ns1].z = areteideale();//( mnpxyd[ns1], direction );
379 //le numero i du point interne
384 //nombre de sommets de la frontiere et internes
387 // creation de l'arbre-4 des te (tableau letree)
388 // ajout dans les te des sommets des lignes et des points internes imposes
389 // =======================================================================
390 // premiere estimation de mxtree
393 NEWTREE: //en cas de saturation de l'un des tableaux, on boucle
394 MESSAGE( "Debut triangulation avec mxsomm=" << mxsomm );
395 if( mntree != NULL ) delete [] mntree;
397 mntree = new Z[motree*(1+mxtree)];
398 if( mntree==NULL ) goto ERREUR;
400 //initialisation du tableau letree et ajout dans letree des sommets 1 a nbsomm
401 teajte( mxsomm, nbsomm, mnpxyd, comxmi, aretmx, mxtree, mntree, ierr );
407 //saturation de letree => sa taille est augmentee et relance
410 MESSAGE( "Nouvelle valeur de mxtree=" << mxtree );
416 MESSAGE( "Temps de l'ajout arbre-4 des Triangles Equilateraux=" << d << " secondes" );
417 if( ierr != 0 ) goto ERREUR;
418 //ici le tableau mnpxyd contient les sommets des te et les points frontaliers et internes
420 // homogeneisation de l'arbre des te a un saut de taille au plus
421 // prise en compte des tailles d'aretes souhaitees autour des sommets initiaux
422 // ===========================================================================
423 // reservation de la queue pour parcourir les te de l'arbre
424 if( mnqueu != NULL ) delete [] mnqueu;
426 mnqueu = new Z[mxqueu];
427 if( mnqueu==NULL) goto ERREUR;
429 tehote( nutysu, nbarpi, mxsomm, nbsomm, mnpxyd,
431 mntree, mxqueu, mnqueu,
436 MESSAGE("Temps de l'adaptation et l'homogeneisation de l'arbre-4 des TE="
437 << d << " secondes");
440 //destruction du tableau auxiliaire et de l'arbre
445 MESSAGE( "Redemarrage avec la valeur de mxtree=" << mxtree );
453 // trianguler les triangles equilateraux feuilles a partir de leurs 3 sommets
454 // et des points de la frontiere, des points internes imposes interieurs
455 // ==========================================================================
456 tetrte( comxmi, aretmx, nbarpi, mxsomm, mnpxyd,
457 mxqueu, mnqueu, mntree, mosoar, mxsoar, n1soar, mnsoar,
458 moartr, mxartr, n1artr, mnartr, mnarst,
461 // destruction de la queue et de l'arbre devenus inutiles
462 delete [] mnqueu; mnqueu=NULL;
463 delete [] mntree; mntree=NULL;
468 MESSAGE( "Temps de la triangulation des TE=" << d << " secondes" );
470 // ierr =0 si pas d'erreur
471 // =1 si le tableau mnsoar est sature
472 // =2 si le tableau mnartr est sature
473 // =3 si aucun des triangles ne contient l'un des points internes
474 // =5 si saturation de la queue de parcours de l'arbre des te
475 if( ierr != 0 ) goto ERREUR;
477 //qualites de la triangulation actuelle
478 qualitetrte_( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
479 nbt, quamoy, quamin );
481 // boucle sur les aretes internes (non sur une ligne de la frontiere)
482 // avec echange des 2 diagonales afin de rendre la triangulation delaunay
483 // ======================================================================
484 // formation du chainage 6 des aretes internes a echanger eventuellement
485 aisoar( mosoar, mxsoar, mnsoar, na );
486 tedela( mnpxyd, mnarst,
487 mosoar, mxsoar, n1soar, mnsoar, na,
488 moartr, mxartr, n1artr, mnartr, n );
490 MESSAGE( "Nombre d'echanges des diagonales de 2 triangles=" << n );
493 MESSAGE("Temps de la triangulation Delaunay par echange des diagonales="
494 << d << " secondes");
496 //qualites de la triangulation actuelle
497 qualitetrte_( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
498 nbt, quamoy, quamin );
500 // detection des aretes frontalieres initiales perdues
501 // triangulation frontale pour les restaurer
502 // ===================================================
504 if( mn1arcf != NULL ) delete [] mn1arcf;
505 if( mnarcf != NULL ) delete [] mnarcf;
506 if( mnarcf1 != NULL ) delete [] mnarcf1;
507 if( mnarcf2 != NULL ) delete [] mnarcf2;
508 mn1arcf = new Z[1+mxarcf];
509 if( mn1arcf == NULL ) goto ERREUR;
510 mnarcf = new Z[3*mxarcf];
511 if( mnarcf == NULL ) goto ERREUR;
512 mnarcf1 = new Z[mxarcf];
513 if( mnarcf1 == NULL ) goto ERREUR;
514 mnarcf2 = new Z[mxarcf];
515 if( mnarcf2 == NULL ) goto ERREUR;
517 terefr( nbarpi, mnpxyd,
518 mosoar, mxsoar, n1soar, mnsoar,
519 moartr, mxartr, n1artr, mnartr, mnarst,
520 mxarcf, mn1arcf, mnarcf, mnarcf1, mnarcf2,
523 MESSAGE( "Restauration de " << n << " aretes perdues de la frontiere ierr=" << ierr );
526 MESSAGE("Temps de la recuperation des aretes perdues de la frontiere="
527 << d << " secondes");
529 if( ierr != 0 ) goto ERREUR;
531 //qualites de la triangulation actuelle
532 qualitetrte_( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
533 nbt, quamoy, quamin );
535 // fin de la triangulation avec respect des aretes initiales frontalieres
537 // suppression des triangles externes a la surface
538 // ===============================================
539 // recherche du dernier triangle utilise
540 mn = mxartr * moartr;
541 for ( ndtri0=mxartr; ndtri0<=1; ndtri0-- )
544 if( mnartr[mn] != 0 ) break;
547 if( mntrsu != NULL ) delete [] mntrsu;
548 mntrsu = new Z[ndtri0];
549 if( mntrsu == NULL ) goto ERREUR;
551 if( mnlftr != NULL ) delete [] mnlftr;
552 mnlftr = new Z[nblf];
553 if( mnlftr == NULL ) goto ERREUR;
555 for (n=0; n<nblf; n++) //numero de la ligne fermee de 1 a nblf
558 tesuex( nblf, mnlftr,
559 ndtri0, nbsomm, mnpxyd, mnslig,
560 mosoar, mxsoar, mnsoar,
561 moartr, mxartr, n1artr, mnartr, mnarst,
564 delete [] mnlftr; mnlftr=NULL;
565 delete [] mntrsu; mntrsu=NULL;
569 MESSAGE( "Temps de la suppression des triangles externes=" << d << "ierr=" << ierr );
570 if( ierr != 0 ) goto ERREUR;
572 //qualites de la triangulation actuelle
573 qualitetrte_( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
574 nbt, quamoy, quamin );
576 // amelioration de la qualite de la triangulation par
577 // barycentrage des sommets internes a la triangulation
578 // suppression des aretes trop longues ou trop courtes
579 // modification de la topologie des groupes de triangles
580 // mise en delaunay de la triangulation
581 // =====================================================
582 mnarcf3 = new Z[mxarcf];
583 if( mnarcf3 == NULL )
585 cout << "aptrte: MC saturee mnarcf3=" << mnarcf3 << endl;
588 teamqt_( nutysu, aretmx, airemx,
589 mnarst, mosoar, mxsoar, n1soar, mnsoar,
590 moartr, mxartr, n1artr, mnartr,
591 mxarcf, mnarcf2, mnarcf3,
592 mn1arcf, mnarcf, mnarcf1,
593 nbarpi, nbsomm, mxsomm, mnpxyd, mnslig,
595 if( mnarcf3 != NULL ) {delete [] mnarcf3; mnarcf3=NULL;}
596 if( mn1arcf != NULL ) {delete [] mn1arcf; mn1arcf=NULL;}
597 if( mnarcf != NULL ) {delete [] mnarcf; mnarcf =NULL;}
598 if( mnarcf1 != NULL ) {delete [] mnarcf1; mnarcf1=NULL;}
599 if( mnarcf2 != NULL ) {delete [] mnarcf2; mnarcf2=NULL;}
603 MESSAGE( "Temps de l'amelioration de la qualite de la triangulation=" << d );
604 if( ierr == -13 ) ierr=0; //6/10/2006 arret de l'amelioration apres boucle infinie dans caetoi
605 if( ierr != 0 ) goto ERREUR;
607 //qualites de la triangulation finale
608 qualitetrte_( mnpxyd, mosoar, mxsoar, mnsoar, moartr, mxartr, mnartr,
609 nbt, quamoy, quamin );
611 // renumerotation des sommets internes: mnarst(i)=numero final du sommet
612 // ===================================
613 for (i=0; i<=nbsomm; i++)
616 for (nt=1; nt<=mxartr; nt++)
618 if( mnartr[nt*moartr-moartr] != 0 )
620 //le numero des 3 sommets du triangle nt
621 nusotr( nt, mosoar, mnsoar, moartr, mnartr, nosotr );
622 //les 3 sommets du triangle sont actifs
623 mnarst[ nosotr[0] ] = 1;
624 mnarst[ nosotr[1] ] = 1;
625 mnarst[ nosotr[2] ] = 1;
629 for (i=1; i<=nbsomm; i++)
635 // generation du tableau uvst de la surface triangulee
636 // ---------------------------------------------------
637 if( uvst != NULL ) delete [] uvst;
639 if( uvst == NULL ) goto ERREUR;
642 for (i=0; i<nbsomm; i++ )
647 uvst[nbst].x = mnpxyd[i].x;
648 uvst[nbst].y = mnpxyd[i].y;
650 //si le sommet est un point ou appartient a une ligne
651 //ses coordonnees initiales sont restaurees
658 //retour aux coordonnees initiales dans uvslf
660 n = n - 1000000 * l + nudslf[l-1] - 1;
661 uvst[nbst].x = uvslf[n].x;
662 uvst[nbst].y = uvslf[n].y;
666 //point utilisateur n interne impose
667 //retour aux coordonnees initiales dans uvpti
668 uvst[nbst].x = uvpti[n-1].x;
669 uvst[nbst].y = uvpti[n-1].y;
676 // generation du tableau 'nsef' de la surface triangulee
677 // -----------------------------------------------------
678 // boucle sur les triangles occupes (internes et externes)
679 if( nust != NULL ) delete [] nust;
680 nust = new Z[nbsttria*nbt];
681 if( nust == NULL ) goto ERREUR;
683 for (i=1; i<=mxartr; i++)
685 //le triangle i de mnartr
686 if( mnartr[i*moartr-moartr] != 0 )
688 //le triangle i est interne => nosotr numero de ses 3 sommets
689 nusotr( i, mosoar, mnsoar, moartr, mnartr, nosotr );
690 nust[nbt++] = mnarst[ nosotr[0] ];
691 nust[nbt++] = mnarst[ nosotr[1] ];
692 nust[nbt++] = mnarst[ nosotr[2] ];
696 nbt /= nbsttria; //le nombre final de triangles de la surface
697 MESSAGE( "APTRTE: Fin de la triangulation plane avec "<<nbst<<" sommets et "
698 << nbt << " triangles" );
701 MESSAGE( "APTRTE: Temps total de la triangulation plane=" << tcpu << " secondes" );
703 // destruction des tableaux auxiliaires
704 // ------------------------------------
706 if( mnarst != NULL ) delete [] mnarst;
707 if( mnartr != NULL ) delete [] mnartr;
708 if( mnslig != NULL ) delete [] mnslig;
709 if( mnsoar != NULL ) delete [] mnsoar;
710 if( mnpxyd != NULL ) delete [] mnpxyd;
711 if( mntree != NULL ) delete [] mntree;
712 if( mnqueu != NULL ) delete [] mnqueu;
713 if( mntrsu != NULL ) delete [] mntrsu;
714 if( mnlftr != NULL ) delete [] mnlftr;
715 if( mn1arcf != NULL ) delete [] mn1arcf;
716 if( mnarcf != NULL ) delete [] mnarcf;
717 if( mnarcf1 != NULL ) delete [] mnarcf1;
718 if( mnarcf2 != NULL ) delete [] mnarcf2;
719 if( mnarcf3 != NULL ) delete [] mnarcf3;
723 if( ierr == 51 || ierr == 52 )
725 //saturation des sommets => redepart avec 2 fois plus de sommets
732 MESSAGE( "APTRTE: Triangulation NON REALISEE avec erreur=" << ierr );
733 if( ierr == 0 ) ierr=1;
739 void qualitetrte_( R3 *mnpxyd,
740 Z & mosoar, Z & mxsoar, Z *mnsoar,
741 Z & moartr, Z & mxartr, Z *mnartr,
742 Z & nbtria, R & quamoy, R & quamin )
743 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
744 // but : calculer la qualite moyenne et minimale de la triangulation
745 // ----- actuelle definie par les tableaux mnsoar et mnartr
748 // mnpxyd : tableau des coordonnees 2d des points
749 // par point : x y distance_souhaitee
750 // mosoar : nombre maximal d'entiers par arete et
751 // indice dans mnsoar de l'arete suivante dans le hachage
752 // mxsoar : nombre maximal d'aretes stockables dans le tableau mnsoar
753 // attention: mxsoar>3*mxsomm obligatoire!
754 // mnsoar : numero des 2 sommets , no ligne, 2 triangles de l'arete,
755 // chainage des aretes frontalieres, chainage du hachage des aretes
756 // hachage des aretes = mnsoar(1)+mnsoar(2)*2
757 // avec mxsoar>=3*mxsomm
758 // une arete i de mnsoar est vide <=> mnsoar(1,i)=0 et
759 // mnsoar(2,arete vide)=l'arete vide qui precede
760 // mnsoar(3,arete vide)=l'arete vide qui suit
761 // moartr : nombre maximal d'entiers par arete du tableau mnartr
762 // mxartr : nombre maximal de triangles declarables
763 // mnartr : les 3 aretes des triangles +-arete1, +-arete2, +-arete3
764 // arete1 = 0 si triangle vide => arete2 = triangle vide suivant
767 // nbtria : nombre de triangles internes au domaine
768 // quamoy : qualite moyenne des triangles actuels
769 // quamin : qualite minimale des triangles actuels
770 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
773 Z nosotr[3], mn, nbtrianeg, nt, ntqmin;
783 for ( nt=1; nt<=mxartr; nt++ )
788 //un triangle occupe de plus
791 //le numero des 3 sommets du triangle nt
792 nusotr( nt, mosoar, mnsoar, moartr, mnartr, nosotr );
794 //la qualite du triangle ns1 ns2 ns3
795 qutr2d( mnpxyd[nosotr[0]-1], mnpxyd[nosotr[1]-1], mnpxyd[nosotr[2]-1],
801 //la qualite minimale
802 if( qualite < quamin )
808 //aire signee du triangle nt
809 d = surtd2( mnpxyd[nosotr[0]-1], mnpxyd[nosotr[1]-1], mnpxyd[nosotr[2]-1] );
812 //un triangle d'aire negative de plus
814 MESSAGE("ATTENTION: le triangle " << nt << " de sommets:"
815 << nosotr[0] << " " << nosotr[1] << " " << nosotr[2]
816 << " a une aire " << d <<"<=0");
819 //aire des triangles actuels
826 MESSAGE("Qualite moyenne=" << quamoy
827 << " Qualite minimale=" << quamin
828 << " des " << nbtria << " triangles de surface plane totale="
833 //le numero des 3 sommets du triangle ntqmin de qualite minimale
834 nusotr_( ntqmin, mosoar, mnsoar, moartr, mnartr, nosotr );
835 MESSAGE("Triangle de qualite minimale "<<quamin<<" de sommets:"
836 <<nosotr[0]<<" "<<nosotr[1]<<" "<<nosotr[2]<<" ");
837 for (int i=0;i<3;i++)
838 MESSAGE("Sommet "<<nosotr[i]<<": x="<< mnpxyd[nosotr[i]-1].x
839 <<" y="<< mnpxyd[nosotr[i]-1].y);
843 MESSAGE( "ATTENTION: "<< nbtrianeg << " TRIANGLES d'AIRE NEGATIVE" );